© ESO 2019

1 Introduction

Understanding the first steps in the formation of protostars and protoplanetary disks is a great unsolved problem of modern astrophysics. Observationally, the key to constraining protostar formation models lies in high-resolution studies of the youngest protostars. Class 0 objects, which were originally discovered at millimeter wavelengths, are believed to be the youngest known accreting protostars (André et al. 1993). Because they are observed only t ≲ 4−9 × 10⁴ yr after the formation of a central hydrostatic protostellar object (Evans et al. 2009; Maury et al. 2011) while most of their mass is still in the form of a dense core/envelope (M_env > M_⋆), Class 0 protostars are believed to be representative of the main accretion phase of protostellar evolution and are likely to retain detailed information on the initial conditions of protostellar collapse (see review by André et al. 2000; Dunham et al. 2014).

High-resolution studies of Class 0 protostars are also key to constraining theoretical models for the formation of protostellar disks. At the simplest level, the formation of circumstellar disks is a natural consequence of the conservation of angular momentum during the collapse of rotating protostellar envelopes in the course of the main accretion phase (Cassen & Moosman 1981; Terebey et al. 1984). Hydrodynamic simulations show that in the absence of magnetic fields, rotationally supported disks form and quickly grow to reach large radii > 100 au after a few thousand years (Yorke & Bodenheimer 1999). These hydrodynamical disks are often massive enough to be gravitationally unstable, and their fragmentation has been suggested to contribute to the formation of brown dwarfs and multiple stellar systems (Stamatellos et al. 2007; Vorobyov 2010). On the other hand, early ideal magnetohydrodynamics (MHD) numerical simulations describing the protostellar collapse of magnetized envelopes had difficulties to form rotationally supported disks at scales r > 10 au because of strong magnetic braking: the increase in magnetic energy as field lines are dragged inward during protostellar collapse (Galli et al. 2006; Hennebelle & Ciardi 2009) reduces the envelope rotation and delays the formation of large disks. This so-called “magnetic braking catastrophe” was quickly mitigated, however, by including non-uniform initial conditions such as a magnetic field that is misaligned with the core rotation axis, transonic turbulent cores, or the treatment of radiative transfer (Joos et al. 2012; Seifried et al. 2012; Bate et al. 2014; Machida et al. 2014). More realistic numerical simulations including non-ideal MHD physics that allows dissipating and/or decoupling magnetic fields from the inner protostellar environment have recently been developed by several groups (Machida et al. 2011; Dapp et al. 2012; Masson et al. 2012; Li et al. 2014). Most numerical MHD studies now agree that including one or several resistive or dissipative effects (e.g., ambipolar diffusion, Ohmic dissipation, or the Hall effect) allows small (10 < r < 100 au) rotationally supported disks to be formed during magnetized protostellar collapse (Machida et al. 2014; Tsukamoto et al. 2015; Masson et al. 2016). The exact ingredients responsible for early disk properties during the main accretion phase remain widely debated, however.

From the observational point of view, young stellar objects (YSOs) have been studied in great detail in recent years, showing that large disks with radii ≳100 au are common in Class II (e.g. Andrews et al. 2009; Isella et al. 2009; Ricci et al. 2010; Spezzi et al. 2013) and Class I objects (Wolf et al. 2008; Jørgensen et al. 2009; Takakuwa et al. 2012; Eisner 2012; Lee et al. 2016; Sakai et al. 2016). However, we still lack good constraints on the properties of the progenitors of these disks during earlier phases of evolution: it has been difficult to observationally characterize the first stages of disk formation around Class 0 protostars because emission from the protostellar envelope dominates at most scales that are probed by single-dish telescopes (Motte & André 2001) or early interferometric observations (Looney et al. 2000). Long-wavelength observations are required to image deeply embedded disks and peer through dense protostellar envelopes. Subarcsecond resolution is needed to match the disk sizes at the typical distances of nearby star-forming clouds (d ~ 100–400 pc), as is high sensitivity to detect the weak fluxes of the youngest disks. Until recently, the small number of Class 0 protostars in nearby clouds and their relatively weak emission on small scales has restricted the millimeter interferometric studies that are required to reach subarcsecond (or ≲100 au) resolution to the most extreme objects. For example, a survey of bright Class 0 and Class I protostars with the Sub-Millimeter Array (SMA) by Jørgensen et al. (2009) attributed the detection of compact dust continuum emission components, all unresolved at the ~2′′ scales that these data probe, to the potential presence of disks with masses between 0.002 and 0.5 M_⊙ during the Class 0 phase. This simple interpretation was questionable since modeling of the millimeter continuum emission from Class 0 protostars sometimes indicates that an irregular density structure (e.g., complex envelope substructure or radial density enhancements at small scales) can lead to additional compact continuum emission in Class 0 protostars without the need of a disk component (Chiang et al. 2008; Maury et al. 2014).

Only the recent advent of powerful interferometric facilities with kilometer baselines at (sub)millimeter (submm) wavelengths has allowed the inner envelopes of Class 0 protostars to be explored at resolutions and sensitivities that are sufficient to distinguish envelope emission from resolved disk emission at the relevant scales (20–200 au). A pilot high-resolution study of five Class 0 protostars in Taurus and Perseus was carried out by Maury et al. (2010) with the IRAM Plateau de Bure Interferometer (PdBI). No large r > 100 au disks or protobinaries with separations 50 < a < 500 au were detected. Maury et al. (2010) concluded that the formation of protostellar disks and the fragmentation of dense cores into multiple systems at scales 50–500 au might be largely modified by magnetic fields during the main accretion phase. The apparent lack of large r > 100 au Class 0 disks could only be reproduced when magnetic-braking effects were included (Hennebelle & Teyssier 2008), while pure hydrodynamical models produced too many large r > 100 au disks. Subsequently, Maury et al. (2014) showed that any disk component in the NGC1333-IRAS2A protostar would need to be ≲ 40 au to reproduce the radial profile of the millimeter dust continuum emission observed with the PdBI. Maret et al. (2014) modeled the kinematics observed in methanol emission lines toward the same source and found that no significant rotation pattern was detected on similar scales. This confirmed the absence of a large disk in this Class 0 protostar.

The fast improvement of interferometric facilities such as ALMA has recently allowed a few other high-resolution studies to be made. A handful of Class 0 protostars have been proposed to have resolved disk-like rotation in their inner envelopes: while VLA 1623, HH212, and L1448-NB have been suggested to harbor Keplerian-like kinematics at scales 40 < r < 100 au (Murillo et al. 2013; Codella et al. 2014a; Tobin et al. 2016a), the most convincing case for a resolved protostellar disk was found in the L1527 Class 0/I protostar, where Ohashi et al. (2014) found a transition from a rotating infalling envelope to Keplerian motions in a disk at radii 50−60 au. On the other hand, Yen et al. (2015b) concluded that the disk in the Class 0 protostar B335 must have a radius r ≲ 10 au to reproduce the absence of a Keplerian pattern from the velocity field observed with ALMA at subarcsecond scale in this edge-on source.

To summarize, disks in young Class 0 protostars have largely remained elusive up to now. The intrinsic difficulty in distinguishing individual components in envelope-dominated objects has precluded any statistical observational constraints on the distribution of disks sizes and masses. Only a survey providing high angular resolution observations for a large sample of Class 0 protostars can shed light on the controversy about the pristine characteristics of protostellar disks and ultimately on the importance of magnetic fields in regulating disk formation during protostellar formation (Li et al. 2014).

The Continuum And Lines in Young ProtoStellar Objects (CALYPSO¹) IRAM Large Program is a survey of 16 Class 0 protostars, carried out with the IRAM Plateau de Bure (PdBI) interferometer in three spectral setups (around 94, 219, and 231 GHz). CALYPSO was crafted as an effort to make progress in our understanding of the angular momentum problem for star formation through high angular resolution observations (0.3′′, i.e., ≲ 100 au) of a significant sample of the youngest protostars. The main goals of the PdBI observations are to improve our understanding of (1) the formation of accretion disks and multiple systems during protostellar collapse, (2) the role of Class 0 jets and outflows in angular momentum extraction, and (3) the kinematics and structure of the inner protostellar environment. The 16 Class 0 objects of the CALYPSO sample are among the youngest known solar-type protostars (André et al. 2000), with M_env ~ 0.5–10 M_⊙, and internal luminosities L_int ~ 0.03–30 L_⊙ (see Table 1). This sample comprises most of the pre-Herschel confirmed Class 0 protostars located in nearby (d < 420 pc) clouds that could be observed in shared tracks from the IRAM-PdBI location. Several papers have been published based on the CALYPSO survey and discuss remarkable properties of individual protostars: IRAS2A (Maret et al. 2014; Codella et al. 2014b; Maury et al. 2014), IRAS4A (Santangelo et al. 2015), L1157 (Podio et al. 2016), or the Class I protostar SVS13A (Lefèvre et al. 2017). Two analyses of the molecular line emission in CALYPSO subsamples have also been carried out (Anderl et al. 2016; De Simone et al. 2017).

This paper is the first of a series of statistical analyses (Belloche et al., in prep.; Maret et al., in prep.; Podio et al., in prep.; Gaudel et al., in prep.) that are carried out for the whole CALYPSO sample. They address the three cornerstone questions described above that lie at the heart of the scientific motivations of CALYPSO. This paper focuses on the inner density structure(s) of the CALYPSO Class 0 protostars, analyzing the dual-frequency dust continuum emission visibilities to probe the structure of protostellar envelopes down to radii ~ 30 au (for the Taurus sources) to ~90 au (for the Serpens Main sources), with special emphasis on characterizing candidate protostellar disks and multiple systems. We show the dust continuum emission datasets, obtained at 1.37 and 3.18 mm, for the whole CALYPSO sample (see Table 1) in Sect. 2 and the dust continuum sources detected in our maps in Sect. 3. Section 4 is dedicated to describing the analysis of the dust continuum visibility dataset we performed to test whether candidate disk components are detected in each of the primary targets of our sample. In Sect. 5.1 we discuss the occurrence and properties of the candidate embedded disks in Class 0 protostellar envelopes from our sample and the literature. In Sect. 5.2 we explore current predictions from theoretical models of protostellar collapse, and the constraints our results bring to these models. Finally, we briefly discuss in Sect. 5.3 the multiplicity of Class 0 protostars in our sample and in the literature and compare them to the multiplicity properties of YSOs at more advanced evolutionary stages.

2 Observations and data reduction

2.1 PdBI observations and calibration

Observations of the CALYPSO sources were carried out with the IRAM PdBI between September 2010 and March 2013 (see details in Table A.1). We adopted a multiple configuration strategy, using the six-antenna array in the most extended configuration (A array) and in an intermediate antenna configuration (C array). This provided a fairly dense coverage of the uv-plane with 30 baselines ranging from 16 to 760 m. We used the WideX² correlator to cover a 3.8 GHz broadband spectral window for each of the three spectral setups, which were observed separately: at 1.29 mm for observations with a WideX central frequency of 231 GHz, at 1.37 mm for observations with a WideX central frequency of 219 GHz, and at 3.18 mm for observations with a central WideX frequency of 94 GHz. A higher spectral resolution correlator was placed onto a handful of molecular emission lines, but we present here only the analysis of the continuum emission extracted from the WideX dataset. The proximity of some of the sources on the sky allowed us to use common gain calibrators for several groups of sources and therefore to time-share a total of 37 tracks of ~8 h on the 16 sources. Each track was divided unequally among the sources, roughly in inverse proportion to the source luminosities so as to obtain the most homogeneous signal-to-noise ratios (S/N) in the sample. For each track, nearby amplitude and phase calibrators were observed to determine the time-dependent complex antenna gains. The calibration was performed in the GILDAS/CLIC³ environment. The correlator bandpass was calibrated on strong quasars (e.g., 3C273 and 3C454.3), while the absolute flux density scale was usually derived from observations of MWC349 and 3C84. The absolute flux calibration uncertainty is estimated to be ~10% at 94 GHz and ~15% at 231 GHz.

For sources where the peak of the dust continuum emission is detected at >40 σ (roughly equivalent to all sources with peak flux > 80 mJy beam⁻¹ at 231 GHz; see Table 2), we also carried out self-calibration of the continuum emission visibility dataset to improve the S/N at the longest baselines.

2.2 PdBIcontinuum emission maps

After a cross-check of the absolute flux calibration consistency of the observing dates, the visibility datasets obtained with the A and C configuration were combined, for each of the three frequency setups independently. The continuum visibilities were generated from the WideX units, avoiding channels in which line emission was detected at a level > 5 σ in the spectrum integrated over 2′′ around the peaks of the continuum emission reported in Table 3. Examples of wide-band spectra are shown in Maury et al. (2014), while the analysis of the molecular content of the whole sample will be presented in Belloche et al. (in prep.). In order to produce images and visibility curves with optimum sensitivity and best uv-coverage at 1.3 mm, we merged the continuum visibility data that were independently obtained at 219 and 231 GHz: we used the spectral index computed from the shortest common baseline of our PdBI observations at 94 and 231 GHz (α_20kλ; see Col. 6 of Table 4) for each individual source to scale the visibility amplitudes obtained at 219 GHz to 231 GHz. In the following, all 231 GHz fluxes, maps, and visibility profiles stem from the combined data.

Imaging of the continuum visibility tables was carried out using a robust scheme for weighting that allowed us to improve the resolution and lower the side lobes without exceedingly degrading the overall sensitivity. The robust (Briggs 1995) threshold we adopted was moderate (r = 0.3⁴) and allowed us to reach a good compromise between imaging quality and sensitivity while enhancing the contribution of the high spatial frequencies. This resulted in typical full widths at half-maximum (FWHMs) of the synthesized beams ≲ 0.5′′ at 231 GHz. Exceptions were made for the low-luminosity sources IRAM 04191, L1521F, and GF9-2, where natural weighting was used to maximize the sensitivity for these faint objects (producing synthesized beams ~ 1.0′′ at 231 GHz). The maps were subsequently cleaned using the Hogbom CLEAN algorithm provided in the GILDAS/MAPPING software, with a cleaning threshold set to twice the rms noise in the map. The properties of the resulting CLEANed maps (synthesizedHPBWs, rms noises) are reported in Table 2. The dust continuum emission maps obtained at 231 and 94 GHzare shown in Figs. 1–3 for two of the Perseus sources (L1448-2A and SVS13B) and a Taurus source (L1527), while the remaining dust continuum maps for the 13 other fields in the CALYPSO sample are shown in Figs. B.1–B.13.

3 Dust continuum sources in the CALYPSO maps

Dust continuum emission was successfully detected at 231 and 94 GHz toward all the Class 0 protostars we targeted in our survey: the positions and flux densities of the continuum sources detected in the fields observed with the PdBI are reported in Table 3. The target Class 0 protostars are labeled as “primary protostars” in the “source nature” column.

Multiple dust continuum emission components are sometimes detected in our PdBI maps. For each field that was observed with CALYPSO, we built a first sample of millimeter continuum sources from our 231 GHz PdBI dataset and selected components whose dust continuum emission was detected above the 5σ level in the maps (within the 21′′ FWHM primary beam). This detection threshold was chosen so as to limit contamination from spurious sources created by deconvolution of data sampling only partially the uv-plane. We note that throughout this paper, we always report and use rms noises at phase center, but the maps we show and the fluxes we report are corrected for the primary beam attenuation (the phase centers used in our observations are given in Table A.1). At the positions of the 231 GHz continuum sources, we checked the PdBI 94 GHz maps for independent detection of a counterpart at lower frequency. All the sources detected in the CALYPSO 231 GHz maps have a counterpart detected at > 3σ in the 94 GHz maps when it is possible to resolve them (the synthesized beam of the 94 GHz data is 2.5 times larger than the beam of the 231 GHz data).

Moreover, since the PdBI 94 GHz maps have a larger primary beam than the 231 GHz maps, four sources detected at > 11′′ from the main protostar in the 94 GHz map have no 231 GHz counterpart since they fall outside the primary beam of the higher frequency observations. In these cases (e.g., SVS13-C; see Fig. 2), the 94 GHz continuum source is reported only if the source was detected at another wavelength in the literature. Applying this method to all fields observed in the CALYPSO survey, we detected a total of 30 dust continuum emission sources. We report the coordinates and peak flux densities in Table 3 and discuss their nature in Sect. 5.3.

Fig. 1 1.3 mm (231 GHz, left panel) and 3.3 mm (94 GHz, right panel) PdBI dust continuum emission maps of L1448-2A, corrected for primary beam attenuation. The ellipses in the bottom left corner show the respective synthesized beam sizes. The contours are levels of −3σ (dashed), 5σ, and 10σ, then from 20σ in steps of 20σ (rms noise computed in the map before primary beam correction, reported in Table 2). The blue and red arrows show the direction of the protostellar jet(s) either as found in our CALYPSO dataset or from the literature (see Table 1).

Fig. 2 Same as Fig. 1 for the SVS13 region. The dashed pink contour shows the PdBI primary beam at each frequency. The contours show levels of −3σ (dashed), 5σ, and 10σ, then from 20σ in steps of 20σ (see Table 2). The blue and red arrows show the direction of the protostellar jets from our CALYPSO data for SVS13A/B (Podio et al., in prep.), while the SVS13C outflow PA stems from the CARMA map (Plunkett et al. 2013).

4 Analysis of dust continuum source structures

4.1 Building a proper visibility dataset for each primary target protostar

Each secondary source reported in Table 3 was removed from the visibility data by subtracting a Gaussian source model with FWHM at most twice the synthesized beam size (to remove only compact components and avoid affecting the envelope emission), or a point source if separated from the primary continuum component by less than two synthesized beams. During this subtraction process, the coordinates and peak fluxes of secondary components were fixed to the values reported in Cols. 2 and 3 of Table 3. Subtracting the visibilities that are due to secondary components allowed us to perform a focused analysis of the continuum emission that originates from each main protostar that was targeted by our observing program. We stress, however, that this operation does not preclude the possible presence of tighter multiple components at scales that cannot be resolved with PdBI. For example, the real parts of the visibilities of L1448-NB, IRAS2A, and L1157 show oscillations at long baselines in circularly averaged data, suggesting that additional structure or components (at a < 0.4′′) that are not centered on the phase center may remain in these three sources. From now on, we discuss and analyze the millimeter continuum emission originating solely from the primary protostar in each field.

For all primary protostars targeted with CALYPSO, we report the integrated flux densities (above a 5σ level) that we recovered with our PdBI observations in Table 4. Since the largest scale sampled by the PdBI observations depends on the frequency ν_obs and the shortest baseline in the array B_min (maximum recoverable scale MRS ~ 0.6c∕(ν_obs × B_min), i.e., 14′′ at 94 GHz and 6′′ at 231 GHz), we also report in Table 4 the matching visibility flux values, obtained at the shortest common B = 20 kλ baseline (value averaged in a 10 kλ bin). We note that the shortest baseline sampled at 94 GHz is ~8 kλ, so the integrated fluxes from the 94 GHz maps can sometimes be significantly higher than the 20 kλ visibility flux values at 94 GHz. For the continuum data at 231 GHz, and especially for strong sources, the flux recovered in the CLEANed maps is sometimes significantly lower than the flux measured at the 20 kλ baselines. This illustrates the difficulty of reliably reconstructing interferometric maps from data with incomplete uv-coverage, and it justifies our approach to carry out modeling in the uv-plane rather than in the image plane.

Fig. 3 Same as Fig. 1 for L1527. The contours show levels of −3σ, 5σ, and 10σ, and then from 20σ in steps of 40σ (see Table 2).

Fig. 4 Examples of radial intensity profiles from single-dish maps obtained at the IRAM-30 m telescope that are used as large-scale constraint for the envelope modeling. These profiles were obtained from the dust continuum emission maps acquired using the IRAM-30 m and are published in Motte & André (2001) and Kaas et al. (2004).

4.2 Large-scale constraints from single-dish observations of the dust continuum emission

We used a combination of PdBI configurations to observe the continuum emission down to baselines of 17 kλ at 231 GHz (8 kλ at 94 GHz), probing spatial scales up to ~6′′ at 231 GHz (14′′ at 94 GHz). Protostellar envelopes typically extend to larger angular scales at the distance of our sources (Motte & André 2001), and therefore cannot always be completely sampled using our PdBI observations alone (see Fig. 4 for an example of IRAM-30 m intensity profiles for four of the protostars in our sample). We therefore collected information from the literature on single-dish continuum observations of the target protostars that probe envelope scales ≳ 10′′, so as to better constrain the outer envelope density profiles through (i) dust continuum fluxes that are integrated at envelope scales reported in Col. 8 of Table 4 and (ii) peak flux densities obtained with single-dish telescopes that trace the material at the scales of the single-dish beam size, reported in Col. 10 of Table 4.

We can extrapolate the dust continuum fluxes found in the literature, which are often obtained at slightly different frequencies than our PdBI setups, when we assume a simple power-law dependence of the flux density on frequency by scaling the flux by an average spectral index at envelope scales. We used the spectral index computed from the shortest common baseline of our dual frequency PdBI observations at 20 kλ (see Col. 6 of Table 4) for each individual source to extrapolate the single-dish fluxes to the frequency of our PdBI observations. When secondary components were present in the field (see Table 3), we subtracted their flux densities as estimated in the PdBI u–v plane from the extrapolated single-dish flux of the primary sources. The resulting extrapolated total envelope fluxes and peak flux densities at 231 GHz are reported in Cols. 13 and 14 of Table 4, respectively, for the sources that are resolved by single-dish observations. Values for the total fluxes at 94 GHz were obtained by scaling the 231 GHz fluxes assuming the spectral index indicated in Col. 6.

Because single-dish continuum observations are broadband and the spectral index of the dust continuum emission in the envelope at larger scales is not well constrained by our PdBI dual-frequency observations (since the minimum common baseline at 94 and 231 GHz only probes scales up to ≲6′′), the fluxes at baselines <10 kλ extrapolated from single-dish data are uncertain. When we used them in our modeling of the envelope emission, we therefore allowed the total envelope flux to vary within a range ± 30%.

None of the four sources in NGC 1333 (IRAS2A, IRAS4A, IRAS4B, and SVS13B) is resolved by single-dish observations: the MAMBO bolometer-array studies by Motte & André (2001) and Chini et al. (1997) measured peak flux densities that are roughly equal to the integrated fluxes, integrated in areas twice the beam size. We therefore constrained the envelope fluxes of these sources to match the single-dish peak fluxes (± 40%) in an outer envelope radius smaller than or equal to the single-dish beam (± 40%). For L1521F, IRAM-30 m/MAMBO observations suggest a total flux 500–1000 mJy in a radius 30′′ at 243 GHz (see Motte & André 2001; Tokuda et al. 2016) and a peak flux 120 mJy in a 14′′ beam (similar to the value reported by Crapsi et al. 2004 in the nominal IRAM-30 m beam of 11′′). Although the L1521F envelope was resolved with MAMBO, the total flux is somewhat uncertain as the observations were affected by relativelystrong sky noise. We therefore only fixed the peak flux in our model fit, allowing the total integrated flux to vary by up to ± 50% around the single-dish value. For GF9-2 we used the IRAM-30 m fluxes from Wiesemeyer (1997): 60 mJy beam⁻¹ peak flux and 315 mJy integrated up to a radius 35′′ at 240 GHz, extrapolated to our frequencies. Using a spectral index α_{20 kλ} = 2.1 between 94 and 231 GHz (see Table 4), we expect an integrated flux of 291 mJy at 231 GHz (55 mJy peak flux) and 44 mJy at 94 GHz (8 mJy peak flux). Since these extrapolated envelope-scale fluxes are quite uncertain, we let them quite loose during the fitting process with an error bar of 50% at both frequencies.

Fig. 5 Top panels: 1.3 mm (231 GHz) PdBI dust continuum emission visibility real parts as a function of baseline length (circularly averaged in 20 kλ bins) for L1448-2A. The left panel uses a linear scale, while the right panel shows the same data and models in logarithmic scale to enhance the visibility of the long baseline points. The equivalent scales probed by the baselines in the image space are indicated at the top of each panel (computed as 0.6 × λ∕b). In both panels, the plain lines show the best-fit Plummer-only envelope (Pl, black) and Plummer + Gaussian (PG, red) models we found to reproduce the visibility profile of the dust continuum emission in this source. The dashedlines show the two components included in the best-fit PG model: the Plummer-envelope component in dashed lightblue and the additional Gaussian component in dot-dashed dark blue (see Table 6 for more information on the two models). The Plummer + Gaussian (PG) model is not statistically better than the Plummer-only (Pl) for L1448-2A. Bottom panels: same as the top panels for the 94 GHz visibility profiles.

4.3 Comparison to protostellar envelope models

From the continuum visibility dataset constructed for each primary protostar (see Sect. 4.1), we extracted the flux density of the protostellar dust continuum emission as a function of spatial frequency by averaging the real parts of the observed visibilities in bins of baseline lengths (bin widths of 20 kλ at 231 GHz, and bin widths of 10 kλ at 94 GHz). We obtained visibility curves sampling baselines [20–590] kλ at 231 GHz and [10–240] kλ at 94 GHz. Error bars on visibility amplitudes were estimated from the individual weights (w_i) of the interferometric visibility amplitudes (y_i), which were then averaged in bins of uv-distance: the error bar is the error on the mean value in an individual bin $(y_{\text{err}}=\Sigma (w_i^2\times {(y_i-w_{\text{mean}})}^2)/{(\Sigma w_i)}^2$ and w_mean = Σ(y_i × w_i)∕Σw_i). The continuum visibility curves for two of the Perseus sources (L1448-2A and SVS13B) and a Taurus source (L1527) are shown in Figs. 5–7, while the visibility curves of the other sources of the sample can be found in Figs. C.3–C.17. Our PdBI visibility profiles probe the radial distribution of the dust continuum emission from spatial scales θ = 0.35′′ to θ = 12′′. The envelope emission at larger scales is constrained by the single-dish dust continuum profiles, which are mostly obtained with the IRAM-30 m telescope (beam FWHM 11′′ at 1.3 mm).

We compared these dust continuum emission visibility curves to simple models of protostellar envelopes in a way similar to what was carried out in Maury et al. (2014) to analyze the preliminary CALYPSO data for the NGC 1333 IRAS2A protostar. The envelope models considered here and the results of their comparison with the dust continuum observations are described below.

4.3.1 Parameterized envelope models

For protostars in which a hydrostatic core has formed at the center, self-similar collapse solutions (Shu 1977; Whitworth & Summers 1985) suggest that the radial density profile of the envelope at large radii r > R_i is expected to range between ρ(r)∝ r^−3∕2 in the inner dynamical free-fall region (in spherical similarity solutions) and ρ(r) ∝ r⁻² in the outerregion where the initial conditions found in prestellar cores would have been conserved, where R_i can be the centrifugal radius if a disk develops or the radius of an inner cavity in the envelope. Since the radius of the expansion wavefront grows at the sound speed (~0.2 km s⁻¹ at 10 K), reaching ~2000 au in about 5 × 10⁴ yr (Class 0 lifetime), it seems reasonable to adopt a single power-law envelope model to describe the inner envelope probed by the PdBI observations (the largest scales probed are 500–2000 au). Since our sample consists of objects at different evolutionary stages from very early Class 0 to borderline Class 0/I sources (e.g., L1527) that are embedded in different environments, we chose to explore a range of density profiles from ρ(r)∝ r^−2.5 to the more shallow ρ(r) ∝ r⁻¹, in agreement with published millimeter continuum studies of resolved protostellar envelopes (Motte & André 2001; Looney et al. 2003; Shirley et al. 2002).

We adopted a simple spherically symmetric parameterized model with a Plummer-like density profile (Plummer 1911; Whitworth & Ward-Thompson 2001) to describe the protostellar density structure. In this model, the inner envelope has a nearly constant central density inside R_i, while the outer envelope hasa density profile approaching a power law r^−p, where p is a constant index: $$\rho (r)\hspace{0.17em}=\hspace{0.17em}{\rho }_{\text{i}}{\left[\frac{R_{\text{i}}}{{\left(R_{\text{i}}^2+r^2\right)}^{1/2}}\right]}^p\equiv \hspace{0.17em}\frac{{\rho }_{\text{i}}}{{\left(1+{\left(r/R_{\text{i}}\right)}^2\right)}^{p/2}}.$$(1)

The model has three free parameters: a truncation radius at R_out to account for the finite size of the envelope, the index p, which fixes the logarithmic slope of the density profile at radii R_out > r > R_i, and the inner radius R_i, which marks the end of the approximately uniform density of the inner region (0 < r ≤ R_i). In practice, R_i may be interpreted as the radius at which rotational motions start to dominate envelope infall motions.

The dust temperature in the inner region of a protostellar envelope is governed by infrared radiation from the release of gravitational energy of the material that accretes onto the central object. We ignore here the contribution of the interstellar radiation field, which mostly contributes to heating the outer layers of the envelope, although this contribution may be significant in the very low luminosity objects IRAM04191 and L1521F. We note that the dust continuum emission tracing the inner envelopes in our sample is optically thin at the two wavelengths we probed, as checked from the peak and integrated brightness temperatures from our PdBI maps (see Cols. 4 and 5 of Table 4), except for IRAS4A where the emission is partially optically thick at ~ 0.5′′ scales. Following Butner et al. (1990) and Terebey et al. (1993), the dust temperature distribution T(r) in the power-law envelope at radii r > R_i that is due to the central heating of the protostellar envelope may be approximated as follows in the case of optically thin dust continuum emission: $$T(r)=60{\left(\frac{r}{13\hspace{0.17em}400 \text{au}}\right)}^{-q}{\left(\frac{L_{\text{int}}}{{10}^5 L_{\odot }}\right)}^{q/2} \text{K},$$(2)

where r represents the radius from the central source of luminosity L_int. We varied the index q of the power-law temperaturedependence (T(r) ∝ r^−q) between q = 0.3 and q = 0.5 in our model fits.

The emerging intensity distribution of the envelope results from the combination of the dust density and temperature distributions of Eqs. (1) and (2). For instance, in the Rayleigh–Jeans regime, the radial intensity distribution of a power-law envelope is expected to scale as I(r)∝ r^−(p+q−1) (cf. Adams 1991). In practice, we therefore considered the following model intensity distribution: $$I(\overline{r})\hspace{0.17em}=\hspace{0.17em}\frac{I_0}{{\left(1+{\left(\overline{r}/R_{\text{i}}\right)}^2\right)}^{(p+q-1)/2}},$$(3)

where I₀ is the intensity from radii r < R_i and $\overline{r}$ represents the projected radius. Considering the combination of likely density and temperature distributions presented above, we let our model span a range from p + q = 1.3 to p + q = 2.9 between R_i and R_out.

Fig. 6 Same as Fig. 5 for SVS13B. The red Plummer + Gaussian (PG) model that includes a marginally resolved 0.2′′ FWHM Gaussian component is as good as the black Plummer-only (Pl) model in reproducing the 231 GHz visibility profile, but it provides a more reasonable value of the p + q index. Bottom panels: same as the top panels for the 94 GHz visibility profiles. The PG model does not perform any better than the Pl model for the 94 GHz visibility profile (see Table 6).

Fig. 7 Same as Fig. 5 for L1527. The red Plummer + Gaussian (PG) model that includes a 0.4′′ -FWHM Gaussian component (see Table 6 for more information on the two models) is statistically better than the black Plummer-only (Pl) model. Bottom panels: same as the top panels for the 94 GHz visibility profiles.

4.3.2 Fitting the dust continuum emission visibilities with envelope models

To model our interferometric observations, we converted the intensity distribution of the Plummer envelope model into a 1D visibility curve as a function of baseline $b=\sqrt{u^2+v^2}$ using a Hankel transform (corresponding to the 2D Fourier transform of a circularly symmetric function, see Bracewell 1965; Berger & Segransan 2007): $$V(b)=2\text{π}{\displaystyle {\int }_0^{\infty }{I}_{\text{ν}}}(r_{\text{b}})J_0\left(2\text{π}{r}_{\text{b}}b\right)r_{\text{b}}\text{d}{r}_{\text{b}},$$(4)

where J₀ is the zeroth-order Bessel function, $$J_0(z)=\frac{1}{2\text{π}}{\displaystyle {\int }_0^{2\text{π}}\exp }\left(-iz\cos \text{θ}\right)\text{dθ}.$$(5)

For example, this interferometric transform turns a spherically symmetric power-law intensity distribution I(r)∝ r^−(p+q−1) into a power-law visibility distribution as a function of u-v distance, V (b) ∝ b^p+q−3, solely determined by the power-law indices of the temperature and density profiles p and q in the envelope at r> R_i.

We performed a least-squares fit of the observed visibility data with the interferometric transform of the parameterized envelope model described by Eqs. (1) and (2). The envelope power-law index (p + q) and the inner cavity radius R_i were left as free parameters within the previously mentioned ranges, when we adjusted the observed profiles. The estimated parameter uncertainties were computed from the covariance matrix (see the description of the mpcurvefit in Markwardt 2009). The best-fit set of envelope parameters that reproduces the PdBI continuum visibility distribution for each of the 16 CALYPSO Class 0 sources is reported in Table 6. The modeled envelope visibility profiles are shown for three sources of our sample as a black curve overlaid on the observed visibilities in Figs. 5–7, while the models for other sources of the sample are shown in Figs. C.3–C.17.

The radial distributions of the millimeter dust continuum emission of 10 of the 16 Class 0 protostars in our sample (L1448-2A, IRAS2A, SVS13B, IRAS4A1, IRAM04191, L1521F, SerpM-S68N, SerpS-MM18, SerpS-MM22, and L1157) can be satisfactorily reproduced with models of Plummer-like protostellar envelopes with reduced χ² values ≤ 3 (see the model parameters and reduced χ² values reported in the first two lines of Table 6 for each source) in the two continuum bands that are probed in CALYPSO observations. Moreover, the model envelope sizes and power-law indices p + q for these 10 sources that we report in Table 6 show a good agreement in the two continuum bands (within 15%). This validates our choice of not fitting the two visibility profiles at the two frequencies jointly, and it shows that when a Plummer envelope model is satisfactory, there is no need to fix the envelope parameters to reproduce the dust continuum emission distribution at the longest baselines that are only probed by the 231 GHz data.

4.4 Comparison with envelope models including a disk-like component

For three protostars in our sample (L1448-C, L1448-NB1, and L1527; see Table 6), no single Plummer envelope model was able to satisfactorily reproduce the continuum emission data at both 231 and 94 GHz (reduced χ² > 3). At these two frequencies, the curvature of the radial intensity profile for these sources does not follow a simple power-law trend, or at uv-distances longer than 100 kλ, their PdBI visibility curve does not decrease quickly enough to be reproduced by one envelope model alone. In addition, while the 94 GHz visibility profiles of IRAS4B and GF92 are satisfactorily described by single-envelope models, their 231 GHz visibility profiles are inconsistent with single Plummer-like envelope models. This suggests that emission at baselines probed only by the 231 GHz data (250–550 kλ) differs from the best-fit envelope model we found to reproduce the 94 GHz continuum emission. Finally, while the 231 GHz profile can be modeled quite satisfactorily by a Plummer model (reduced χ² ≤ 3) for SerpM-SMM4, this does not hold for the 94 GHz profile.

The failure to reproduce one or both of the visibility profiles using circularly symmetric Plummer-like envelope models may reflect the presence in the source of either (i) asymmetric features in the envelope structure or (ii) an additional density or temperaturecomponent at small scales that is not properly modeled, such as a protostellar disk that is embedded in the envelope⁵. To check whether the visibility profiles of these six objects could include an additional continuum emission component originating from a disk that is embedded in the envelope, we also fit the visibility profile of each CALYPSO protostar with a two-component model consisting of a Plummer-like envelope (cf. Eq. (1)) and an additional compact circular Gaussian component located at phase center, whose flux F_Gauss and FWHM Θ_Gauss were left as free parameters (hereafter PG model for Plummer + Gaussian model). For these sources where the Gaussian component is spatially resolved, the value of R_i is enforced to be higher than the value of Θ_Gauss during the minimization process. The visibility function of the Gaussian disk component was calculated as $$V(b)\propto \exp (-\text{π}\times {\Theta }_{\text{Gauss}}\times b^2/4\text{ln}\hspace{0.17em}2).$$(6)

The parameters of the best-fit PG models are reported in Table 6. Figures C.3–C.17 show the Pl and PG models over the 231 and 94 GHz visibility curves of all the sources in our sample.

We performed a Fisher⁶ statistical test (later referred to as an F-test) to decide which of the best-fit Plummer-only model (Pl) and the best Plummer + Gauss (PG) models we obtained for each source in our sample was better. The goal was to avoid overmodeling the visibility profiles of sources with complex structure using a model that reproduced the observations better only coincidentally, because it includes two additional free parameters, such as the PG model. The results of the F-test (computing P values from the F distribution to test if the PG models are statistically better than the Pl models despite their two additional free parameters), which compared the best-fit Pl and PG models we obtained for each protostar are reported in Col. 10 of Table 6: we indicate a “yes” if the F-test suggests that the PG model is statistically better than the Pl model, and a “no” otherwise.

For all 10 sources that have satisfactory (reduced χ² < 3) Pl models at the two frequencies (231 and 94 GHz), none of the PG models provided a better fit for the two frequencies, except for SerpS-MM22. However, for five out of the six sources with unsatisfactory Pl models, the F-test suggests at either 231 GHz (L1448-C, L1527, IRAS4B, and GF92) or 94 GHz (L1448-C, L1527 and SerpM-SMM4) that the PG models reproduce the millimeter dust continuum visibility profiles significantly better than the Pl models.

For three protostars in our sample that are better reproduced with PG models at the two frequencies (L1527, SerpM-SMM4, and SerpS-MM22), our modeling suggests an additional Gaussian component is detected, that is well resolved in our 231 GHz PdBI data. While the inclusion of a Gaussian component significantly improved the minimization at both frequencies for two other protostars (L1448-C and GF9-2), the FWHMs of these components are very similar to the scale probed by our longest baseline (570 kλ, i.e., probing radii down to 0.26′′). This indicates that the additional Gaussian components suggested by our modeling are at best marginally resolved considering the increase in phase noise at the longest baselines, and that their sizes are therefore highly uncertain (see Table 6). For IRAS4B, the inclusion of a large and bright Gaussian component significantly improves the minimization for the 231 GHz visibility profile, which was unsatisfactorily modeled with a single Plummer envelope, but does not improve the modeling of our 94 GHz data. For five protostars that are well modeled at both frequencies using Pl models (SVS13B, IRAM04191, L1521F, SerpS-MM18, and L1157), the inclusion of unresolved Gaussian components in the PG models allows obtaining a significantly better minimization at 231 GHz, where we reach the best angular resolution, but does not improve the modeling of the 94 GHz visibility profiles. When we consider the null flux level dependence at the longest baselines in our PdBI observations, the decreasing slope of such low fluxes with increasing baseline at the longest baselines is questionable, and we therefore conservatively assume that these Gaussian components are detected but unresolved.

Finally, none of the models we attempted performed very well for L1448-NB1 (with reduced χ² values always higher than or equal to 3), probably because the emission at long baselines (> 200 kλ) shows wiggles that cannot be reproduced with a simple model such as we adopted here. In the specific case of this protostar, other models centered either on L1448-NB2 or in between the two millimeter sources L1448-NB1 and L1448-NB2 were attempted. They tentatively suggested the presence of an additional circumbinary structure (see Appendix C.2).

Figures 5–8 show examples of Plummer-only (Pl, black curve) and Plummer + Gauss (PG, red curve) model fits for four sources of our sample: L1448-2A, L1157, and SVS13B, which are satisfactorily described by an envelope-only model, and L1527, for which the inclusion of an additional central Gaussian source allows us to reproduce the visibility profile better than an envelope-only model.

4.5 Characterization of disk-like components in the sample

For a 2D Gaussian distribution, most of the radiation (90%) is emitted from within a radius r < FWHM, thus we conservatively defined the candidate disk radius to be the estimated FWHM size of the additional Gaussian component when the PG model provided a better model than the Pl model (as suggested by the F-test, see previous section and Col. 10 of Table 6).

For sources where the Pl model reproduces the visibility profile better at both frequencies (L1448-2A, IRAS4A, SerpM-S68N, and SerpS-MM18), we performed a new minimization of the PG model with fixed envelope parameters from the best-fit Pl model (PGf models in Appendix C). The Gaussian parameters of this PGf model were compared to the Gaussian parameters of the best-fit PG model (with free envelope parameters): the highest values provide upper limitsto the disk size and flux in the source. For sources for which the PG model reproduces the two the visibility profiles better (L1448-C, L1527, SerpM-SMM4, SerpS-MM22, and GF9-2), the parameters of the Gaussian component at 231 GHz (best angular resolution) are taken as the candidate disk size and disk flux in the source (or upper limitsif the Gaussian component in the PG model is unresolved). For sources for which only the 231 GHz visibility profile (probing the smallest spatial scales) is better reproduced by a PG model (SVS13B, IRAS4B, IRAM04191, L1521F, and L1157), we used the properties of the Gaussian component of the PG231 model (see Table 6) as candidate disk size and flux, or upper limits if the Gaussian component is unresolved. IRAS2A is the only source in the sample that is better reproduced by the PG model at 94 GHz, but not at 231 GHz. However, the Pl model at 94 GHz is also satisfactory (reduced chi square of 0.9) and the Gaussian component of the PG model is unresolved: hence, we consider it likely that an unresolved candidate disk is detected in our 94 GHz data and report the candidate disk as unresolved in Table 5, with the upper limits on its size and flux stemming from the best-fit PG model at 231 GHz (consistent with the parameters of the best-fit PG model at 94 GHz). Following this method, we report in Table 5 the detection of candidate protostellar disk-like structures in our CALYPSO data together with their radii and dust continuum fluxes at 231 GHz.

We are limited in the disk sizes we can probe (set by the spatial resolution of our data) and the surface brightness of disks we can detect (set by the sensitivity of our data at the longest baselines). In our PG modeling, the Gaussian size and flux are both free parameters: hence, if a large but faint disk component were to emerge above our sensitivity threshold, we would be able to model it with this description. Table 5 shows that when the visibility profile can be well modeled by a Pl model, the maximum flux that can be added in a disk-like component is generally low. Although it is not impossible that disks fainter than these upper limits would extend to larger radii than the upper-limit radii we report, it seems rather unlikely because (i) we would then detect their emission at our shorter baselines where our sensitivity is significantly better, and (ii) from a physical point of view, it seems unrealistic that < 1% of the envelope mass could reside in an extremely low-mass thin disk that is rotationally supported up to large radii at the center of an infalling envelope 100–200 times its mass. Thus, the maximum parameters of dust disk-like components provided by our simple models can be used robustly as upper limits on the sizes and fluxes of candidate disks.

We notethat the two Gaussian parameters (flux and FWHM) are somewhat degenerate when no Gaussian component is clearly detected (sensitivity and resolution wise): for the sources where the Pl model is best, we also performed additional PG models in which we tied the disk flux F_disk to 10% the envelope flux F_env. This ratio issimilar to those observed in the resolved disk candidates in our sample. Minimizing these PGt models to the visibility profiles (the disk size and envelope parameters were let free to vary) allowed us to set an upper limit to the radius of the disk-like component that can be added in each protostar that is well described by an envelope-only model. Because these upper-limit radii stem from a strong hypothesis (obtained when we fixed F_disk∕F_env = 0.1), we do not report them in Table 5. We do report their parameters (see PGt models) for each source individually in Appendix C, however. We stress that these PGt models are always worse than either the Pl or PG models that are performed with free parameters only, and the size of the Gaussian component is always found to be smaller than the best-fit PG model (where we did not constrain the ratio F_disk∕F_env). This exercise only points out that the protostellar disks in the sources that are well reproduced by the Pl model are in any case much less massive and much smaller than the resolved disk candidates found in the sources that are better reproduced with the PG model.

Moreover, we also carried out an analysis of the visibility profiles obtained in sectors of the uv-coverage to model only the dust continuum emission in the equatorial plane (at 90° from the outflow axis reported in Table 1) for all the strong sources in our sample. The parameters we found for this specific modeling (models Pleq and PGeq) are reported in Appendix C for each source with F_peak > 80 mJy beam⁻¹ at 231 GHz: the equatorial plane modeling produces results similar to the modeling using the full uv-coverage.

The choice of a Gaussian model to describe a possible disk contribution is also the simplest model that can be used to describe additional emission of unknown nature at small scales. We carried out tests using disks models with truncated power-law surface density profiles added to our Plummer envelope model for the two sources where extended or strong candidate-disk emission is detected (Serp-SMM4 and L1527), and we found very similar candidate disk properties as had been derived using Gaussian components (size and flux compatible within the error bars). Since such power-law disk models introduce two more free parameters in an already quite degenerate modeling, they are not well adapted to our PdBI data but may be used to perform 2D modeling in the uv-plane with future observations that will provide better resolution and sensitivity.

Fig. 8 Same as Fig. 5 for L1157. The black Plummer-only (Pl) model is statistically better than the red Plummer +Gaussian (PG) model that includes an unresolved Gaussian component (see Table 6 for more information on the two models). Bottom panels: same as the top panels for the 94 GHz visibility profiles.

5 Discussion

Here, we discuss the occurrence of large Class 0 disks from our analysis of the CALYPSO sample. We show that less than 25% of the CALYPSO Class 0 protostars include a resolved disk-like component with radii > 60 au. Taking into account all available literature on resolved Class 0 disks (or upper limits obtained with interferometers), we also show that a similar fraction (< 28%) of Class 0 protostars studied so far may harbor large disks with radii >60 au. We also discuss the multiplicity of the CALYPSO protostars and compare with the literature on Class 0, Class I, and Class II multiplicity. Finally, we argue that only magnetized protostellar collapse models can reproduce the disk size distribution that is found in theCALYPSO sample.

5.1 Occurrence of large Class 0 disks

5.1.1 In the CALYPSO sample

Our sample of Class 0 protostars was observed with sufficient resolution and sensitivity, at wavelengths that are sensitive to the bulk of dense circumstellar material, so that we probe the pristine properties of the progenitors of protoplanetary disks and multiple systems at the typical scales where they are observed at later stages of evolution. Based on the visibility analysis, we can detect continuum disk-like structures down to radii 0.15′′−0.2′′ at 231 GHz depending on the S/N at our longest baselines. We find the following (see Table 5): three low-mass Class 0 protostars have a well-resolved disk-like continuum structure at radii > 60 au (L1527, Serp-SMM4, and SerpS-MM22) that is detected at both 94 and 231 GHz. Two protostars (GF9-2 and L1448-C) are better described at these two frequencies by a model that includes a marginally resolved candidate disk component at scales 30–50 au. Six protostars have indications of a disk-like component that is only detected at 231 GHz (IRAS4B, SVS13B, IRAM04191, L1521F, SerpS-MM18, and L1157): of these six sources, only the disk candidate in IRAS4B is resolved by our 231 GHz observations. IRAS2A is marginally better described with an unresolved disk-like component at 94 GHz, and models with and without a Gaussian component are equally good at 231 GHz. In this case, we consider that an unresolved disk component might be present, although it is not formally detected. We report the maximum size (unresolved) and flux of the disk-like component from the Plummer + Gaussian model at 231 GHz. Finally, four protostars (L1448-2A, L1448-NB1, IRAS4A1, and SerpM-S68N) show no indication of a disk-like component at our sensitivity and spatial resolution.

In the case of the L1448-NB1/NB2 system, the individual protostars are not found to harbor individual disk-like structures at radii larger than ≳50 au. However, a structure of radius ~200−250 au is tentatively detected centered either on NB2 or in between the two millimeter sources (see Appendix C.2): this emission could trace the circumbinary structure reported in ALMA observations by Tobin et al. (2016a).

Furthermore, the case of IRAS4B is quite peculiar: if we do not constrain the total envelope flux to contain at least half the single-dish peak flux, then the PdBI visibilities of this source seem well described by a Gaussian component alone. Considering that this source is also unresolved by single-dish studies (and our modeling suggests a very compact envelope structure, if any), its Class 0 nature may be questioned and it may be more typical of the Class I stage (as, e.g., SVS13A) still embedded in its native cloud (resolved out by our interferometric observations, hence the Gaussian-like profile). It is also possible that the 231 GHz dust emission at PdBI scales is optically thick, but this is not suggested by either (i) the temperature brightness of the 231 GHz PdBI dust continuum emission peak or (ii) the similarly small envelope size found when modeling the optically thin 94 GHz dust continuum profile. Finally, we stress that the remarkably large size of the disk-like component in Serp-SMM4 (290 au) is quite intriguing and calls for further analysis of the kinematics at these scales to characterize the nature of the dust continuum emission, which does not follow a standard envelope power-law radial distribution.

When we assume that all resolved continuum structures we detected in our sample trace protostellar disks, the fraction of candidate protostellar disks with radii >60 au is 4 out of 16 (IRAS4B, L1527, SerpS-MM22, and SerpM-SMM4), i.e. ~ 25%. Moreover, our modeling shows that the candidate disk radii are <100 au in 14 out of 16 Class 0 protostars. Since these disk-like continuum structures need to be confirmed kinematically, the analysis of the dust continuum emission in the CALYPSO sample provides an upper limit to the occurrence of large (with radii ≳ 60 au) disks of ≤ 25% during the main accretion phase. Our results show that most Class 0 protostars have embedded disks (75% of the sources in our sample are better described when a disk component is included), but most of these young protostellar disks are small in size and flux.

5.1.2 Comparison with other works

Several groups have recently tried to characterize the properties of disks around Class 0 protostars. For example, Segura-Cox et al. (2016, 2018) used the VANDAM 8 mm continuum VLA survey of 43 Class 0 and Class 0/I protostars in Perseus. They found that only 15 of these Class 0 and 0/I sources are resolved at 15 au, and the candidate disks of only 10 of the VANDAM Class 0 protostars are resolved at radii > 12 au. All of these 10 disk Class 0 candidates have radii r < 45 au, and the disk component is not resolved in 67% of the VANDAM protostars at a 12 au scale. When the VANDAM and CALYPSO samples are combined, only 4 out of 52 Class 0 and 0/I protostars, that is, ≤ 8%, appear to have dust continuum circumstellar disks with radii ≳60 au. We did not include these disk sizes in our statistics for Class 0 protostars: as stated by the authors themselves, they may only represent lower limits because the 8 mm dust continuum emission size could be biased by a population of large dust grains that drift inward. We note that our CALYPSO results, which are obtained at shorter wavelengths, are less strongly subject to dust property limitations⁷, but we nevertheless also find that most Class 0 disks have small radii r < 60 au, with radius values consistent for the sources in common to both samples.

Other recent results for the quest of rotationally supported circumstellar disks in Class 0 protostars include L1527 (Ohashi et al. 2014, ALMA), which has a centrifugal radius of ~ 60 au, and VLA 1623 (Murillo et al. 2013), which was suggested to show Keplerian disk motions up to r ~ 189 au. The latter requires confirmation, however, because the model that was used is complex. In the HH212 protostar, Lee et al. (2018a) found a candidate Keplerian disk at radius ~ 40 au that is embedded in a dusty disk-like structure at radius ~ 60 AU. Recent SMA observations of the Class 0 protostar BHR7 suggest a compact dust continuum component at r ~ 120 au, although the kinematics obtained within the same dataset do not confirm Keplerian rotation at similar scales (Tobin et al. 2018). Gerin et al. (2017) analyzed ALMA 0.8 mm dust continuum observations of two Class 0 protostars, B1b-N and B1b-S. They found that compact disk-like components could only contribute to the intensity profiles at scales ≲50 au. The 0.87 mm ALMA observations by Tokuda et al. (2017) enabled the detection of a small 10 au radius disk in L1521F. For the B335 protostar, observations of the kinematics in the inner envelope with ALMA provided an upper limit to the disk radius of B335 < 15 au (Yen et al. 2015b). The two Class 0 protostars IRAS16253 and IRAS15398 have been observed with ALMA by Yen et al. (2017), who found a candidate disk radius of 20 au in IRAS15398 and an upper-limit disk radius of 6 au in IRAS16253. The Class 0 L1455 IRS1 was observed with the SMA, and the size of its protostellar disk was constrained to < 200 AU (Chou et al. 2016). The HH211 Class 0 protostar was observed with ALMA, and its disk radius was constrained to ≲ 10 au (Lee et al. 2018b). Finally, in an SMA survey of envelope kinematics toward a sample of 14 Class 0 protostars, Yen et al. (2015a) found that three sources (L1448-NB, Per-emb 9, and IRAS 03292+3039) harbored kinematics that were consistent with rotational motions dominating infalling motions in their inner envelopes at scales ~ 100−200 au. Two Class 0 protostars of their sample, L1527 and HH212, show rotationally dominated motions close to Keplerian at radii 50–80 au: rotationally supported disks with radii ~ 40−60 au have indeed be confirmed by other groups in these two sources. The remaining 8 Class 0 protostars of the Yen et al. (2015a) sample do not show any indication of rotationally dominated motions at scales sampled by the SMA observations (1.5′′–5′′ synthesized beams, and heterogeneous source distances). Our CALYPSO analysis triples the number of Class 0 protostars that are analyzed at high angular resolution to seek for disk-like structures at wavelengths that probe the bulk of circumstellar material (< 3 mm). When we combine the CALYPSO sample with the recent ALMA results that probe radii ≲60 au, only 5–7 out of 26 Class 0 protostars, that is, <28%, show either confirmed or candidate disks at radii >60 au. The disk radii as a function of the protostellar bolometric luminosity for all Class 0 disks that are characterized with millimeter-wavelength observations that probe scales ≳60 au from CALYPSO and the literature are shown as red symbols in Fig. 9.

The main caveat in characterizing candidate disk structures from continuum visibilities stems from the uncertain estimate of the large-scale envelope contribution. We stress that we left the envelope parameters free to vary in our modeling to account for the variety of envelopes that is observed in Class 0 protostars. We also emphasize that all of our best-fit models suggest envelope parameters that are physically reasonable. If the large-scale envelope deviates from the spherical symmetry assumed here, the estimated properties of the candidate disks may change slightly. However, the spherical envelope model has the advantage of simplicity and accurately reproduces the observed visibilities at short baselines. Our additional modeling, which used only the equatorial dust continuum profiles, suggests that the possible asymmetry of the envelope structure does not significantly affect our results.

We also acknowledge that the full extent of rotationally supported disks may be larger than the extent seen in millimeter dust continuum emission. We stress, however, that this is not the case in L1527, for example, where our estimate of the disk radius from dust continuum emission (56 ± 10 au) agrees with the radius that was determined from the molecular gas (centrifugal radius ~54–74 au; Ohashi et al. 2014; Aso et al. 2017) based on an analysis of the CO kinematics at an angular resolution similar to our CALYPSO data. We note, however, that Keplerian motions were tentatively detected up to r_c ~ 90 au in the SO emission of L1527 (CALYPSO data; Maret et al., in prep.).

Our study suggests that many Class 0 disks still remain uncharacterized because they lie at small scales that are not yet probed by current surveys in the (sub)millimeter regime, where the structure of Class 0 objects is best characterized. In the future, it will be of paramount importance to obtain accurate disk size and disk mass distributions at scales that are yet unexplored by most observational studies to make progress in our understanding of the formation and early evolution of stars and protoplanetary disks. Such studies will allow us to probe the disk kinematics and ultimately constrain the central protostellar mass and track its growth against the evolution of the envelope + disk system.

Fig. 9 Protostellar disk radii Rdisk as a functionof the protostellar bolometric luminosity Lbol for all Class 0 (red circles) and Class I (black stars) protostars observed between 0.7 and 2.7 mm with a spatial resolution better than 50 au from the CALYPSO sample and theliterature. The Class 0 disks shown here are from the CALYPSO sample (filled circles) and those characterized in the literature (contour circles, see text for names and references). When only upper limits on disk sizes could be estimated, the upper-limit values are indicated with arrows pointing down. The dotted red box shows the location of the ten Class 0 disk radii from the 8mm continuum observations of the VANDAM sample (Segura-Cox et al. 2018). The individual Class I disks shown are those described in the text of Sect. 5.1.3. The dashed lines show the median disk radii: in red from the sample of 25 Class 0 disks with $\widetilde{R_0}<$50 au (including upper limit radii), and in black the median radius from the sample of 25 Class I disks with $\widetilde{R_I}=115$ au.

5.1.3 Comparison with disk properties at later YSO stages

Early millimeter interferometric observations of T Tauri stars suggested that they are surrounded by rather large Keplerian protoplanetary disks with radii between 20 and 450 au (mean size of 165 au; Isella et al. 2009; Andrews et al. 2010; Guilloteau et al. 2011) and disk masses in the range 0.004–0.055 M_⊙ around 0.2–2.5 M_⊙ stars (Andrews et al. 2010). Surveys of the dust continuum emission from protoplanetary disks also suggest a trend that more massive T Tauri stars tend to host more massive disks (Andrews et al. 2013; Ribas et al. 2017). While the distribution of disk radii is not well characterized for Class I protostars, large (≥ 100 au) disks are generally found toward more evolved sources where M_⋆∕M_tot > 0.65 (where M_tot = M_⋆ + M_env + M_disk) (Barenfeld et al. 2017). Individual studies of Class I disks that are characterized by either dust or molecular lines observations between 0.7 mm and 2.7 mm with a spatial resolution better than 50 au include SVS13A, R-CrA, L1551NE, L1551-IRS1, L1455-IRS1, Lupus3-MMS, L1489, IRS43, IRS63, TMR1, TMC1, TMC1A, and L1536 (Lindberg et al. 2014; Takakuwa et al. 2014; Harsono et al. 2014; Aso et al. 2015; Yen et al. 2017), Elias29, WL12 (Miotello et al. 2014), and the 10 Taurus Class I from Sheehan & Eisner (2017). They are shown in Fig. 9 as black stars. We stress that some radii are estimated based on the kinematics of molecular lines (centrifugal radius), while some are determined from continuum-size measurements. Although there is no trend of systematic radius differences between the samples characterized by the two methods, we refer to the individual papers quoted for each source for more information on individual disk properties. The radii of these 25 Class I disks suggest that the median radius of Class I disks (115 au) is slightly larger than the median radius of Class 0 disks (< 50 au).

Recent observations of the dust continuum emission from large samples of T Tauri stars with new interferometric facilities (PdBI, ALMA) suggest that a population of faint and small (radii < 60 au) disks might be the bulk of circumstellar disks (Pascucci et al. 2016; Tripathi et al. 2017). This population of small dusty disks may have remained hidden so far because of earlier detectability limits that skewed the disk size distribution of Class II YSOs toward the largest disks. This is the case, for example, in the study by Tripathi et al. (2017), who measured the sizes of 50 Class II disks (Taurus and Ophiuchus star-forming regions): they found a size distribution ranging from 19 to 182 au, with a median of 48 au. Moreover, Tazzari et al. (2017) observed 22 Class II disks in Lupus, measuring effective radii ranging from 18 to 129 au, with a median of 55 au. Finally, in Chamaeleon I, Pascucci et al. (2016) showed that they were able to resolve only half of the disks in their sample (ALMA observations at 100 au resolution), with an indication that smaller dust disks are found around lower-mass stars.

The disk size distribution around Class 0 protostars remains to be sampled properly, but our results for Class 0 protostars suggest that most disks might start very small (<60 au, see Fig. 9). Although the viscous evolution of disks over a time span of a few Myr is highly unconstrained, if the pristine disk size distribution is indeed dominated by a large population of small disks, it would provide a natural explanation for the revised estimates of disk size distribution in more evolved YSOs. However, the dusty extent of T Tauri disks might also reflect the evolution of the dust during the pre-main-sequence phase (dust drift and coagulation), while their gas extent might have been affected by viscous evolution with time (Ansdell et al. 2018). We examine in the following possible scenarii to explain the lack of large Class 0 disks that is suggested by our CALYPSO observations.

5.2 Implications for theoretical models of protostellar disk formation

The collapseof rotating dense cores naturally produces rotationally supported disks as a result of angular momentum conservation (Terebey et al. 1984; Shu et al. 1987). For low-mass protostars, disks are expected to grow quickly in mass and size with time, reaching rather large size (r ~ 100 au; Cassen & Moosman 1981; Basu 1997; Matsumoto & Hanawa 2003; Bate 2018) on timescales of < 10⁴ yr (the centrifugal radius growswith t³ in the Shu model).

Star-forming clouds are permeated by large-scale magnetic fields, as recently pointed out by the analysis of the polarized dust emission obtained with the Planck spatial observatory (Planck Collaboration XIII 2016). Moreover, multiple results in the literature (Matthews et al. 2009; Girart et al. 2009; Hull et al. 2014) reported the detection of magnetic fields in protostellar cores. A recent SMA survey of low-mass protostars (Galametz et al. 2018) shows that magnetic fields are detected in all low-mass protostellar cores: early non-detections of the magnetization of cores probably were due to sensitivity limitations and the intrinsic difficulty of properly observing polarized dust emission.

Even modest values of the magnetic field strength significantly modify the collapse models because the radial and toroidal components of the magnetic field are amplified by the differential motions within the collapsing core. This is especially true in models with low turbulent energy, which are representative of low-mass cores: in the following, we discuss relevant magnetized models that include rotation (rotational over gravitational energy β ~ 1−5%) and low levels of turbulence (subsonic). In moderately magnetized models with a ratio of mass-to-flux over critical mass-to-flux ratio, $\mu =(M/\Phi )/{(M/\Phi )}_c<10$ (with ${(M/\Phi )}_c=c_1/3\text{ π}\times {(5/G)}^{1/2}$; see Mouschovias & Spitzer 1976), the collapse primarily occurs along the field lines, thus reducing the angular momentum that is transported into the inner envelope, whereas at the same time, strong magnetic braking occurs (Galli et al. 2006; Krasnopolsky et al. 2010). As a consequence, no centrifugally supported disks forms in these ideal MHD numerical models (Hennebelle & Fromang 2008; Mellon & Li 2008): only when the magnetic braking drops efficiency because the surrounding envelope dissipates do these models predict the rapid growth of disks at up to a few hundred au near the end of the main accretion phase (Machida et al. 2011). To mitigate this effect, magnetized numerical simulations of core collapse were carried out using different initial conditions in an attempt to weaken the magnetic braking effect, notably by introducing an angle between the core rotation axis and the magnetic field direction at core scale. These models showed that for a low mass-to-flux ratio (μ < 3), magnetic braking still completely inhibits the formation of rotationally supported disks, while for weakly magnetized cases μ > 5, the formation of small Keplerian disks would become possible if the difference in angle ⟨ j, B⟩ is large enough (typically >45°; see Joos et al. 2012; Krumholz et al. 2013; Ciardi & Hennebelle 2010). Under some specific conditions, models that include turbulence also form large disks because turbulence can induce some misalignment that reduces the magnetic braking efficiency (Santos-Lima et al. 2013; Seifried et al. 2015).

Recently, non-ideal MHD models were carried out. They included the effects of different mechanisms that allowed the field to diffuse under some conditions (for a review, see Li et al. 2014). Generally, the inclusion of diffusive physics in numerical models weakens the magnetic braking in the inner envelope and might allow the formation of rotationally supported disks, although these are systematically much smaller than disks that formed in purely hydrodynamical models. The exact disk properties from non-ideal MHD simulations of protostellar collapse are still heavily debated, however, see for example the diverging results regarding ambipolar diffusion: Mellon & Li (2009), Li et al. (2011), and Dapp et al. (2012) concluded that AD fails to enable the formation of rotationally supported disks, while Zhao et al. (2016) proposed that the depletion of small grains could increase the ambipolar diffusivity and hence allow the formation of such disks. A recent analytical model by Hennebelle et al. (2016), which has been compared with a series of numerical simulations, argued that AD in combinationwith magnetic braking during magnetized collapse leads to the formation of small disks, whose sizes (~ 20 au) are remarkably constant and only weakly depend on the initial rotational energy and magnetization of the cores (if β and μ are not too large). This study thus suggests that magnetic self-regulation during the collapse may produce a rather narrow disk size distribution during the protostellar phase.

At the current state of our knowledge, our PdBI results that suggest the paucity of large r ≳ 60 au disks duringthe Class 0 phase may either point to

• (i)

  the extreme youth of these objects (which may not yet have developed a disk),

• (ii)

  an overestimate of the initial angular momentum of the cores in hydrodynamical collapse models, or

• (iii)

  a magnetically regulated disk formation scenario.

The first possibility (i) is mostly ruled out by the presence of the large-scale outflows that are observed to be powered by every protostar in our sample (Podio et al., in prep.). This testifies that accretion onto the central protostar has been proceeding for at least a few thousand years, similar to the typical timescale at which hydrodynamical models of the collapse of rotating dense cores predict large (r ~ 100 au) rotationally supported disks as a consequence of angular momentum conservation (Goodwin et al. 2004).

The possibility that rotation is overestimated in hydrodynamical models of protostellar collapse (ii) can be addressed since rotation (expressed as the ratio of rotational over gravitational energy β) has been observed to be β ~ 2−10% in prestellar cores (Goodman et al. 1993; Caselli et al. 2002b; Belloche 2013). The typical angular momentum values observed in protostellar cores are on the order of a few 10⁻³ km s⁻¹.pc (Belloche 2013; Yen et al. 2015a). A possible systematic over-estimate of the angular momentum has been reported from the analysis of numericalsimulations of the collapse of protostellar cores (Dib et al. 2010), although this was not confirmed by a recent numerical study (Zhang et al. 2018). Moreover, the rotation that is observed is usually interpreted as the value for a uniform core ${\beta }_{\text{obs}}=\frac{1}{3}\times \frac{R^3{\Omega }^2}{GM}\times \text{sin}(i),$ with i the inclination angle, R and M the radius and mass of the core, respectively. Protostellar dense cores are centrally concentrated, however, which increases the gravitational energy and decreases the rotational energy, so that β could actually be smaller (a factor of ~3; see Belloche et al. 2002; Belloche 2013). Hydrodynamical collapse, even when the cores have only a small fraction of the observed β, quickly develops much larger disks (e.g., Goodwin et al. 2004) than are observed in our sample because of angular momentum conservation if no mechanism is found to dissipate or redistribute the envelope angular momentum. For example, for a centrally condensed (ρ ∝ r⁻²) protostellar dense core of diameter 0.1 pc that is rotating uniformly with β = 0.02, gravitational collapse would lead to the formation of a disk at r_c ∝ β × R_core = 200 au when half the envelope mass has collapsed into a protostar-disk system (i.e., the end of the main-accretion phase). To produce a disk size that is consistent with the bulk of our observations, we would have to remove > 75% of the angular momentum from the infalling material at scales r > r_c in order to shrink the radius of the disk to r_c ≲ 50 au. Although it is unclear whether we have robust quantitative constraints on the angular momentum enclosed in star-forming cores, an overestimate of β from observations by one order of magnitude, as would be requested to explain our observations, seems highly unlikely.

Magnetic torques acting in the envelope are, as far as we are aware, the only mechanism that is able to greatly reduce the transport of angular momentum during protostellar collapse (Li et al. 2014), as was shown by recent results that pointed toward a magnetically regulated disk formation scenario in the B335 protostar (Maury et al. 2018). Finally, if the magnetic field is a key player during the protostellar collapse that leads to the formation of solar-type stars, as suggested by our observations, we stress that identifying an embedded protostellar disk requires spectral line observations that show Keplerian rotation (such as L1527) before they are considered robust rotationally supported disks. Magnetized models of protostellar formation generically produce flattened inner envelope structures because the collapse is favored along the main direction of the magnetic field at core scale. Such disk-like structures (pseudo-disks; Galli et al. 2006) are not rotationally supported but may be confused with disks in the absence of kinematical information.

5.3 Multiplicity properties of Class 0 protostars

5.3.1 In the CALYPSO sample

The tentative nature of the multiple dust continuum components we detected in our maps is reported in Table 3. Only the sources that are detected at both frequencies (with separations a > of the synthesized beam scale) are when possible claimed to be robust protostellar candidates (secondary protostar in the sixth column of Table 3). The PdBI maps have a synthesized beam at 231 GHz that can separate systems in the maps with separations larger than ~ 60 au in Taurus (three sources) and~90 au in Perseus and Serpens South (nine sources). For L1157 and GF9-2, although the distance is subject to more uncertainties, we are able to probe systems down to separations ~80−100 au. We exclude the Serpens Main sources from the multiplicity analysis because the larger distance precludes us from probing systems closer than 160 au. The sample of 14 protostars we used for the multiplicity analysis allows us to detect companions down to separations 100 au. The S/N of the 231 GHz dust continuum emission PdBI maps allows us to detect companions of fluxes that are comparable to the flux of the primary protostar in the low-luminosity Taurus sources, or 20 times lower than the flux of the primary protostar for the Perseus and Serpens South sources (except for IRAS2A, for which our high-resolution data are taken from the pilot study of Maury et al. 2010, see Appendix C.4). We are sensitive to multiples in an area around the targeted primary protostar of radius 1500 au in Taurus and 2800 au in Perseus and Serpens South. Of the 14 Class 0 protostars at distances <260 pc in our sample, 8 are single in our dust continuum maps at envelope scales (IRAM04191, L1521F, L1527, L1157, GF9-2, SerpS-MM22, SVS13B, and IRAS2A). However, in the case of IRAS2A, our sensitivity at the longest baselines is much lower than for the rest of the sample (see Appendix C.4), and we do not detect the secondary continuum source found with VLA (Tobin et al. 2015a) and ALMA observations (Maury et al., in prep.) at a separation 140 au. This means that 7 of the 14 closest CALYPSO protostars are single protostars at scales 100−2000 au, within the nominal flux ratio limits given above. The secondary continuum sources detected in the remaining sources (L1448-2A, L1448-N, L1448-C, IRAS4A, IRAS4B, and SerpS-MM18) may be candidate protostellar companions, although in some cases, their nature has not been reliably assessed so far (see Appendix C for details on individual sources). Assuming that all continuum sources detected around the primary sources are genuine protostellar companions, we find that about half of the CALYPSO protostars within 260 pc are located in multiple systems, which leads to a multiplicity fraction ${\text{MF}}_{\text{Class0}}^{1\text{mm}}=57\%$ (with the multiplicity fraction defined as $\text{MF}=\frac{B+T+Q+\cdot }{S+B+T+Q+\cdot }$). Because of the limited sample, excluding even one source would affect the overall multiplicity fraction by ≳ 5%, therefore we did not build distributions of multiplicity fractions with separations.

The multiplicity properties of SVS13B and IRAS4B can be discussed since there are no secondary sources detected within their individual envelopes, but they are associated with widely separated sources that are usually considered as individual protostars in a clustered environment (“separate envelope system” following the classification proposed by Looney et al. 2000). It is unclear whether these sources, which are embedded in the same parsec-scale filamentary structure, will end up in a bound system similar to “ultra-wide pairs” (Joncour et al. 2017) since the future evolution of their proper motions cannot be extrapolated with the current data. We therefore consider that SVS13B and IRAS4B are single at envelope scales. When only the common-envelope systems are considered (excluding SVS13B and IRAS4B from the list of multiples), the Class 0 multiplicity fraction in the CALYPSO sample drops to ${\text{MF}}_{\text{Class0}}^{1\text{mm}}=43\%$.

Regarding the “close” multiple systems, we count three multiple systems at separations a < 210 au: L1448-2Ab, whose nature remains to be confirmed, the secondary source in IRAS2A that we did not detect but that is seen with both VLA (Tobin et al. 2015a) and ALMA (Maury et al., in prep.), and L1448-NB2, which may itself be a binary according to Tobin et al. (2016a). Hence, in our sample we find that only 21% of the protostars within 250 pc have candidate companions at scales 100 < a < 210 au (detection limited to companions of similar flux for the Taurus sources). This paucity of multiple systems at ~ 100 au scales has also been found by Tobin et al. (2016b). We suggest that it may be explained by the physicalconditions in Class 0 envelopes at those scales, which limits the formation of massive ≳ 100 au disks and therefore the fragmentation of gravitationally unstable disks to form binaries at these ~ 100 au scales. The magnetized scenario, proposed in the previous section to explain the paucity of large disks, may explain the paucity of fragments at similar scales as well.

5.3.2 Comparison with other works

Large interferometric surveys of the molecular line emission in complete populations should be carried out to lift the current uncertainties on the multiplicity fraction of protostars, but unfortunately, no complete analysis that would allow unambiguously characterizing the nature of the multiple sources when they are detected is available in a large sample of protostars. When all continuum sources within 0.04 pc in Perseus are counted, as is done in the Tobin et al. (2016b) VANDAM study, many separate-envelope systems are included: whether these will indeed end up as bound systems is debatable, especially because the multiplicity of Class I protostars at similar separations is found to be much lower in their sample (MF = 0.23). We therefore focus on the candidate Class 0 multiple systems in the VANDAM sample that are separated by 50 < a < 5000 au: the multiplicity fraction in VANDAM sources in this separation range is ${\text{MF}}_{\text{Class0}}^{8\text{mm}}=45\%$. This fraction is slightly lower than but consistent with the ${\text{MF}}_{\text{Class0}}^{1\text{mm}}=43-57\%$ that we obtain in the CALYPSO sample at separations 100–5000 au.

5.3.3 Comparison to the multiplicity properties at later stages

A large fraction of stars on the main sequence are observed within multiple systems: the frequency of multiple systems is ${\text{MF}}_{0.7-1.3M_{\odot }}^{\text{MS}}=44\%$ for solar-type stars (see Duchêne & Kraus 2013, for a review). T-Tauri stars and Class II objects have a higher multiplicity fraction than their main-sequence descendants (Kraus et al. 2011). For the direct progenitors of T Tauri stars, Class I YSOs, the multiplicity properties are less clear since these systems are still partly embedded and therefore should be observed in the near- to mid-infrared where spatial resolution is usually lower, and more confusion arises from various types of circumstellar emission. Studies of the multiplicity of Class I protostars in the infrared suggest a lower companion frequency (from ~ 47% ± 8% at separations 14–1400 au to ~26% at separations 110–1400 au and ~16% at separations 300–1400 au; Duchêne et al. 2004, 2007) than in Class II YSOs. This might be due to the narrower range of separations where these studies are complete (between 50 and 200 au). The infrared study by Connelley et al. (2009) found an MF${}_{\text{ClassI}}^{\text{IR}}=44\%$ at separations 50–25 000 au, and a distribution of Class I system separations consistent with a flat distribution. The fact that multiple systems are observed in Class I YSOs suggests that fragmentation processes have taken place, which means that they would have occurred during the earliest stages of star formation (prestellar cores and Class 0 protostars).

We compared our results in the CALYPSO Class 0 protostars with infrared studies of the multiplicity of Class I protostars. Our Class 0 sample is well matched to the scales probed in the “restricted sample” of Connelley et al. (2009), which consists of 32 Class I YSOs at d < 500 pc. The authors discarded binary companions with a projected separation smaller than 50 au for completeness arguments. Sixteen of the 32 Class I protostars they observed have a companion, either wide or close, at separations 50–5000 au (i.e., 50% of the sample). Our multiplicity statistics in the CALYPSO sample (7 single and 8 multiples) at the Class 0 stage, although on an admittedly small sample, is not significantly different from the restricted sample of Connelley et al. (2009). It has been suggested (Reipurth 2000; Sadavoy & Stahler 2017) that a large portion of Class 0 protostars form in non-hierarchical multiple systems that dynamically decay during the Class 0 phase, in which case the overall multiplicity fraction of Class 0 protostars should be higher than that of Class I protostars. Although our conclusions should be viewed as preliminary at this stage because of the small sample and inhomogeneous distances of our sources, our results do not seem to be consistent with this scenario, at least at the scales common in the CALYPSO sample and the Connelley sample since the overall multiplicity is found to be very similar in both samples.

6 Conclusions

In the framework of the CALYPSO survey, we have observed the dust continuum emission in a sample of 16 Class 0 protostars with synthesized beams ~0.4′′ using the Plateau de Bure interferometer at IRAM (50–150 au at the distances of our sources). We performed an analysis of the dust continuum visibility data using models of Plummer envelopes and explored the range of parameters for envelope properties. Six of these 16 protostars cannot be satisfactorily described by a single circularly symmetric envelope model over the whole range of spatial scales (40–1000 au) sampled by our observations at either 1.3 mm or 3 mm. We thus modeled the visibility profiles including anadditional Gaussian component for all sources in our sample in order to test whether adding a disk-like component would explain the properties of some of the CALYPSO protostars at long baselines. The main results of our study are listed below.

• 1.

  For 11 of the 16 protostellar profiles we analyzed, the inclusion of disk-like components improves the description of the 231 GHz continuum visibility profiles at the longest baselines we probed.

• 2.

  However, only four protostars in our sample require a candidate disk structure resolved with radii >60 au to reproduce their dust continuum visibility profiles (IRAS4B, L1527, Serpens South MM22, and Serpens Main SMM4). Two other protostars (L1448-C and GF9-2) are better reproduced when an additional disk-like component is included; they are marginally resolved by our data at radii ≲60 au. In one close multiple system (L1448-NB1/2), an additional circumbinary structure with radius ~200 au may be present, while the two individual protostars do not show evidence of individual disks resolved by our observations.

• 3.

  Four out of the 16 protostars in the CALYPSO sample are likely to harbor large, well-resolved individual protostellar disks at radii >60 au. Our observations suggest that while most Class 0 protostars show evidence of embedded disks, most (>75%) Class 0 continuum disks are small. When we combined all the recent 0.8–1 mm interferometric observations of Class 0 protostars that probed the bulk of the envelope emission down to radii ~50 au, only 5–7 out of 26 Class 0 protostars, that is <28%, show either confirmed or candidate large disks at radii ≳60 au. This upper limit on the occurrence of large disks at the Class 0 stage confirms our earlier results (Maury et al. 2010) and suggests that Class 0 disks are small on average.

• 4.

  From a theoretical point of view, if Class 0 protostars contain similar rotational energy as currently estimated in prestellar cores, only magnetized models of protostellar collapse can reproduce such a large population of small disks during the main accretion phase.

• 5.

  The multiplicity fraction in the CALYPSO sample is found to be ~43−57% ± 5% at scales 100–5000 au, which is slightly larger than but in general agreement with previous studies of Class 0 multiplicity.

Our results suggest that the formation of disks and multiple systems during the Class 0 phase could occur at smaller scales than predicted by hydrodynamical models of rotating protostellar collapse. However, we stress that confirming the properties of the embedded protostellar structures requires additional spectral line analysis that either traces rotationally supported motions to robustly identify disk components, or confirms the protostellar nature of the systems at ~50 au scales.

Acknowledgements

We thank the IRAM staff for their support in carrying out the CALYPSO observations, and the anonymous referee for a detailed reading that significantly helped us to improve the paper. We thank Kees Dullemond, Ralf Klessen, Sandrine Bottinelli, Benoit Commercon, and Ralf Launhardt for their participation to the preparation of the CALYPSO project. We thank the INSU “Action Spécifique ALMA” for their financial support to the CALYPSO collaboration. This work has benefited from the support of the European Research Council under the European Union’s Seventh Framework Programme (Advanced Grant ORISTARS with grant agreement No. 291294 and Starting Grant MagneticYSOs with grant agreement No. 679937).

Appendix A Details of CALYPSO observations

Appendix B CALYPSO dust continuum maps

Fig. B.1 1.3 mm (231 GHz) and 3.3 mm (94 GHz) PdBI dust continuum emission maps of L1448-N. The ellipses in the bottom left corner show the respective synthesized beam sizes. The contours are levels of −3σ (dashed) and 5σ, then in stepsof 10σ from 10σ to 100σ, and finally in steps of 20σ beyond (rms noise computed in the map before primary beam correction as reported in Table 2). The maps have been corrected for primary beam attenuation. The blue and red arrows show the direction of the protostellar jet(s) associated with the millimeter sources in the field, see Table 1 for further details.

Fig. B.2 Same as Fig. B.1 for L1448-C. The contours show levels of −3σ (dashed), 5σ, and 10σ and then usesteps of 40σ from 20σ on.

Fig. B.3 Same as Fig. B.1 for IRAS2A. The contours show levels of −3σ (dashed), 5σ, and 10σ and then usesteps of 20σ from 20σ on.

Fig. B.4 Same as Fig. B.1 for IRAS4A. The contours show levels of −3σ (dashed), 5σ, and 10σ and then usesteps of 20σ from 20σ on. IRAS4B is detected to the southwest in the 94 GHz map, outside the area shown here, so that the structure of IRAS4A is well distinguished in the figure.

Fig. B.5 Same as Fig. B.1 for IRAS4B. The contours show levels of −3σ (dashed), 5σ, and 10σ and then usesteps of 30σ from 20σ on. The fluxdensity in the maps has been corrected for primary beam attenuation (the primary beam FWHM is shown as a pink dashed circle in the 94 GHz map).

Fig. B.6 Same as Fig. B.1 for IRAM04191. The contours show levels of −3σ and 5σ and then usesteps of 10σ from 10σ on.

Fig. B.7 Same as Fig. 1 for L1521F. The contours show levels of −3σ and 5σ and then usesteps of 10σ from 10σ on.

Fig. B.8 Same as Fig. B.1 for SerpM-S68N. The contours show levels of −3σ and 5σ and then usesteps of 10σ from 10σ on. The dashed pink contour in the 94 GHz continuum emission map shows the PdBI primary beam.

Fig. B.9 Same as Fig. B.1 for SerpM-SMM4. The contours show levels of −3σ, 5σ, and 10σ and then usesteps of 20σ from 20σ on.

Fig. B.10 Same as Fig. B.1 for SerpS-MM18. The contours show levels of −3σ and 5σ and then usesteps of 20σ from 10σ on. The dashed pink contour shows the PdBI primary beam.

Fig. B.11 Same as Fig. B.1 for SerpS-MM22. The contours show levels of −3σ and 5σ and then usesteps of 10σ from 10σ on.

Fig. B.12 Same as Fig. B.1 for L1157. The contours show levels of −3σ and 5σ and then usesteps of 30σ from 20σ on.

Fig. B.13 Same as Fig. B.1 for GF9-2. The contours show levels of −3σ, 5σ, and 10σ and then usesteps of 20σ from 20σ on.

Appendix C Comments on individual sources

C.1 L1448-2A

Single-dish constraints

We used the 268 GHz Bolocam fluxes reported in Enoch et al. (2006), and subtracted the 15 mJy from L1448-2Ab, which was removed from the PdBI visibilities. We used the dual-frequency PdBI spectral index at 20 kλ (see Table 4) to extrapolate the flux from 268 to 231 GHz. The uncertainties on the 231 GHz extrapolated fluxes are ± 30%. The same method was followed to extrapolate the 94 GHz envelope fluxes (integrated flux and peak flux), with estimated uncertainties of ± 40%.

Multiplicity

L1448-2A is the only millimeter source detected in the field, see Fig. 1, but an extension is also detected west of the millimeter peak position. The coordinates of this secondary source, L1448-2Ab, are reported in Table 3. L1448-2Ab was detected with the VLA (Tobin et al. 2016b) with a peak flux density 0.151 Jy at 8mm, but was never previously detected at higher frequencies, including in the 1.3 mm CARMA map at 0.3′′ resolution published by Tobin et al. (2015b). Its location in the equatorial plane of L1448-2A (as estimated from the jet position angle from L1448-2A), detection in independent datasets at different millimeter wavelengths, and 1.3 mm peak flux half of the flux toward the primary protostar make it a robust protostellar companion candidate. Moreover, we tentatively detect two distinct jets from L1448-2A in our CALYPSO maps of the high-velocity ¹²CO(2-1) emission (blue and red jets are misaligned; Podio et al., in prep.). The exact nature of L1448-2Ab needs to be further characterized with observations at higher angular resolution, but we consider it a candidate protostellar companion in our multiplicity analysis.

Candidate disk

The visibility profiles at 231 and 94 GHz are shown in Figs. 5 and C.1 together with the best Plummer-only (Pl) and Plummer + Gauss (PG) models. When all the dust continuum visibilities are used, the best-fit model to reproduce the continuum emission visibility profiles of L1448-2A is the Plummer-only model for both frequencies (see main text and Table 6). Thebest-fit PG model for the 231 GHz visibility profile includes a 12 mJy Gaussian source with an unresolved FWHM <0.15′′. The upper limits on the parameters for this Gaussian component are reported as upper limits on the disk properties in Table 5. Similarly, the best-fit PG model that reproduces the 94 GHz visibilities includes the minimum (size and flux) Gaussian component allowed in the fitting procedure (see Table 6). When only the equatorial visibilities are used (i.e., only the visibilities at PA − 80° ± 30° are selected), we find that the best-fit model is still the Plummer-only model (see models Pleq and PGeq in Table C.1). It is therefore clear from our analysis of the CALYPSO data that the continuum structure around L1448-2A traces the inner part of the envelope, and no resolved continuum disk-like emission is detected at scales 50–500 au.

Fig. C.1 Visibility profiles of the 231 GHz (black dots) and 94 GHz (gray triangles) continuum emission from L1448-2A. In each visibility curve, we show the best-fit Plummer-only model (Pl, black curve) and the best-fit Plummer + Gauss (PG, red curve) model. In the case of L1448-2A, the F-test suggests that the PG models are not statistically better than the Plummer-only (Pl) models, although the presence of an unresolved (radius < 0.15′′) disk-like component cannot be robustly excluded.

C.2 L1448-NB

Single dish constraints

We used the data from IRAM-30 m MAMBO observations reported in Motte & André (2001), and subtracted the fluxes from the sources that were removed from the PdBI visibilities, that is, L1448-NA and the secondary western source NB2. To extrapolate the MAMBO 243 GHz fluxes to our observing frequencies, we used the large-scale spectral index computed directly from our PdBI dual-frequency visibility amplitudes at 20 kλ (see Table 4). While the integrated flux over 4200 au is 2.1 Jy in Motte & André (2001; MA01), we would expect a total envelope flux 1.79 Jy in a 14′′ radius (the distance adopted in MA01 was 300 pc). However, the modeled R_out of 10 000 au in Motte & André (2001) for L1448-NB suggests a larger radius (20′′ at 300 pc) than the radius in which they report this integrated flux. From their map, we computed the flux in an area of radius 20′′, obtaining an integrated flux 3.7 Jy, which translates into a 3.24 Jy total envelope flux at 231 GHz. Similarly, we used the MA01 peak flux to complete our visibility profile at 9 kλ (for a 20′′ radius source in a 11′′ beam). The estimated uncertainties on the 231 GHz extrapolated fluxes are ± 30%. The same method was followed to extrapolate the 94 GHz envelope fluxes (integrated and peak fluxes), and the estimated uncertainties are ± 40%.

Multiplicity

L1448-NB is the strongest millimeter source in the field, see Fig. B.1. The secondary source L1448-NA, located 1500 au away, was previously detected at both infrared and millimeter wavelengths and is classified as a Class I protostar (Ciardi et al. 2003; Kwon et al. 2006). It drives an outflow whose redshifted lobe was detected by Lee et al. (2015) at a PA 218°. The source L1448-NW was also previously detected in the BIMA and SMA surveys (Looney et al. 2000; Lee et al. 2015), although it is sometimes referred to as L1448-IRS3C: it is associated with two 8 mm sources, might not be gravitationally bound to the NA/NB system, and seems to drive an outflow (Lee et al. 2015), but its nature is not precisely determined.

At smaller scales, NB is resolved into two components, NB1 and NB2. NB1 is the strongest millimeter source, and the western secondary component L1448-NB2 was previously detected in the 1.3 mm dust continuum emission with the PdBI and ALMA maps (Maury et al. 2015; Tobin et al. 2016a) and with the VLA at 8 mm (Lee et al. 2015). It is tentatively resolved into two components called IRS3B-a and IRS3B-b with ALMA and VLA. Its location along the axis of the collimated jet driven by L1448-NB1 (projected onto the plane of the sky) and its flat spectral index (possibly due to free-free emission) would naively suggest that this source is due to interaction between the jet and the surrounding material. However, the presence of a small-scale structure (at scales 200–300 au, see our modeling of the visibility profile for this source described below) surrounding NB1 and NB2 can be used to argue that NB2 does not trace a jet feature, but possibly a protostellar companion in a circumbinary disk (this point is further discussed below). A robust assessment of the nature of this millimeter continuum component will require analyzing a larger set of data, including molecular lines and tracers of protostellar nature, which is beyond the scope of this paper. To build robust estimates of the upper-limit MF from our sample, we consider it a candidate protostellar companion in our multiplicity analysis.

Candidate disk

• Centered on L1448-NB1: Our analysis is based on the visibilities that only contain the NB1 primary protostar, after removing the secondary sources (NA and NB2 in the 231GHz data, and also NW in the 94 GHz map: as an illustrating example, Fig. C.2 shows the maps we obtained after removing these components, to be compared with the maps shown in Fig. B.1). We show that adding a Gaussian structure to the Plummer envelope model centered on NB1 does not improve the modeling of the 231 GHz visibility profile (see Table 6 and Fig. C.3). Regarding the 94 GHz visibility profile centered on NB1, our analysis suggests that an additional Gaussian component improves the modeling: the F-value computed between the best-fit Pl and PG models is 15, while the critical value (corresponding to a probability 0.3% that the Plummer + Gaussian model is better only because it includes two more free parameters) is 8 in this case. None of the models for the 94 GHz visibility profile is considered satisfactory, however, since none of their reduced χ² values drop below 3. The discrepancy between the 231 GHz and 94 GHz models might either arise because the 94 GHz observations do not separate the two millimeter continuum components L1448-NB1/NB2 and detect the additional component associated with L1448-NB2 as an additional Gaussian, while we could remove L1448-NB2 from the 231 GHz observations that do separate the two components. It might also suggest that an additional Gaussian component is indeed present in the system, but is not centered on the main continuum source L1448-NB1. If it is slightly shifted toward the western source L1448-NB2, this shift would not greatly affect the 94 GHz visibility profile, which is insensitive to asymmetries at scales <1′′, while our analysis of the 231 GHz visibility profile centered on the main millimeter source L1448-NB1 would see this circumbinary structure as an asymmetry. Strengthening this hypothesis, the 231 GHz PdBI visibility profile shows oscillations in the circular bins in the 231 GHz visibility profile at baselines >100 kλ. This suggests that there are scales at which the spatial distribution of the emission is not circularly symmetric around the chosen phase center, here L1448-NB1. However, even if we use the equatorial plane visibilities, we still find that the best-fit models are the Plummer-only models.

• Centered on L1448-NB2: Recent ALMA observations by Tobin et al. (2016a) have shown a continuum structure that resembles a disk structure. The authors suggested that it might be centered on L1448-NB2 (also called IRS3B-a in their paper): they argued that although NB2 is not the strongest millimeter source in the system, its mass dominates, and hence the Keplerian motions are centered around the NB2 secondary source. Following these results, we wished to test the hypothesis of an additional disk component centered on L1448-NB2 that would be missed by our analysis, which is centered on the main millimeter source (NB1). We shifted the phases of the PdBI visibilities so that they were centered on the L1448-NB2 component and modeled the continuum visibilities toward L1448-NB2. To do this, we had to remove the brightest millimeter source L1448-NB1 from the visibility tables because otherwise its flux dominated (we subtracted a model point-like source with a flux of the peak flux of L1448-NB1 at the position of L1448-NB1). Centering the Plummer envelope model on NB2 improves the minimization, but for both the 94 GHz and the 231 GHz data, the best-fit models are the Plummer-only models (see models Pl and PG in Table C.3). The dust continuum emission in L1448-NB is not circularly symmetric: if we only use the equatorial plane visibilities, however, we find that the best-fit models are still the Plummer-only models.

• Centered at the barycenter: As an ultimate test of our modeling, we tested the possibility that a circumbinary disk-like additional structure might surround the two main millimeter sources L1448-NB1 and NB2. We removed the point-like contributions from the three protostars L1448-NA, L1448-NB1, and L1448-NB2, which produced profiles of the visibility amplitudes of the remaining underlying structure that we phase-shifted to be centered in between the two millimeter sources. Overall, the continuum emission around the L1448-NB1/NB2 system is only satisfactorily modeled at 94 GHz by either a Plummer envelope or an envelope with an additional Gaussian component (see Table C.4). The PG model is better in this case, however, and points toward the presence of a ~1.1′′ (260 au radius) additional structure with a flux of 50 mJy at 94 GHz, which centered in between the NB1 and NB2 protostars. The strong asymmetries of the continuum emission at 231 GHz when the two strongest sources in the map are excluded preclude us from performing a more advanced modeling of the PdBI continuum visibilities to further explore the nature of the circumbinary structure that emits in the millimeter continuum.

Fig. C.2 1.3 mm (231 GHz) and 3.3 mm (94 GHz) PdBI dust continuum emission maps of L1448-NB1 after subtracting the secondary components. The contours are the same as in Fig. B.1.

We conclude that while the individual protostars NB1 and NB2 do not harbor individual disk-like structures that can be resolved with our PdBI observations, they might both be embedded within a candidate circumbinary disk-like structure at scales ~200 au. The nature of this additional structure needs to be investigated further, especially since its disk nature, suggested in Tobin et al. (2016a), is questioned because it does not seem to be peaking on any of the two bright millimeter sources (that are assumed to contain most of the mass), and no indication that the material in this structure is rotationally supported has been found so far (a clear velocity gradient is detected in CALYPSO observations by Maret et al., in prep., and Gaudel et al., in prep., but it is not well reproduced by Keplerian rotation). It is possible that the structure traces tidal arms that are created by a differential gravitational potential due to relative motions of the multiple components.

C.3 L1448-C

Single-dish constraints

From their IRAM-30 m MAMBO observations, Motte & André (2001) found the following fluxes for the L1448-C envelope: F_peak = 620 mJy and F_int = 910 mJy in a 4200 au radius (assuming d = 300 pc, which translates into a 14′′ radius on the map). We scaled down the IRAM-30 m fluxes, which are obtained at an observing frequency 243 GHz, to 94 and 231 GHz using the PdBI spectral index at short baselines (see Table 4). We removedone additional component from the PdBI visibilities: L1448-CS, but it is out of the single-dish beam, therefore wedid not remove its flux from the extrapolated flux. For the 231 GHz profile, we used an integrated flux $(910\times {(231/240)}^{2.5})$ = 827 mJy, and a peak flux F_peak = 563 mJy, at 8 kλ. We used a 20% uncertainty on both extrapolated fluxes (this sets the upper- and lower-limit values that the total flux parameter is allowed to take in the fitting procedure). At 94 GHz, the peak flux extrapolated from Motte & André (2001) is consistent with the Looney et al. (2003) BIMA flux at short baselines (65 mJy at 2 kλ). We used a 40% uncertainty on both the integrated flux and the peak flux we used to model the 94 GHz profile.

Fig. C.3 Visibility profiles of the 231 GHz (black dots) and 94 GHz (gray triangles) continuum emission from L1448-NB1.Open symbols are used when the visibility real part has a negative value (since absolute values are shown in the log–log plot). In each visibility curve, we show both the best Plummer-only model (Pl, black curve) and the best Plummer + Gauss (PG, red curve) model. In the case of L1448-NB1, the F-test suggests that the PG model for the 231 GHz profile does not perform statistically better than the Plummer-only (Pl) model (see Table 6 for more information on the two models). For the 94 GHz profile, the PG model that includes a 0.94′′ Gaussian component performs better than the Pl model, although none of the two models does satisfactorily reproduce the visibility profile (reduced chi square > 3 in both cases).

Multiplicity

Our CALYPSO dust continuum emission maps are shown in Fig. B.2. The primary protostar is well detected at the center of the field, and we detect continuum emission associated with L1448-CS, 8′′ (2000 au) southeast, at both frequencies. This source has previously been detected at millimeter wavelengths (Jørgensen et al. 2007; Maury et al. 2010; Hirano et al. 2010), and is associated with a mid-infrared source seen with Spitzer (Jørgensen et al. 2006). This southern source is brighter than L1448-C in the mid-infrared, but much weaker in the millimeter and submillimeter bands: hence it is probably a more evolved source (Class I or older), with a separate envelope.

Fig. C.4 Visibility profiles of the 231 GHz (black dots) and 94 GHz (gray triangles) continuum emission from L1448-C. In each visibility curve, we show both the best Plummer-only model (Pl, black curve) and the best Plummer + Gauss (PG, red curve) model. The PG model that includes a marginally resolved 0.16′′ -FWHM Gaussian component reproduces the 231 GHz profile better than the Plummer-only (Pl) model (see Table 6 for more information on the two models). Moreover, the PG model that includes an unresolved additional Gaussian component reproduces the 94 GHz visibility profile better.

Candidate disk

The 231 GHz visibility profiles and models for L1448-C are shown in Fig. C.4. When all the dust continuum visibilities are used, the best-fit model to reproduce the continuum emission visibility profiles of L1448-C is the Plummer + Gaussian model for both frequencies (see Table 6). The best-fit PG model for the 231 GHz visibility profile includes a 130 mJy Gaussian source with an FWHM 0.16′′, which is marginally resolved by our observations. The parameters for this Gaussian component are reported as the candidate disk properties in Table 5. Similarly, the best-fit PG model to reproduce the 94 GHz visibilities includes an unresolved (FWHM < 0.3′′) Gaussian component with a flux of 18 mJy (see Table 6). When only the equatorial visibilities are used (in a direction orthogonal to the jet axis position angle, i.e., selecting only the visibilities at PA − 107° ± 30° for L1448-C), we find that the best-fit model to reproduce the 231 GHz profile is still the Plummer + Gaussian model with a size and flux similar to the best- fit model value using all visibilities (see models Pleq and PGeq in Table C.5). Hence, our analysis of the CALYPSO data suggests that a candidate disk is detected in L1448-C, and that it is marginally resolved at radii ~ 40−50 au.

C.4 NGC1333 IRAS2A

Single-dish constraints

In theirMAMBO observations, Motte & André (2001) found a peak flux 875 mJy beam⁻¹ in the 11′′ beam of the IRAM-30 m, and they claimed that the source is almost unresolved. Hence, we used this peak flux as the envelope-integrated flux for IRAS2A. Using the 20 kλ spectral index from our dual-frequency continuum visibilities (see Table 4), we scaled this value down to our frequencies to obtain the total envelope fluxes at 231 and 94 GHz. We left the source size quite loose in our fitting and allowed values for R_out ~ 6 ± 4′′ since it is unresolved by the single-dish observations. The uncertainty on the total flux was set to 30%.

Multiplicity

In our CALYPSO maps, shown in Fig. B.3, a single protostar is detected. The two continuum sources IRAS2A2 and IRAS2A3, reported in Maury et al. (2014) and Codella et al. (2014b), are now shown to originate from an envelope structure that is detected by the interferometer: our improved dataset after self-calibration (which improved the rms noise and the imaging fidelity) allows us to recover the lower surface brightness emission that surrounds them, and they are found to reconnect with the main source envelope emission. This suggests that, similar to what was proposed in Santangelo et al. (2015) for IRAS4A, these sources represent the dust continuum emission from envelope structures that are detected by the interferometer when only a restricted range of spatial scales is sampled. They are therefore no longer considered as robust compact continuum emission components. Recent VLA (Tobin et al. 2015b) and ALMA observations (not yet published, Maury et al., in prep.) of IRAS2A have revealed the presence of a secondary millimeter dust continuum source located 0.4′′ from the main millimeter source. Although the nature of this secondary component remains to be investigated, it is likely that IRAS2A is a close (≤ 100−200 au) binary system (two separate jets are detected; Podio & CALYPSO, in prep.). This source is not separated from the main millimeter source in our PdBI maps, although a slight extension is detected at its location. Our highest angular resolution data for IRAS2A (configuration A of the PdBi at 231 GHz) was obtained as part of the pilot R068 project (Maury et al. 2010): the lower sensitivity of these earlier observations (obtained before the installation of the WideX correlator) explains the non-detection of a separated compact source at the secondary position. To build robust estimates of the upper-limit MF from our sample, and considering that a secondary source is detected with ALMA and the VLA, we count IRAS2A as a binary system in our multiplicity analysis.

Fig. C.5 Same as Fig. C.4 for IRAS2A. The Plummer-only model is the best-fit model to reproduce the 231 GHz profile. The PG model, which includes an unresolved additional Gaussian component, performs marginally better than the Pl model for the 94 GHz visibility profile.

Candidate disk

The 231 GHz visibility profiles and models are shown in Fig. C.5 for IRAS2A. When all the dust continuum visibilities are used, the best-fit model to reproduce the continuum emission visibility profiles of IRAS2A is the Plummer-only model for the 231 GHz profile (see main text and Table 6). The best-fit PG model for the 231 GHz visibility profile that satisfactorily reproduces the visibility profile but is not better than the Plummer-only model includes a 52 mJy Gaussian source that is unresolved with an FWHM 0.01′′ (the smallest size allowed in our fitting procedure). The upper limits on the parameters for this Gaussian component are reported as upper limits on the disk properties in Table 5. The PG model performs marginally better than the Pl model at reproducing the 94 GHz visibilities (a comparison of the two models produces an F-value of 11, while the critical value of F is 8. Above this, the probability that the PG model is better only by chance is <0.3%). It includes an unresolved Gaussian with a flux of 9 mJy (see Table 6). When only the equatorial visibilities are used (in a direction orthogonal to the jet axis position angle, i.e., selecting only the visibilities at PA − 65° ± 30°), we find that the best-fit model is still the Plummer-only model for the 231 GHz profile (see models Pleq and PGeq in Table C.6). Our analysis of the CALYPSO data therefore suggests that no resolved disk-like emission is detected at scales 50–500 au.

We also tested the effect of removing a point source from the flux of the secondary source that is resolved with the ALMA 1.3 mm observations (14 mJy at 230 GHz; Maury et al., in prep.) from the PdBI visibilities and shifting the phase center to the position of the main millimeter source at (03:28:55.569;31:14:36.952) to model the 231 GHz PdBI visibilities anew. For the sake of clarity and brevity, we do not report this model here, but the best-fit model is still the Plummer-only model, with very similar parameters as were found as best fits for the whole visibility dataset, which is reported in Table 6.

C.5 SVS13B

Single-dish constraints

Fluxes obtained using the MAMBO bolometer array on the IRAM-30 m telescope are reported in Chini et al. (1997): F_peak = 900 mJy and F_int = 1180 mJy. However, Lefloch et al. (1998) reported a peak flux for SVS13B of only 320 mJy beam⁻¹, with the same instrument and telescope. We used the Chini et al. (1997) peak flux as the envelope integrated flux in a source the size of the MAMBO beam (11′′-FWHM), and we used the Lefloch et al. (1998) value as the lower limit allowed for the total envelope flux for the minimization. We used the peak flux as an integrated flux in the MAMBO beam because contamination from the surrounding filament, and more especially, from the nearby SVS13A, precluded obtaining robust integrated fluxes for SVS13B alone at radii larger than 6′′. We scale these fluxes down using the 20 kλ spectral index from our dual-frequency continuum visibilities (see Table 4). We also let the source size quite loose in the minimization procedure, allowing radii from 4′′ to 14′′. The uncertainty on the total envelope flux was set to 40%.

Multiplicity

SVS13B is in a large-scale multiple system (separate envelopes) with SVS13A, a Class I protostar located 3500 au away, and SVS13C, whose nature is as yet undetermined, and which lies 4500 au away. These three sources are detected in our 94 GHz dust continuum emission map shown in Fig. 2. A fourth source, VLA3 detected previously at centimeter wavelengths with the VLA by Rodríguez et al. (1999), is also detected in both our maps (although not at the 10σ level in our 231 GHz data). A millimeter counterpart to VLA3 was also previously tentatively detected in the BIMA maps of Looney et al. (2000) as an extension from the SVS13A source at 110 GHz (they named it A2). The primary Class 0 protostar in the field, SVS13B, does not exhibit any sign of further multiplicity down to the scales of ~50 au that are probed by our data. A single protostellar jet from SVS13B is detected in our CALYPSO maps of molecular line emission at a PA + 167° (Podio & CALYPSO, in prep.).

Fig. C.6 Same as Fig. C.4 for SVS13B. The Plummer-only models are the best-fit models to reproduce the two profiles, although at 231 GHz, a more reasonable (p + q) index is found in the PG model that includes a 0.19′′ FWHM Gaussian component (flux 80 mJy).

Candidate disk

When all the dust continuum visibilities are used, the best-fit model to reproduce the continuum emission visibility profiles of SVS13B is the Plummer-only model for both frequencies (see main text and Table 6). However, the best Plummer-only model to reproduce the 231 GHz visibility profile is found to have p + q = 2.9 which seemsexcessively high. Although we do not exclude that such a value might be physically possible, we stress that including a Gaussian component allows reducing this slope index to p + q = 2.5, which is a more standard value that is expected in protostellar envelopes (Looney et al. 2003; also found p + q = 2.4 at short < 40 kλ baselines in their 3 mm BIMA data). This best-fit PG model for the 231 GHz visibility profile includes a 80 mJy Gaussian source with an FWHM 0.19′′: although this model is not statistically better than the Plummer-only model, we have to keep these facts in mind about SVS13B. In Table 5 we report these values as upper limits for the candidate disk component. The F-test shows that the best-fit model to reproduce the 94 GHz visibilities is also the Plummer-only model (see Table 6 and Figs. 6 and C.6), with a satisfactorily reduced χ² value that isonly slightly lower than that of the Plummer + Gaussian model, but has two free parameters less. When only the equatorial visibilities are used (in a direction orthogonal to the jet axis position angle, i.e., selecting only the visibilities at PA − 103° ± 30° from the jet axis found in the CALYPSO data; Podio et al., in prep.), we find that the best-fit model is still the Plummer-only model at 231 GHz with a high index p + q = 2.9 (see models Pleq and PGeq in Table C.7). Our analysis of the CALYPSO data therefore suggests that the continuum emission from SVS13B is better reproduced by an envelope model down to scales 50 au, but we cannot exclude the presence of a ≲ 60 au candidate disk which in addition produces more reasonable p + q values.

C.6 NGC1333 IRAS4A1

Single-dish constraints

Motte & André (2001) reported that the NGC1333 IRAS4A envelope is unresolved in their IRAM-30 m observations (12′′-FWHM beam). The peak flux is F_peak = 4.1 Jy in the IRAM-30 m beam, which we used as the total envelope flux in a 6′′ source. We removed the contribution from IRAS4A2 and rescaled the MAMBO flux that was obtained at a central frequency of 243 GHz to obtain extrapolated fluxes at both 94 and 231 GHz (using the spectral index at large scales from our dual-frequency PdBI data, see Table 4). We used a 20% uncertainty on the single-dish flux, and let the outer envelope radius vary between 3′′ and 9′′.

Multiplicity

The IRAS4A system is resolved by our CALYPSO observations: see Fig. B.4. The secondary protostar IRAS4A2 was extensively discussed by Santangelo et al. (2015) and drives its own high-velocity jet. The continuum emission associated with the component IRAS4A3 reported in Santangelo et al. (2015) disappears when robust weighting is used. This confirms the suggestion made in Santangelo et al. (2015) that this traces a structure of dust continuum emission that is produced by the outflow interaction with the envelope, and is not a compact component associated with a true protostellar source.

Fig. C.7 Same as Fig. C.4 for IRAS4A1. The Plummer-only models are the best-fit models to reproduce the two visibility profiles.

Candidate disk

The models and visibility profiles of IRAS4A1 are shown in Fig. C.7. When all the dust continuum visibilities are used, the best-fit model to reproduce the continuum emission visibility profiles of IRAS4A1 is the Plummer-only model for both the 231 and 94 GHz visibility profiles (see main text and Table 6). The best-fit PG model for the 231 GHz visibility profile includes a 348 mJy Gaussian source with an FWHM 0.3′′, but does not perform better than the Plummer-only model (slightly larger reduced χ², but two additional free parameters). The upper limits on the parameters for this Gaussian component are reported as upper limits on the disk properties in Table 5. When only the equatorial visibilities are used (in the direction of the C¹⁷O(3–2) velocity gradient that is shown in the SMA map by Ching et al. (2016) at the 4′′ core scale, i.e., selecting only the visibilities at PA − 62° ± 30°), we find thatthe best-fit model is still the Plummer-only model for the 231 GHz profile (see models Pleq and PGeq in Table C.6). It is therefore clear from our analysis of the CALYPSO data that the continuum structure around NGC 1333 IRAS4A1 can be well described with the inner part of the envelope: no resolved disk-like emission is detected at scales 50–500 au around A1. Finally, we note that the brightness temperature obtained from the 231 GHz flux density of IRAS4A1 is 44K in the ~ 0.5′′ beam: this is about the temperature expected from $T_{\text{dust}}=38\times L_{\text{int}}^{0.2}\times {(r/100\hspace{0.17em}\text{au})}^{-0.4}~60$ K at r ~ 100 au. This suggests that the dust continuum emission might be partially optically thick at scales of about the size of the synthesized beam (30–80) au. Using the optically thin 3.2 mm flux density in the ~ 1′′ beam, we deduce a column density 0.1−4 × 10²⁶ cm⁻² (using κ_3.2mm = 0.017 from the spectral index computed at 20 kλ between the 1.3 and 3.2 mm visibility amplitudes). We acknowledge that an optically thick disk could show a shallow flux decline at long baselines that might be confused with the inner part of an optically thin envelope. In our sample, the 231 GHz dust continuum emission is optically thin at all scales probed in all sources except for IRAS4A: the possible confusion caused by the partial optical thickness of the 231 GHz emission in IRAS4A is mitigated by the fact that (i) both visibility profiles at 231 GHz and 94 GHz show a steep decline at long baselines and (ii) the 94 GHz dust continuum emission is optically thin and is well reproduced by an envelope model down to scales 50 au. We therefore argue that despite the partial optical thickness of the 231 GHz emission in IRAS4A1, our non-detection of a candidate disk structure in IRAS4A1 is robust.

We note that the neighboring IRAS4A2 probably includes a small candidate disk structure. When we removed the A2 component to produce clean visibilities for A1 (with the aim of modeling the A1 structure), we realized that our PdBI data are best modeled when a circular disk component of diameter 0.7′′ and flux 430 mJy is used at both 231 and 221 GHz independently for A2: this suggests that a disk-like component of radius 0.3−0.5′′ is associated with IRAS4A2. This is in agreement with observations by Choi et al. (2007, 2010), who showed maps of the blueshifted and redshifted emission of the NH₃(3,3) line around IRAS4A2, which they interpreted as tracing the rotation of a disk at PA ~ 109°, while no such rotation is detected around IRAS4A1.

C.7 NGC1333 IRAS4B1

Single-dish constraints

The IRAM-30 m observations of the protostellar envelope with the MAMBO camera at 243 GHz (Motte & André 2001) suggest that the IRAS4B envelope is unresolved in the IRAM-30 m 11′′ beam, with a peak flux 1470 mJy beam⁻¹. We scaled this flux density down from 243 to 231 GHz and 94 GHz using the PdBI spectral index at 20 kλ (see Table 4). However, the 94 GHz extrapolated flux is lower than the flux obtained with our PdBI observations at 10 kλ, so that we used the Looney et al. (2000) BIMA flux at 108 GHz (200 mJy at 2–8 kλ) to extrapolate an envelope flux at 94 GHz. We also note that these authors find a best-fit envelope size of ~2000 au using a distance of 350 pc (i.e., 6′′ radius). This confirms that the envelope is mostly unresolved in the IRAM-30 m beam. We recovered the entire single-dish flux at 20 kλ in our observations at 231 GHz, which suggests that the envelope of IRAS4B is very compact and is fully probed by the CALYPSO observations We used a 40% uncertainty on the single-dish fluxes, and let the envelope outer radius vary between 3′′ and 9′′.

Multiplicity

The secondary source IRAS4B2 (also called IRAS 4BE, IRAS 4B′, IRAS4C, and IRAS 4BII; see Looney et al. 2000) is detected 10′′ east (see Fig. B.5). It is outside of the IRAS4B1 envelope. Our PdBI CALYPSO data do not detect COM emission in this secondary source (De Simone et al. 2017), and IRAS4B2 is not detected at 70 μm in the Herschel HGBS maps (Sadavoy et al. 2014). This suggests that it is either a very low luminosity source or a very young protostar (candidate first hydrostatic core). The SMA MASSES survey (Lee et al. 2016) and CARMA TADPOL survey (Hull et al. 2014) have recently suggested that some high-velocity ¹²CO emission could trace outflowing gas associated with this secondary millimeter source. Moreover, our CALYPSO observations detect high-velocity blueshifted SO emission originating from IRAS4B2 (Podio & CALYPSO collaboration, in prep.) at PA ~ −99°. Since B2 is not detected in the infrared, it cannot currently be characterized well enough to firmly establish a robust protostellar nature. To build robust estimates of the upper-limit MF from our sample, we consider that IRAS4B is a separate-envelope system in our multiplicity analysis.

Candidate disk

When all the dust continuum visibilities are used, the best-fit model to reproduce the 231 GHz continuum emission visibility profiles of IRAS4B1 is the Plummer + Gaussian model (see main text and Table 6). It includes a 645 mJy Gaussian source with an FWHM 0.53′′. Although the modeling is not satisfactory because the power-law index (p + q) is unrealistically high (2.9), we report the parameters for this Gaussian component as the size and flux of the candidate disk in IRAS4B in Table 5. The PG and Pl models can both reproduce the 94 GHz visibilities, but they also feature an unrealistically high power-law index (2.9; see Table 6). When only the equatorial visibilities are used (in a direction orthogonal to the jet axis position angle, i.e., selecting only the visibilities at PA − 103° ± 30° for IRAS4B1), we find the Plummer + Gaussian model performs slightly better at reproducing the 231 GHz visibility profile, with a Gaussian size and flux similar to the best-fit model values found using all visibilities (see models Pleq and PGeq in Table C.9). We stress that the envelope of IRAS4B is very compact at both frequencies, and its spatial extent seems fully probed by our CALYPSO observations (see Fig. C.8). The IRAS4B visibility profiles are very flat, especially at 94 GHz: such a steep slope has previously been noted in the BIMA observations by Looney et al. (2003) ((p + q) ~ 2.8 at 108 GHz between 5 and 90 kλ). We also find it striking that half of the envelope dust continuum emission flux is included in a Gaussian-like structure: when a dust temperature of 50 K is used, such a 231 GHz flux translates into a very high “disk” mass of 0.2−0.4 M_⊙ depending on the assumptions made on the dust emissivity. This raises questions on the nature of the dust continuum emission in IRAS4B, and we stress that were it a candidate disk, it should show some evident rotational signature that currently is not observed at ~arcsecond scales (Yen et al. 2015b). We conclude that our analysis of the CALYPSO data suggests that a candidate disk-like structure might be detected at radii ~120 au in IRAS4B1, but our current modeling does not allow us to robustly conclude on the nature of the dust continuum emission recovered by the PdBI in this source.

C.8 IRAM04191

Single-dish constraints

In the IRAM-30 m MAMBO observations by Motte & André (2001) at 240 GHz, the envelope is resolved with an integrated flux 650 mJy in a 4200 au radius. The peak flux is 110 mJy in the 11′′ IRAM-30 m beam. We scaled these fluxes down from 243 to 94 GHz and 231 GHz using the PdBI spectral index at 20 kλ (see Table 4). We used a 20% uncertaintyon the single-dish fluxes, and let the envelope outer radius vary between 20′′ and 40′′.

Fig. C.8 Same as Fig. C.4 for IRAS4B1. The PG model that includes a 0.53′′ Gaussian source (flux 645 mJy) reproduces the 231 GHz visibility profile better, while the 94 GHz visibility profile is satisfactorily reproduced by the Plummer-only model.

Fig. C.9 Same as Fig. C.4 for IRAM04191. Open symbols are used when the visibility real part has a negative value (since absolute values are shown in the log-log plot). The red Plummer + Gaussian (PG) model that includes an unresolved 0.17′′ -FWHM Gaussian component (see Table 6) is only statistically better than the black Plummer-only (Pl) model at 231 GHz. The 94 GHz profile can be well modeled with a Plummer-only model. The oscillations of the models are due to the Hankel transform of the power-law envelope model with a low index (p + q = 1.4).

Multiplicity

The secondary source reported by Chen et al. (2012) is not detectedin our dust continuum maps, and is probably due to a deconvolution artifact in their maps: IRAM04191 is single at all scales probed with the CALYPSO data (50–3000 au: see Fig. B.6).

Candidate disk

The lowS/N of the binned visibilities from our PdBI dataset for IRAM04191 makes it difficult to robustly model them in detail. However, the envelope intensity radial distribution of IRAM04191 is quite well characterized at scales > 2′′ (see, e.g., Motte & André 2001; Belloche & André 2004), and it is possible to extrapolate the outer envelope properties to establish constraints on the maximum disk-like component (i.e., maximum radius and flux) that can be added while reproducing the dust continuum interferometric fluxes obtained in IRAM04191. We show that the PG model that includes an unresolved Gaussian component (FWHM ~ 24 au) is better than the Pl model at 231 GHz (see Table C.10), while the 94 GHz profile can be well modeled with a Plummer-only model. The flat profile of the 231 GHz dust continuum emission observed at baselines > 300 kλ strongly suggest that an unresolved candidate disk is responsible for the dust continuum emission at scales < 50 au, but our limited S/N cannot provide a robust characterization of the protostellar disk size in IRAM04191 (all error bars on the fluxes overlap in the baseline range 300–600 kλ). Although the current CALYPSO observations do not allow us to firmly establish the size of the disk-like component in IRAM04191, the 231 GHz visibility curve (see Fig. C.9) shows that we should be able to detect an additional disk-like component that would have a half-maximum flux of 4 mJy at 300 kλ: we report this maximum size and flux for the additional disk-like component, for example, < 57 au and flux <4 ± 1 mJy, in Table 5 and report the disk candidate as unresolved at these scales in Table 5.

Fig. C.10 Same as Fig. C.4 for L1521F. Open symbols are used when the visibility real part has a negative value (since absolute values are shown in the log–log plot). The red Plummer + Gaussian (PG) model that includes an unresolved 0.13′′ FWHM Gaussian component (see Table 6) is statistically better than the black Plummer-only (Pl) model at 231 GHz. The 94 GHz profile can be well modeled with a Plummer-only model.

C.9 L1521F

Single-dish constraints

In the IRAM-30 m MAMBO observations by Tokuda et al. (2016) at 240 GHz, the envelope is resolved with an integrated flux 1.0 Jy in a 30′′ radius. The peak flux is ~90 mJy in the 11′′ IRAM-30 m beam (Crapsi et al. 2004). We scaled these fluxes down from the native 243 GHz observing frequency to 94 and 231 GHz using the PdBI spectral index at 20 kλ (see Table 4). We used a 30% uncertainty on the single-dish fluxes, and let the envelope outer radius vary between 20′′ and 40′′.

Multiplicity

In our dust continuum maps, L1521F is a single source (see Fig. B.7). We detect a southeast extension in the 231 GHz map (see Fig. B.7, also called MMS2 in Tokuda et al. 2016), however, this emission cannot be modeled with a compact component from our dust continuum emission visibilities, and it does not have an infrared counterpart. We therefore suggest that MMS2, lying in the outflow cavity, is probably a dust emission feature from a structured cavity.

Candidate disk.

The lowS/N of the binned visibilities from our PdBI dataset for L1521F make it difficult to robustly model them in detail, but it is possible to extrapolate the outer envelope properties to establish constraints on the maximum disk-like component (i.e., maximum radius and flux) that can be added while reproducing the dust continuum interferometric fluxes obtained in L1521F. The 231 GHZ profile is best modeled using a PG model that includes an unresolved Gaussian component (FWHM 0.13′′; see Fig. C.10): this suggests that a candidate disk of radius ~ <20 au and flux ~1.3 ± 0.4 mJy may be detected in L1521F. However, as in the case of IRAM04191, we are limited by the S/N < 3 at long (>300 kλ) baselines forL1521F, so that the size of the Gaussian component that can be added to the Plummer envelope model is subject tohigh uncertainties. To remain robust in our characterization of a potential protostellar disk in L1521F, we report in Table 5 that a candidate disk-like component might contribute a flux 1.3 mJy at scales < 57 au (300 kλ). We also note that 0.87 mm ALMA observations by Tokuda et al. (2017) reported the detection of a small 10 au disk in L1521F, which is consistent with our upper-limit size and flux at 1.3 mm.

Fig. C.11 Same as Fig. C.4 for L1527. The two visibility profiles are better reproduced with models that include an additional Gaussian source whose FWHM is 0.4′′ at 231 GHz.

C.10 L1527

Single-dish constraints

In the IRAM-30 m observations by Motte & André (2001) with MAMBO at 240 GHz, the L1527 envelope is resolved with an integrated flux 1500 mJy in a 4200 au radius. The peak flux is 375 mJy in the 11′′ IRAM-30 m beam. We scaled these fluxes down from 243 to 94 GHz and 231 GHz using the PdBI spectral index at 20 kλ (see Table 4). We used a 20% uncertainty on the extrapolated single-dish fluxes, and let the envelope outer radius vary between 25′′ and 35′′.

Multiplicity

In our dust continuum maps, L1527 is a single source (see Fig. 3).

Candidate disk

The models and the visibility profiles of L1527 are shown in Fig. 7 in the text and in Fig. C.11 here. When all the dust continuum visibilities are used, the best-fit model to reproduce the continuum emission visibility profiles of L1527 is the Plummer + Gaussian model for both frequencies (see main text and Table 6). The best-fit PG model for the 231 GHz visibility profile includes a 215 mJy Gaussian source with FWHM 0.4′′. The parameters for this Gaussian component are reported as the candidate disk properties in Table 5. Similarly, the best-fit PG model to reproduce the 94 GHz visibilities includes a marginally resolved Gaussian component with FWHM = 0.3′′ and flux 23 mJy (see Table 6). When only the equatorial visibilities are used (in a direction orthogonal to the jet axis position angle, i.e., selecting only the visibilities at PA 19° ± 30° for L1527), we find that the Plummer + Gaussian model performs better at reproducing the 231 GHz visibility profile, witha Gaussian size and flux similar to the best-fit model values found using all visibilities (see models Pleqand PGeq in Table C.12). Our analysis of the CALYPSO data therefore suggests that a candidate disk-like structure is detected at radii ~ 50−70 au in L1527.

C.11 Serpens Main S68N

Single-dish constraints

SerpM-S68N (McMullin et al. 1994) is also known as SMM 9 (Casali et al. 1993) or Ser-emb 8 (Enoch et al. 2011), and is located in the Serpens Main cluster. Two protostars are located within the FWHM of the Bolocam core Ser-emb 8 associated with S68N, but the envelope associated with S68N can be clearly identified in the IRAM-30 m observations by Kaas et al. (2004) at better angular resolution. Based on the MAMBO map, we estimate an integrated flux 1030 mJy in a 15′′ radius area, and a peak flux 550 mJy in the 11′′ IRAM-30 m beam. We scaled these fluxes down from the MAMBO central frequency of 243 to 94 GHz and 231 GHz using the PdBI spectral index at 20 kλ (see Table 4). We used a 15% uncertainty on the extrapolated single-dish fluxes, and let the envelope outer radius vary between 8′′ and 22′′.

Multiplicity

SerpM-S68N (S68N in the following) is single in the 231 GHz PdBI primary beam area (21′′ FWHM; see Fig. B.8). Our PdBI 94 GHz map indicates a compact source at the position of S68N, and two additional sources ~12′′ and ~20′′ to the northeast. This is consistent with the CARMA 230 GHz sources found by Enoch et al. (2011). While the region is crowded and it is difficult to robustly associate sources at different wavelengths with different resolutions, it seems that the secondary source S68Nb has no Spitzer MIPS counterpart at 24 μm, while S68Nc has an associated Spitzer IRAC source at 8 μm (Enoch et al. 2009; Harvey et al. 2007): it is classified as a Class I protostar. The three sources (S68N, S68Nb, and S68Nc) are embedded in a common parsec-scale filamentary structure seen in the single-dish maps of the dust continuum emission, but they have separate envelopes.

Candidate disk

When all the dust continuum visibilities are used, the Plummer model is as good a model as the PG model to reproduce the continuum emission visibility profiles of S68N at both frequencies (see main text and Table 6). The best-fit PG model for the 231 GHz visibility profile includes an unresolved Gaussian source with an FWHM 0.11′′ and flux 28 mJy. The parameters for this Gaussian component are reported as the upper-limit values for any disk in S68N in Table 5, although our observations do not hint at the presence of a disk in S68N. When only the equatorial visibilities are used (in a direction orthogonal to the jet axis position angle, i.e., selecting only the visibilities at PA 45° ± 30° for S68N), we find that the Plummer model performs as well as the Plummer + Gaussian model to reproduce the 231 GHz visibility profile, with a Gaussian size and flux similar to the best-fit model values found using all visibilities (see models Pleq and PGeq in Table C.12). Our analysis of the CALYPSO data therefore suggests that no disk-like structure is resolved in S68N at our spatial resolution, and any disk can only be present at radii <50 au.

Fig. C.12 Same as Fig. C.4 for S68N. The Plummer-only models are the best-fit models to reproduce both visibility profiles.

C.12 Serpens Main SMM4

Single-dish constraints

The envelope associated with SerpM-SMM4 (SMM4 in the following) can be clearly identified in the IRAM-30 m observationsby Kaas et al. (2004): based on the MAMBO map, we estimate an integrated flux 2350 mJy in a 20′′ radius area (background subtracted), and a peak flux 1000 mJy in the 11′′ IRAM-30 m beam. We removed the contribution from MM4b that is included in the 11′′ beam, and scaled these fluxes down from the MAMBO central frequency of 243 to 94 GHz and 231 GHz (using the PdBI spectral index at 20 kλ, see Table 4). We used a 20% uncertainty on the extrapolated single-dish fluxes, and let the envelope outer radius vary between 15′′ and 25′′.

Multiplicity

In both our dust continuum maps at 231 and 94 GHz, we detect a secondary component, SMM4b, about 7′′ southwest of the strongest millimeter source (see Fig. B.9). This is the first time that this secondary component is detected, and its nature is therefore unclear. However, it seems that the Herschel/PACS emission peaks toward SMM4b and not toward the main protostar. Moreover, a water maser was detected at (18:29:56.51, 01:13:11.6, equ. J2000) with the VLA, see Furuya et al. (2003): this position is closer to the position of MM4b than to the position of MM4a. Finally, the methanol outflow detected in Kristensen et al. (2010) was not found to peak on the millimeter emission peak reported in the literature. Based on these pieces of evidence, we suggest that SMM4b is probably the driving source of the jet or outflow, while the more quiescent SMM4a dominates the millimeter dust continuum emission. A more detailed analysis of the jet properties and chemical content for these two sources will be proposed byPodio & CALYPSO (in prep.) and Belloche & CALYPSO (in prep.).

Candidate disk

When all the dust continuum visibilities are used, the best-fit model to reproduce the continuum emission visibility profiles of SMM4a is the Plummer + Gauss model for both frequencies (see main text and Table 6). The best-fit PG model for the 231 GHz visibility profile includes a Gaussian source with an FWHM 0.7′′ and flux 595 mJy. The parameters for this Gaussian component are reported as the values for the properties of the candidate disk in SMM4a in Table 5. Hence, our analysis of the CALYPSO data suggests that continuum emission that is not accounted for by circular-symmetric Plummer-like envelope models may trace a disk-like structure in SMM4a that is resolved with our observations with a radius ~ 300 au. If this continuum emission indeed traces a disk structure, the disk flux at 231 GHz accounts for ~ 30% of the total envelope flux from single-dish observations (obtained within a 20′′ radius), which is very unusual for Class 0 protostars. Based on the internal luminosity of SMM4a in the HGBS (2 L_⊙; see Table 1), we would expect a dust temperature ~ 20−30 K at the 200 kλ scale, turning the 231 GHz disk flux into a disk mass >0.3 M_⊙. Such a high disk mass is surprising: a more complete analysis of SMM4a with 2D modeling of interferometric data that cover a larger spatial dynamic range is needed to better understand the nature of this source.

Fig. C.13 Same as Fig. C.4 for SMM4. Open symbols are used when the visibility real part has a negative value (since absolute values are shown in the log-log plot). The PG models are the best-fit models to reproduce the two CALYPSO visibility profiles: at 231 GHz, the additional Gaussian component has an FWHM 0.7′′ and flux 595mJy.

C.13 Serpens South MM18

Single-dish constraints

The envelope associated with Serpens South MM18 (SerpS-MM18 in the following) has been mapped with MAMBO on the IRAM-30 m telescope (Maury et al. 2011): based on the MAMBO map, we estimate an integrated flux 2505 mJy in a 18′′ FWHM source, and a peak flux 1376 mJy in the 11′′ IRAM-30 m beam. We removed the contribution from MM18b and scaled these fluxes down from the MAMBO central frequency of 243 to 94 GHz and 231 GHz (using the PdBI spectral index at 20 kλ; see Table 4). We used a 20% uncertainty on the extrapolated single-dish fluxes, and let the envelope outer radius vary between 10′′ and 20′′.

Multiplicity

SerpS-MM18 is separated into two sources in our CALYPSO maps: MM18a the primary protostar (dominatingthe millimeter continuum emission) and a secondary source MM18b, weaker and more compact, 10′′ to the southwest (see Fig. B.10). At a distance of 250 pc for Serpens South, these two sources are separated by 2600 au. We note that Ortiz-León et al. (2015) argued that Serpens Main and W40 are at a same distance, about 430 pc, based on their VLA parallax measurements of seven sources in both clouds. While W40 and Serpens South seem to belong to the same extinction wall (Straižys et al. 2003), it is not yet clear at which distance the Serpens South filament is located (Könyves et al. 2015). For consistency with previous studies, we use here the value of 250 pc but acknowledge that the distance might be twice larger, and the physical separation between MM18a and MM18b could thus be up to 5000 au. MM18b was originally detected in Maury et al. (2011), then with CARMA by Plunkett et al. (2013) at 3 mm (CARMA-6) and with the VLA by Kern et al. (2016), who classified this source as a Class I (source VLA_13). More recently, it was detected with ALMA in Band 6 (Plunkett et al. 2018; source serps33). MM18a drives a collimated jet (PA + 188° for the blue lobe; Podio et al., in prep.), while MM18b is associated with outflowing gas that follows a cavity with a rather large opening angle (PA of the redshifted emission ~ +8°, from our CALYPSO data).

Candidate disk

When all the dust continuum visibilities are used, the best-fit model to reproduce the 231 GHz continuum emission visibility profile of SerpS-MM18a is the Plummer model that includes a Gaussian (see main text and Table 6), with an FWHM 0.13′′ and flux 76 mJy. At 94 GHz, the Plummer-only model satisfactorily reproduces the visibility profile. Our modeling therefore suggests that a marginally resolved candidate disk might be detected in SerpS-MM18a, and the parameters of the Gaussian component that are included in the best-fit PG model for the 231 GHz profile are reported as the upper-limit disk size and disk flux in Table 5.

C.14 Serpens South MM22

Single-dish constraints

The envelope associated with Serpens South MM22 (SerpS-MM22 in the following) has been mapped with the MAMBO camera on the IRAM-30 m (Maury et al. 2011): based on the MAMBO map, we estimate an integrated flux 261 mJy in a 20′′ FWHM, and a peak flux 129 mJy in the 11′′ IRAM-30 m beam. We scaled these fluxes down from the MAMBO central frequency of 243 to 94 GHz and 231 GHz (using the PdBI spectral index at 20 kλ, see Table 4). We used a 40% uncertainty on the extrapolated single-dish fluxes, and let the envelope outer radius vary between 10′′ and 25′′.

Multiplicity

SerpS-MM22 is single in our dust continuum maps (see Fig. B.11).

Candidate disk

When all the dust continuum visibilities are used, the Pl and PG models are satisfactory for both frequencies (see main text and Table 6), but the Plummer + Gauss models are statistically better (at 231 GHz, the F value is 19, compared to a critical value of 10 above which the probability of a better minimization that is only due to the use of a model containing two more free parameters is <0.3%). The PG model at 231 GHz includes a Gaussian source with an FWHM 0.25′′ and flux 31 mJy. At 94 GHz, the PG model includes a marginally resolved Gaussian component similar to the one found in the 231 GHz visibility profile (FWHM 0.31′′ and flux 3.2 mJy). The parameters of the Gaussian component included in the best-fit model for the 231 GHz profile are reported as properties of the candidate disk in SerpS-MM22, in Table 5. Hence, our analysis suggests that a disk-like structure might be present in SerpS-MM22 at radii ~ 62 au. However, we stress that the low S/N of our CALYPSO visibility profiles for this low-luminosity source precludes us from concluding firmly: additional deeper observations of SerpS-MM22 are needed to fully characterize this candidate disk structure.

Fig. C.14 Same as Fig. C.4 for SerpS-MM18. The PG model is the best-fit model to reproduce the 231 GHz visibility profile: it includes a marginally resolved additional Gaussian component with an FWHM 0.13′′ (flux 76 mJy). The 94 GHz profile can be well modeled with a Plummer-only envelope model.

Fig. C.15 Same as Fig. C.4 for SerpS-MM22. Open symbols are used when the visibility real part has a negative value (since absolute values are shown in the log–logplot). The red Plummer + Gaussian (PG) models are the best-fit models for the two CALYPSO visibility curves. The PG model at 231 GHz includes a 0.25′′ FWHM Gaussian component (see Table 6).

C.15 L1157

Single-dish constraints

Motte & André (2001) IRAM 30 m MAMBO observations at 240 GHz show an unresolved envelope with a peak flux 630 mJy beam⁻¹ in the IRAM-30 m beam 11′′ beam. We used the dual-frequency PdBI spectral index at 20kλ to extrapolate the flux from 240 to 231 GHz. The extrapolated flux is 569 mJy: since the envelope is unresolved in the IRAM-30 m 11′′ beam, we used the IRAM-30 m peak flux as envelope-integrated flux, and constrained the core angular size at 8′′−16′′. The uncertainties on the 231 GHz extrapolated fluxes are ± 30%. For the 94 GHz envelope fluxes, we used an integrated flux extrapolated from the shortest baseline CARMA flux in Kwon et al. (2015), 70 mJy, and a peak flux extrapolated from Motte & André (2001) flux at 240 GHz, 53 mJy at 4kλ, with estimated uncertainties ± 40%.

Fig. C.16 Same as Fig. C.16 for L1157. The red Plummer + Gaussian (PG) model that includes an unresolved Gaussian component (see Table 6 for moreinformations on both models) is statistically better than the black Plummer-only (Pl) model at 231 GHz. At 94 GHz, the best-fit model is the Plummer-envelope model.

Multiplicity

We find no evidence of multiplicity in L1157 down to the 0.4′′ scales probed by our longest baselines. Kwon et al. (2015) suggested that multiple jets originate from the L1157 inner envelope, but the detection of a unique resolved high-velocity jet with CALYPSO observations (Podio et al. 2016) rather suggests L1157 is indeed a single protostar at scales 70−1000 au.

Candidate disk

The models that reproduce the visibility profile for L1157 are shown in Fig. 8. When all the dust continuum visibilities are used, the best-fit model to reproduce the continuum emission visibility profiles of L1157 is the PG model for the 231 GHz data (see main text and Table 6): it includes an unresolved 56 ± 6 mJy Gaussian source. The F-value obtained when the best-fit Pl and best-fit PG model are compared is 9 for the 231 GHz dataset, which is slightly above the critical value (the critical value of F above which the probability of a coincidental better minimization due to two additional parameters is 8). The parameters for this Gaussian component are reported in Table 5. The best-fit model to reproduce the 94 GHz visibilities is the Pl model.

When only the equatorial visibilities are used (in a direction orthogonal to the jet axis position angle, i.e., selecting only the visibilities at PA + 73° ± 30°), we find that the PG model performs slightly better than the Pl model as well (see models Pleq and PGeq in Table C.17). Our analysis of the CALYPSO data therefore suggests that the continuum structure around L1157 is dominated by emission from the inner part of the envelope, and an unresolved candidate disk only contributes at scales < 50 au to the long-baseline 231 GHz fluxes.

C.16 GF9-2

Single-dish constraints

The envelope associated with GF9-2 is quite unconstrained, and this core has been little studied so far. It is also known as L1082C (Bontemps et al. 1996; Caselli et al. 2002a). Located at 200 pc, it has been shown to be associated with an infrared source (Furuya et al. 2006, 2014), and a maser was detected to be associated with the core (Furuya et al. 2001). In Helmut Wiesemeyer’s Ph.D. thesis (Wiesemeyer 1997), where the NH₃ emission maps are used as a temperature probe at an arcminute scale, the envelope mass is estimated to be ~ 0.3 M_⊙. We note that a virial analysis that uses the C¹⁸O line width at similar scales suggests an envelope mass 1 M_⊙, while Furuya et al. (2006) found an envelope mass 0.6 M_⊙ based on the dust continuum maps at 350 μm within a 5400 au area.

Fig. C.17 Same as Fig. C.4 for GF9-2. The 231 GHz visibility profile is shown with black dots and the 94 GHz visibility profile with gray triangles. The red Plummer + Gaussian (PG) model that includes a Gaussian component (0.18′′ FWHM for the 231 GHz and unresolved at 94 GHz) is statisticallybetter than the black Plummer-only (Pl) model. See Tables 6 and C.18 for more information on the models.

We used the IRAM-30 m observations by Wiesemeyer (1997) with the MPIfR bolometer to constrain the large-scale envelope properties. Based on these data, we estimate an integrated flux 315 mJy in an area 20′′ in radius, and a peak flux 60 mJy beam⁻¹ in a 12′′ beam. We scaled these fluxes down from the original central frequency of 243 to 94 GHz and 231 GHz (using the PdBI spectral index at 20 kλ; see Table 4). We used a 40% uncertainty on the extrapolated single-dish fluxes, and let the envelope outer radius vary between 20′′ and 40′′.

Multiplicity

We find no evidence of multiplicity in GF9-2 at the scales probed by our CALYPSO observations (50–5000 au).

Candidate disk

For the two visibility profiles at 231 and 94 GHz, models that include an additional Gaussian component are statistically better at reproducing the intensity radial distribution. Therefore, we conclude that the GF9-2 core seems to include a small-scale disk-like component whose size is ~ 36 ± 10 au and whose flux is 12 mJy at 231 GHz. We stress that the fluxes in the two visibility profiles obtained with the PdBI are surprisingly flat (see Fig. C.17), even at the shortest baselines we probed (10 kλat 94 GHz): no clear envelope-like profile is detected in our data even at the relatively large scales we probe (2−10′′). The envelope emission might either be completely filtered out at these scales, or GF9-2 might be a compact object that is embedded in a large-scale filamentary structure that is detected with the single-dish observations. The OVRO 3 mm continuum emission was found to be shifted to the east with respect to the peak of the N₂H⁺ core, which led Furuya et al. (2006) to propose that the millimeter continuum emission was dominated by a protostellar object while the molecular core traces a younger object that might be prestellar in nature. Here the peaks of our PdBI dust continuum emission coincide with the OVRO 3 mm continuum peak position, which trace the same object. The possibility that GF9-2 is a very young object (e.g., first hydrostatic core, see discussion about this possibility in Furuya et al. 2014) seems to be ruled out by our detection of a jet that originates from the continuum source (Podio et al. in prep.), but the current observations do not allow us to firmly establish the nature of the PdBI dust continuum source we detect.

Our modeling relies on the large-scale envelope parameters that are obtained solely from single-dish observations since emission from the envelope, if any, is not detected in our PdBI data. Because for GF9-2 the quality of the IRAM-30 m data is significantly poorer than in other single-dish data that we used for the remaining sources in the sample, we had to use quite loose constraints on the single-dish fluxes that trace the outer scales in our modeling: we stress that the exact properties of this candidate disk in the GF9-2 core need to be confirmed.

References

All Tables

All Figures

Fig. 1 1.3 mm (231 GHz, left panel) and 3.3 mm (94 GHz, right panel) PdBI dust continuum emission maps of L1448-2A, corrected for primary beam attenuation. The ellipses in the bottom left corner show the respective synthesized beam sizes. The contours are levels of −3σ (dashed), 5σ, and 10σ, then from 20σ in steps of 20σ (rms noise computed in the map before primary beam correction, reported in Table 2). The blue and red arrows show the direction of the protostellar jet(s) either as found in our CALYPSO dataset or from the literature (see Table 1).

In the text

Fig. 2 Same as Fig. 1 for the SVS13 region. The dashed pink contour shows the PdBI primary beam at each frequency. The contours show levels of −3σ (dashed), 5σ, and 10σ, then from 20σ in steps of 20σ (see Table 2). The blue and red arrows show the direction of the protostellar jets from our CALYPSO data for SVS13A/B (Podio et al., in prep.), while the SVS13C outflow PA stems from the CARMA map (Plunkett et al. 2013).

In the text

	Fig. 3 Same as Fig. 1 for L1527. The contours show levels of −3σ, 5σ, and 10σ, and then from 20σ in steps of 40σ (see Table 2).
In the text

	Fig. 4 Examples of radial intensity profiles from single-dish maps obtained at the IRAM-30 m telescope that are used as large-scale constraint for the envelope modeling. These profiles were obtained from the dust continuum emission maps acquired using the IRAM-30 m and are published in Motte & André (2001) and Kaas et al. (2004).
In the text

Fig. 5 Top panels: 1.3 mm (231 GHz) PdBI dust continuum emission visibility real parts as a function of baseline length (circularly averaged in 20 kλ bins) for L1448-2A. The left panel uses a linear scale, while the right panel shows the same data and models in logarithmic scale to enhance the visibility of the long baseline points. The equivalent scales probed by the baselines in the image space are indicated at the top of each panel (computed as 0.6 × λ∕b). In both panels, the plain lines show the best-fit Plummer-only envelope (Pl, black) and Plummer + Gaussian (PG, red) models we found to reproduce the visibility profile of the dust continuum emission in this source. The dashedlines show the two components included in the best-fit PG model: the Plummer-envelope component in dashed lightblue and the additional Gaussian component in dot-dashed dark blue (see Table 6 for more information on the two models). The Plummer + Gaussian (PG) model is not statistically better than the Plummer-only (Pl) for L1448-2A. Bottom panels: same as the top panels for the 94 GHz visibility profiles.

In the text

Fig. 6 Same as Fig. 5 for SVS13B. The red Plummer + Gaussian (PG) model that includes a marginally resolved 0.2′′ FWHM Gaussian component is as good as the black Plummer-only (Pl) model in reproducing the 231 GHz visibility profile, but it provides a more reasonable value of the p + q index. Bottom panels: same as the top panels for the 94 GHz visibility profiles. The PG model does not perform any better than the Pl model for the 94 GHz visibility profile (see Table 6).

In the text

	Fig. 7 Same as Fig. 5 for L1527. The red Plummer + Gaussian (PG) model that includes a 0.4′′ -FWHM Gaussian component (see Table 6 for more information on the two models) is statistically better than the black Plummer-only (Pl) model. Bottom panels: same as the top panels for the 94 GHz visibility profiles.
In the text

	Fig. 8 Same as Fig. 5 for L1157. The black Plummer-only (Pl) model is statistically better than the red Plummer +Gaussian (PG) model that includes an unresolved Gaussian component (see Table 6 for more information on the two models). Bottom panels: same as the top panels for the 94 GHz visibility profiles.
In the text

Fig. 9 Protostellar disk radii Rdisk as a functionof the protostellar bolometric luminosity Lbol for all Class 0 (red circles) and Class I (black stars) protostars observed between 0.7 and 2.7 mm with a spatial resolution better than 50 au from the CALYPSO sample and theliterature. The Class 0 disks shown here are from the CALYPSO sample (filled circles) and those characterized in the literature (contour circles, see text for names and references). When only upper limits on disk sizes could be estimated, the upper-limit values are indicated with arrows pointing down. The dotted red box shows the location of the ten Class 0 disk radii from the 8mm continuum observations of the VANDAM sample (Segura-Cox et al. 2018). The individual Class I disks shown are those described in the text of Sect. 5.1.3. The dashed lines show the median disk radii: in red from the sample of 25 Class 0 disks with $\widetilde{R_0}<$50 au (including upper limit radii), and in black the median radius from the sample of 25 Class I disks with $\widetilde{R_I}=115$ au.

In the text

Fig. B.1 1.3 mm (231 GHz) and 3.3 mm (94 GHz) PdBI dust continuum emission maps of L1448-N. The ellipses in the bottom left corner show the respective synthesized beam sizes. The contours are levels of −3σ (dashed) and 5σ, then in stepsof 10σ from 10σ to 100σ, and finally in steps of 20σ beyond (rms noise computed in the map before primary beam correction as reported in Table 2). The maps have been corrected for primary beam attenuation. The blue and red arrows show the direction of the protostellar jet(s) associated with the millimeter sources in the field, see Table 1 for further details.

In the text

	Fig. B.2 Same as Fig. B.1 for L1448-C. The contours show levels of −3σ (dashed), 5σ, and 10σ and then usesteps of 40σ from 20σ on.
In the text

	Fig. B.3 Same as Fig. B.1 for IRAS2A. The contours show levels of −3σ (dashed), 5σ, and 10σ and then usesteps of 20σ from 20σ on.
In the text

	Fig. B.4 Same as Fig. B.1 for IRAS4A. The contours show levels of −3σ (dashed), 5σ, and 10σ and then usesteps of 20σ from 20σ on. IRAS4B is detected to the southwest in the 94 GHz map, outside the area shown here, so that the structure of IRAS4A is well distinguished in the figure.
In the text

	Fig. B.5 Same as Fig. B.1 for IRAS4B. The contours show levels of −3σ (dashed), 5σ, and 10σ and then usesteps of 30σ from 20σ on. The fluxdensity in the maps has been corrected for primary beam attenuation (the primary beam FWHM is shown as a pink dashed circle in the 94 GHz map).
In the text

	Fig. B.6 Same as Fig. B.1 for IRAM04191. The contours show levels of −3σ and 5σ and then usesteps of 10σ from 10σ on.
In the text

	Fig. B.7 Same as Fig. 1 for L1521F. The contours show levels of −3σ and 5σ and then usesteps of 10σ from 10σ on.
In the text

	Fig. B.8 Same as Fig. B.1 for SerpM-S68N. The contours show levels of −3σ and 5σ and then usesteps of 10σ from 10σ on. The dashed pink contour in the 94 GHz continuum emission map shows the PdBI primary beam.
In the text

	Fig. B.9 Same as Fig. B.1 for SerpM-SMM4. The contours show levels of −3σ, 5σ, and 10σ and then usesteps of 20σ from 20σ on.
In the text

	Fig. B.10 Same as Fig. B.1 for SerpS-MM18. The contours show levels of −3σ and 5σ and then usesteps of 20σ from 10σ on. The dashed pink contour shows the PdBI primary beam.
In the text

	Fig. B.11 Same as Fig. B.1 for SerpS-MM22. The contours show levels of −3σ and 5σ and then usesteps of 10σ from 10σ on.
In the text

	Fig. B.12 Same as Fig. B.1 for L1157. The contours show levels of −3σ and 5σ and then usesteps of 30σ from 20σ on.
In the text

	Fig. B.13 Same as Fig. B.1 for GF9-2. The contours show levels of −3σ, 5σ, and 10σ and then usesteps of 20σ from 20σ on.
In the text

Fig. C.1 Visibility profiles of the 231 GHz (black dots) and 94 GHz (gray triangles) continuum emission from L1448-2A. In each visibility curve, we show the best-fit Plummer-only model (Pl, black curve) and the best-fit Plummer + Gauss (PG, red curve) model. In the case of L1448-2A, the F-test suggests that the PG models are not statistically better than the Plummer-only (Pl) models, although the presence of an unresolved (radius < 0.15′′) disk-like component cannot be robustly excluded.

In the text

	Fig. C.2 1.3 mm (231 GHz) and 3.3 mm (94 GHz) PdBI dust continuum emission maps of L1448-NB1 after subtracting the secondary components. The contours are the same as in Fig. B.1.
In the text

Fig. C.3 Visibility profiles of the 231 GHz (black dots) and 94 GHz (gray triangles) continuum emission from L1448-NB1.Open symbols are used when the visibility real part has a negative value (since absolute values are shown in the log–log plot). In each visibility curve, we show both the best Plummer-only model (Pl, black curve) and the best Plummer + Gauss (PG, red curve) model. In the case of L1448-NB1, the F-test suggests that the PG model for the 231 GHz profile does not perform statistically better than the Plummer-only (Pl) model (see Table 6 for more information on the two models). For the 94 GHz profile, the PG model that includes a 0.94′′ Gaussian component performs better than the Pl model, although none of the two models does satisfactorily reproduce the visibility profile (reduced chi square > 3 in both cases).

In the text

Fig. C.4 Visibility profiles of the 231 GHz (black dots) and 94 GHz (gray triangles) continuum emission from L1448-C. In each visibility curve, we show both the best Plummer-only model (Pl, black curve) and the best Plummer + Gauss (PG, red curve) model. The PG model that includes a marginally resolved 0.16′′ -FWHM Gaussian component reproduces the 231 GHz profile better than the Plummer-only (Pl) model (see Table 6 for more information on the two models). Moreover, the PG model that includes an unresolved additional Gaussian component reproduces the 94 GHz visibility profile better.

In the text

	Fig. C.5 Same as Fig. C.4 for IRAS2A. The Plummer-only model is the best-fit model to reproduce the 231 GHz profile. The PG model, which includes an unresolved additional Gaussian component, performs marginally better than the Pl model for the 94 GHz visibility profile.
In the text

	Fig. C.6 Same as Fig. C.4 for SVS13B. The Plummer-only models are the best-fit models to reproduce the two profiles, although at 231 GHz, a more reasonable (p + q) index is found in the PG model that includes a 0.19′′ FWHM Gaussian component (flux 80 mJy).
In the text

	Fig. C.7 Same as Fig. C.4 for IRAS4A1. The Plummer-only models are the best-fit models to reproduce the two visibility profiles.
In the text

	Fig. C.8 Same as Fig. C.4 for IRAS4B1. The PG model that includes a 0.53′′ Gaussian source (flux 645 mJy) reproduces the 231 GHz visibility profile better, while the 94 GHz visibility profile is satisfactorily reproduced by the Plummer-only model.
In the text

Fig. C.9 Same as Fig. C.4 for IRAM04191. Open symbols are used when the visibility real part has a negative value (since absolute values are shown in the log-log plot). The red Plummer + Gaussian (PG) model that includes an unresolved 0.17′′ -FWHM Gaussian component (see Table 6) is only statistically better than the black Plummer-only (Pl) model at 231 GHz. The 94 GHz profile can be well modeled with a Plummer-only model. The oscillations of the models are due to the Hankel transform of the power-law envelope model with a low index (p + q = 1.4).

In the text

Fig. C.10 Same as Fig. C.4 for L1521F. Open symbols are used when the visibility real part has a negative value (since absolute values are shown in the log–log plot). The red Plummer + Gaussian (PG) model that includes an unresolved 0.13′′ FWHM Gaussian component (see Table 6) is statistically better than the black Plummer-only (Pl) model at 231 GHz. The 94 GHz profile can be well modeled with a Plummer-only model.

In the text

	Fig. C.11 Same as Fig. C.4 for L1527. The two visibility profiles are better reproduced with models that include an additional Gaussian source whose FWHM is 0.4′′ at 231 GHz.
In the text

	Fig. C.12 Same as Fig. C.4 for S68N. The Plummer-only models are the best-fit models to reproduce both visibility profiles.
In the text

	Fig. C.13 Same as Fig. C.4 for SMM4. Open symbols are used when the visibility real part has a negative value (since absolute values are shown in the log-log plot). The PG models are the best-fit models to reproduce the two CALYPSO visibility profiles: at 231 GHz, the additional Gaussian component has an FWHM 0.7′′ and flux 595mJy.
In the text

	Fig. C.14 Same as Fig. C.4 for SerpS-MM18. The PG model is the best-fit model to reproduce the 231 GHz visibility profile: it includes a marginally resolved additional Gaussian component with an FWHM 0.13′′ (flux 76 mJy). The 94 GHz profile can be well modeled with a Plummer-only envelope model.
In the text

	Fig. C.15 Same as Fig. C.4 for SerpS-MM22. Open symbols are used when the visibility real part has a negative value (since absolute values are shown in the log–logplot). The red Plummer + Gaussian (PG) models are the best-fit models for the two CALYPSO visibility curves. The PG model at 231 GHz includes a 0.25′′ FWHM Gaussian component (see Table 6).
In the text

	Fig. C.16 Same as Fig. C.16 for L1157. The red Plummer + Gaussian (PG) model that includes an unresolved Gaussian component (see Table 6 for moreinformations on both models) is statistically better than the black Plummer-only (Pl) model at 231 GHz. At 94 GHz, the best-fit model is the Plummer-envelope model.
In the text

	Fig. C.17 Same as Fig. C.4 for GF9-2. The 231 GHz visibility profile is shown with black dots and the 94 GHz visibility profile with gray triangles. The red Plummer + Gaussian (PG) model that includes a Gaussian component (0.18′′ FWHM for the 231 GHz and unresolved at 94 GHz) is statisticallybetter than the black Plummer-only (Pl) model. See Tables 6 and C.18 for more information on the models.
In the text