© M. Gaudel et al. 2020

1 Introduction

Stars form via the gravitational collapse of ≲0.1 pc dense cores, which are embedded within molecular clouds (André et al. 2000; Ward-Thompson et al. 2007; di Francesco et al. 2007). Prestellar cores become unstable and collapse due to their own gravitational potential. One or several stellar embryos form in their center. This is the beginning of the main accretion phase called the protostellar phase. Observations of the molecular line emission from large samples of cores in close star-forming regions revealed that velocity gradients are ubiquitous to prestellar structures at scales of 0.1−0.5 pc (Goodman et al. 1993; Caselli et al. 2002). These were interpreted as slow rotation inherited from their formation process (Goodman et al. 1993; Caselli et al. 2002; Ohashi et al. 1999; Redman et al. 2004; Williams et al. 2006). Assuming these gradients trace organized rotational motions, the observed velocities lead to a typical angular rotation velocity of Ω ~ 2 km s⁻¹ pc⁻¹ and specific angular momentum values of j = v× r ~ 10⁻³−10⁻¹ km s⁻¹ pc.

During the collapse, if the angular momentum of the parent prestellar cores is totally transferred to the stellar embryo, the gravitational force can not counteract the centrifugal force and the stellar embryo fragments prematurely before reaching themain sequence. This is the angular momentum problem for star formation (Bodenheimer 1995). Although observational studies suggest a trend of decreasing specific angular momentum toward smaller core sizes of j ∝ r^1.6, the j measured in prestellar structures at scales of 10 000 au (~ 10⁻³ km s⁻¹ pc, Caselli et al. 2002)is still typically three orders of magnitude higher than the one associated with the maximum rotational energy that a solar-type star can sustain (j_break ~ 10¹⁸ cm² s⁻¹ ~ 3 ×10⁻⁶ km s⁻¹ pc). The physical mechanisms responsible for the angular momentum redistribution before the matter is accreted by the central stellar object have still to be identified.

During the star formation process, disk formation is expected to be a consequence of angular momentum conservationduring the collapse of rotating cores (Cassen & Moosman 1981; Terebey et al. 1984). From observational studies, disks are common in Class II objects (Andrews et al. 2009; Isella et al. 2009; Ricci et al. 2010; Spezzi et al. 2013; Piétu et al. 2014; Cieza et al. 2019). Thus, disk formation has been naturally considered as a possible solution to the angular momentum problem by redistributing the four orders of magnitude of j measured from prestellar cores to the T-Tauri stars (j ~ 2 ×10⁻⁷ km s⁻¹ pc, Bouvier et al. 1993): the disk would store and evacuate the angular momentum of the matter by viscous friction (Lynden-Bell & Pringle 1974; Hartmann et al. 1998; Najita & Bergin 2018) or thanks to disk winds (Blandford & Payne 1982; Pelletier & Pudritz 1992; Pudritz et al. 2007) before the matter is accreted by the central stellar object. However, the spatial distribution of angular momentum during disk formation within star-forming structures at scales between the outer core radius and the stellar surface are still largely unconstrained.

Class 0 protostars are the first (proto)stellar objects observed after the collapse in prestellar cores (André et al. 1993, 2000). Due to their youth, most of their mass is still in the form of a dense, collapsing, reservoir envelope surrounding the central stellar embryo (M_env ≫ M_⋆). Thus, they are likely to retain the initial conditions inherited from prestellar cores, in particular regarding angular momentum. The young stellar embryo mass increases via the accretion of the gaseous and dusty envelope in a short timescale (t < 10⁵ yr, Evans et al. 2009; Maury et al. 2011). During this main accretion phase, most of the final stellar mass is accreted and, at the same time, the infalling gas must redistribute most of its initial angular momentum before reaching the central stellar embryo. Class 0 protostars are therefore key objects to understand the distribution of angular momentum of the material directly involved in the star formation process and constrain physical mechanisms responsible for the redistribution as disk formation.

Clear signatures of rotation (Belloche et al. 2002; Belloche & André 2004; Chen et al. 2007) and infalling gas are generally detected in the envelopes of Class 0 protostars (see the review by Ward-Thompson et al. 2007). Thanks to observations of the dense molecular gas emission, rotational motions where characterized in seven Class 0 or I protostellar envelopes at scales between 3500 and 10000 au (Ohashi et al. 1997a; Belloche et al. 2002; Chen et al. 2007). These envelopes exhibit an average angular momentum of ~10⁻³ km s⁻¹ pc at scales of r < 5000 au, consistent with the j measured in prestellar cores by Caselli et al. (2002) (~10⁻³ km s⁻¹ pc). These studies suggest an angular momentum that is constant with radius in Class 0 protostellar envelopes. These flat profiles are generally interpreted as the conservation of the angular momentum. From hydrodynamical simulations, the conservation of these typical values of angular momentum results in the formation of large rotationally supported disks with radii >100 au in a few thousand years (Yorke & Bodenheimer 1999). However, from observational studies, large disks are rare around Class 0 protostars (r < 100 au, Maury et al. 2010, 2019; Kurono et al. 2013; Yen et al. 2013; Segura-Cox et al. 2016). Only recent numerical simulations including magnetohydrodynamics and non-ideal effects, such as ambipolar diffusion, Ohmic dissipation, or the Hall effect, allow small rotationally supported disks to be formed (< 100 au; Machida et al. 2014; Tsukamoto et al. 2015; Masson et al. 2016).

Very few studies have been able to produce resolved profiles of angular momentum to characterize the actual amount of angular momentum present at the smallest scales (r ≲1000 au) within protostellar envelopes. From interferometric observations, Yen et al. (2015a) derive specific angular momentum values of ~ 2 × 10⁻⁴ km s⁻¹ pc in seven Class 0 envelopes at scales of r ~ 1000 au, values which are below the trend observed by Ohashi et al. (1997a). These studies have put constraints on the angular momentum properties of Class 0 protostellar envelopes and suggest the material at r ~ 1000 au must reduce its angular momentum by at least one order of magnitude from outer envelope to disk scales. Yen et al. (2011, 2017) show specific angular momentum profiles down to ~350 au in two Class 0 protostellar envelopes: j ~ 6 ×10⁻⁴ km s⁻¹ pc at r ~ 1000 au and j ≲10⁻⁴ km s⁻¹ pc at r ~ 350 au. In this case, conservation of angular momentum during rotating protostellar collapse might not be the dominant process leading to the formation of disks and stellar multiple systems. It is therefore crucial to obtain robust estimates of the angular momentum of the infalling material in protostellar envelopes during the main accretion phase by analyzing the kinematics from the outer regions of the envelope (10 000 au, Motte & André 2001) to the protostellar disk (< 50 au, Maury et al. 2019).

2 The CALYPSO survey

The Continuum And Lines in Young ProtoStellar Objects (CALYPSO¹) IRAM Large Program is a survey of 16 nearby Class 0 protostars (d < 450 pc), carried out with the IRAM Plateau de Bure interferometer (PdBI) and IRAM 30-meter telescope (30m) at wavelengths of 1.29, 1.37, and 3.18 mm. The CALYPSO sources are among the youngest known solar-type Class 0 objects (André et al. 2000) with envelope masses of M_env ~1.5 M_⊙ and internal luminosities of L_bol ~0.1−30 L_⊙ (Maury et al. 2019).

The CALYPSO program allows us to study in detail the Class 0 envelope chemistry (Maury et al. 2014; Anderl et al. 2016; De Simone et al. 2017; Belloche et al. 2020), disk properties (Maret et al. 2014, 2020; Maury et al. 2019) and protostellar jets (Codella et al. 2014a; Santangelo et al. 2015; Podio et al. 2016, and in prep.; Lefèvre et al. 2017; Anderl et al. 2020). One of the main goals of this large observing program is to understand how the circumstellar envelope is accreted onto the central protostellar object during the Class 0 phase, and ultimately tackle the angular momentum problem of star formation. This paper presents an analysis of envelope kinematics, for the 12 sources from the CALYPSO sample located at d ≤ 350 pc (see Table 1) and discuss our results on the properties of the angular momentum in Class 0 protostellar envelopes.

We adopt the dust continuum peak at 1.3 mm (225 GHz) determined from the PdBI datasets by Maury et al. (2019) as origin of the coordinate offsets of the protostellar envelopes (see Table 1). We report for each source in Table 1 the outflow axis considered as the rotation axis and estimated by Podio & CALYPSO (in prep.) from high-velocity emission of ¹²CO, SiO, and SO at scales <10′′. We assume the equatorial axis of the protostellar envelopes, namely the intersection of the equatorial plane with the plane of the sky at the distance of the source, to be perpendicular to the rotation axis. The SiO emission in the CALYPSO maps is very collimated, so the uncertainties on the direction of the rotation axis, and thus on the direction of the equatorial axis, are smaller than ±10° (Podio & CALYPSO, in prep.). For L1521-F, IRAM04191, and GF9-2, no collimated SiO jet is detected, thus, the uncertainties are a bit larger (±20°). We also use estimates of the inclination of the equatorial plane with respect to the line of sight from the literature. These estimates, which come from geometric models that best reproduce the outflow kinematics observed in molecular emission, are highly uncertain since we do not have access to the 3D-structure of each source.

3 Observations and dataset reduction

To probe the dense gas in our sample of protostellar envelopes, we use high spectral resolution observations of the emission of two molecular lines, C¹⁸O (2−1) at 219.560 GHz and N₂H⁺ (1−0) at 93.171 GHz. In this section, we describe the dataset properties exploited to characterize the kinematics of the envelopes at radii between r ~50 and 5000 au from the central object.

3.1 Observations with the IRAM Plateau de Bure Interferometer

Observations of the 12 protostellar envelopes considered here were carried out with the IRAM Plateau de Bure Interferometer (PdBI) between September 2010 and March 2013. We used the 6-antenna array in two configurations (A and C), providing baselines ranging from 16 to 760 m, to carry out observations of the dust continuum emission and a dozen molecular lines, using three spectral setups (around 94, 219, and 231 GHz). Gain and flux were calibrated using CLIC which is part of the GILDAS² software. The details of CALYPSO observations and the calibration carried out are presented in Maury et al. (2019). The phase self-calibration corrections derived from the continuum emission gain curves, described in Maury et al. (2019), were also applied to the line visibility dataset (for all sources in the restricted sample studied hereexcept the faintest sources IRAM04191, L1521F, GF9-2, and L1448-2A). Here, we focus on the C¹⁸O (2−1) emission line at 219 560.3190 MHz and the N₂H⁺ (1−0) emission line at 93 176.2595 MHz, observed with high spectral resolution (39 kHz channels, i.e., a spectral resolution of 0.05 km s⁻¹ at 1.3 mm and 0.13 km s⁻¹ at 3 mm). The C¹⁸O (2−1) maps were produced from the continuum-subtracted visibility tables using either (i) a robust weighting of 1 for the brightest sources to minimize the side-lobes, or (ii) a natural weighting for the faintest sources (IRAM04191, L1521F, GF9-2, and L1448-2A) to minimize the rms noise values. We resampled the spectral resolution to 0.2 km s⁻¹ to improve the signal-to-noise ratio (S/N) of compact emission. The N₂H⁺ (1−0) maps were produced from the continuum-subtracted visibility tables using a natural weighting for all sources. In all cases, deconvolution was carried out using the Hogbom algorithm in the MAPPING program of the GILDAS software.

3.2 Short-spacing observations from the IRAM 30-m telescope

The short-spacing observations were obtained at the IRAM 30-m telescope (30m) between November 2011 and November 2014. Details of the observations for each source are reported in Table A.1. We observed the C¹⁸O (2−1) and N₂H⁺ (1−0) lines using the heterodyne Eight MIxer Receiver (EMIR) in two atmospheric windows: E230 band at 1.3 mm and E090 band at 3 mm (Carter et al. 2012). The Fast Fourier Transform Spectrometer (FTS) and the VErsatile SPectrometer Array (VESPA) were connected to the EMIR receiver in both cases. The FTS200 backend provided a large bandwidth (4 GHz) with a spectral resolution of 200 kHz (0.27 km s⁻¹) for the C¹⁸O line, while VESPA provided high spectral resolution observations (20 kHz channel or 0.063 km s⁻¹) of the N₂H⁺ line. We used the on-the-fly spectral line mapping, with the telescope beam moving at a constant angular velocity to sample regularly the region of interest (1′ × 1′ coverage for the C¹⁸O emission and 2′ × 2′ coverage for the N₂H⁺ emission). The mean atmospheric opacity at 225 GHz was τ₂₂₅ ~ 0.2 during the observations at 1.3 mm and τ₂₂₅ ~ 0.5 during the observations at 3 mm. The mean values of atmospheric opacity are reported for each source in Table A.1. The telescope pointing was checked every 2−3 h on quasars close to the CALYPSO sources, and the telescope focus was corrected every 4−5 h using the planets available in the sky. The single-dish dataset were reduced using the MIRA and CLASS programs of the GILDAS software following the standard steps: flagging of incorrect channels, temperature calibration, baseline subtraction, and gridding of individual spectra to produce regularly-sampled maps.

3.3 Combination of the PdBI and 30m data

The IRAM PdBI observations are mostly sensitive to compact emission from the inner envelope. Inversely, single-dish dataset contains information at envelope scales (r ~ 5−40′′) but its angular resolution does not allow us to characterize the inner envelope emission at scales smaller than the beamwidth. To constrain the kinematics at all relevant scales of the envelope, one has to build high angular resolution dataset which recovers all emission of protostellar envelopes. We merged the PdBI and the 30m datasets (hereafter PdBI + 30m) for each tracer using the pseudo-visibility method³ : we generated pseudo-visibilities from the Fourier transformed 30m image data, which are then merged to the PdBI dataset in the MAPPING program of the GILDAS software. This process degrades the angular resolution of the PdBI dataset but recovers a large fraction of the extended emission. The spectral resolution of the combined PdBI + 30m dataset is limited by the 30m dataset at 1.3 mm (0.27 km s⁻¹) and the PdBI one at 3 mm (0.13 km s⁻¹).

We produced the N₂H⁺ (1−0) combined datacubes in such a way to have a synthesized beam size < 2′′ and a noise level < 10 mJy beam⁻¹. As the N₂H⁺ emission traces preferably the outer protostellar envelope, we used a natural weighting to build the combined maps to minimize the noise rather than to maximize angular resolution. The C¹⁸O (2−1) combined maps were produced using a robust weighting scheme, in order to obtain synthesized beam sizes close to the PdBI ones, and to minimize the side-lobes. In all cases, the deconvolution was carried out using the Hogbom algorithm in MAPPING.

3.4 Properties of the analyzed maps

Following the procedure described above, we have obtained, for each source of the sample, a set of three cubes for each of the two molecular tracers C¹⁸O (2−1) and N₂H⁺ (1−0), probing the emission at different spatial scales (PdBI map, combined PdBI+30m map, and 30m map). In order to build maps with pixels that contain independent dataset and avoid oversampling, we inversed visibilities from the PdBI and the combined PdBI+30m datasets using only 4 pixels per synthesized beam, and we smoothed the resulting maps afterwards to obtain 2 pixels per element of resolution. The properties of the resulting maps are reported in Appendix B. The spatial resolution of the molecular line emission maps is reported in Tables B.1 and B.2. The spatial extent of the molecular emission, the rms noise levels, and the integrated fluxes are reported in Tables B.3 and B.4.

4 Envelope kinematics from high dynamic range datasets

4.1 Integrated intensity maps

To identify at which scales of the protostellar envelopes the different datasets are sensitive to, we produced integrated intensity maps by integrating spectra of each pixel for the molecular lines C¹⁸O (2−1) and N₂H⁺ (1−0) from the PdBI, combined, and 30m datasets for each source. For C¹⁸O (2−1), we integrated each spectrum on a velocity range of ± 2.5 km s⁻¹ around the velocity of the peak of the mean spectrum of each source. The 1−0 line of N₂H⁺ has a hyperfine structure with seven components (see Fig. 2). We integrated the N₂H⁺ spectra over a range of 20 km s⁻¹ encompassing the seven components. Figure 1 shows as an example the integrated intensity maps obtained forL1448-2A. The integrated intensity maps of the other sources are provided in Appendix H.

We used the integrated intensity maps to measure the average emission size of each tracer in each dataset above a 5σ threshold. The values reported in Tables B.3 and B.4 are the average of two measurements: an intensity cut along the equatorial axis and circular averages at different radii around the intensity peak position of the source. Only pixels whose intensity is at least 5 times higher than the noise in the map are considered to build these intensity profiles. The FWHM of the adjustment by a Gaussian function allows us to determine the average emission size of the sources. For both tracers and for all sources in our sample, the emission is detected above 5σ in an area larger in the combined datasets than in the PdBI datasets, and smaller than in the 30m ones (see Tables B.3 and B.4). Our three datasets are thus not sensitive to the same scales and allow us to probe different scales within the 12 sampled protostellar envelopes: the 30m datasets trace the outer envelope, the PdBI datasets the inner part and the combined ones the intermediate scales.

The C¹⁸O and N₂H⁺ molecules do not trace the same regions of the protostellar envelope either: Anderl et al. (2016) report from an analysis of the CALYPSO survey that the N₂H⁺ emission forms a ring around the central C18O emission in four sources. Previous studies (Bergin et al. 2002; Maret et al. 2002, 2007; Anderl et al. 2016) show that N₂H⁺, which is abundant in the outer envelope, is chemically destroyed when the temperature in the envelope reaches the critical temperature (T ≳ 20K) at which CO desorbs from dust ice mantles. Thus, while N₂H⁺ can be used to probe the envelope kinematics at outer envelope scales, C¹⁸O can be used as a complementary tracer of the gas kinematics at smaller radii where the embedded protostellar embryo heats the gas to higher temperatures.

The C¹⁸O emission is robustly detected (>5σ) in our PdBI observations for most sources, except for L1521F and IRAM04191 which are the lowest luminosity sources of our sample (see Table 1), and for SVS13-B where the emission is dominated by its companion, the Class I protostar SVS13-A. For most sources, the interferometric map obtained with the PdBI shows mostly compact emission (r <3′′, see Table B.3). However, the C¹⁸O emission from the 30m datasets shows more complex structures (see Appendix H). Assuming that, under the hypothesis of spherical geometry, the emission from a protostellar envelope is compact (r ≲ 40′′, i.e., ≲10000 au, see Table 1) and stands out from the environment in which it is embedded, the 30m emission of L1448-2A, L1448-C, and IRAS4A comes mainly from the envelope.

The N₂H⁺ emission is detected in our combined observations for all sources. In the four sources studied by Anderl et al. (2016), they do not detect the emission at the 1.3 mm continuum peak, but emission rings around the C¹⁸O central emission. From Table B.4, we noticed two types of emission morphologies based on the PdBI dataset: compact (r < 7′′, see Table B.4) or filamentary (r ≥ 9′′). In the same way as the C¹⁸O emission, the N₂H⁺ emission from the 30m datasets shows complex structures with radius r ≳ 40′′ for most sources, except for five sources (IRAM04191, L1521F, L1448-NB, L1448-C, and L1157) where the emission is consistent with the compact emission of the protostellar envelope.

The C¹⁸O emission from the PdBI is not centered on the continuum peak for three sources in our sample: IRAS4A, L1448-NB, and L1448-2A (see Appendix H). For each of these sources, the PdBI 1.3 mm dust continuum emission map resolves a close binary system (< 600 au) with both components embedded in the same protostellar envelope (Maury et al. 2019, see Table 1). The origin of the coordinate offsets is chosen to be the main protostar, secondary protostar, and the middle of the binary system for IRAS4A, L1448-NB, and L1448-2A, respectively, to study the kinematics in a symmetrical way.

4.2 Velocity gradients in protostellar envelopes

To quantify centroid velocity variations at all scales of the protostellar envelopes, we produced centroid velocity maps of each Class 0 protostellar envelope by fitting all individual spectra (pixel by pixel) by line profile models in the CLASS program of the GILDAS software. We only considered the line intensity detected with a S/N higher than 5. We fit the spectra to be able to deal withmultiple velocity components. Indeed, because protostellar envelopes are embedded in large-scale clouds, multiple velocity components can be expected on some lines of sight where both the protostellar envelope and the cloud emit. For example, Belloche et al. (2006) find several velocity components in their 30m of the N₂H⁺ emission of IRAS4A (see Appendix H.6). Except for IRAS4A and IRAS4B for which we fit two velocity components (see details in Appendix H), for most sources we used a Gaussian line profile to model the C¹⁸O (2−1) emission, with the line intensity, full width at half maximum (FWHM), and centroid velocity let as free parameters (see Fig. 2). In the case of N₂H⁺ (1−0), we used a hyperfine structure (HFS) line profile to determine the FWHM and centroid velocity of the molecular line emission (see Fig. 2). Figures 3–14 show the centroid velocity maps obtained for each source of the sample using the PdBI, combined, and 30m datasets for both the C¹⁸O and N₂H⁺ emission.

For most sources in our sample, these centroid velocity maps reveal organized velocity patterns with blue-shifted and red-shiftedvelocity components on both sides of the central stellar embryo, along the equatorial axis where such velocity gradients could be due to rotation of the envelopes. The global kinematics in Class 0 envelopes is a complex combinationof rotation, infall, and outflow motions. The observed velocities are projected on the line of sight and thus, are a mix of the various gas motions. Therefore, it is not straightforward to interpret a velocity gradient in terms of the underlying physical process producing it. In order to have an indication of the origin of these gradients, we performed a least-square minimization of a linear velocity gradient model on the velocity maps following: $$v_{\text{grad}}=v_0+a\Delta \alpha +b\Delta \beta ,$$(1)

with Δα and Δβ the offsets with respect to the central source (Goodman et al. 1993).

This simple model provides an estimate of the reference velocity v₀ called systemic velocity, the direction Θ, and the amplitude G of the mean velocity gradient. One would expect a mean gradient perpendicular to the outflow axis if the velocity gradient was due to rotational motions in an axisymmetric envelope. A mean gradient oriented along the outflow axis could be due to jets and outflows or infall in a flattened geometry. The gradients were fit on the region of the velocity maps shown in Fig. 3, namely 10′′ × 10′′ in the PdBI and combined datasets for the C¹⁸O emission (lower left and central panels), 40′′ × 40′′ for the N₂H⁺ emission from the PdBI and combined datasets (upper left and central panels), and 80′′ × 80′′ and 160′′ × 160′′, respectively for the C¹⁸O and N₂H⁺ emission from the 30m datasets (right panels). Table 2 reports for each source the significant mean velocity gradients detected with an amplitude higher than 2σ. No significant velocity gradient is observed for IRAM04191, L1521F, and SVS13-B in C¹⁸O emission at scales of r <5′′ or for L1448-C and IRAS4B at r > 30′′ (see Table 2).

Seven of the 12 sources in our sample show a mean gradient in C¹⁸O emission aligned with the equatorial axis (ΔΘ < 30°) which could trace rotational motions of the envelope at scales of r < 5′′. At similar scales, four sources (L1448-NB, L1521F, L1157, GF9-2) show gradients with intermediate orientation (30° < ΔΘ < 60°). Finally, L1448-2A shows a mean gradient aligned to the outflow axis rather than the equatorial axis (ΔΘ > 60°). For these last five sources, the gradients observed could be due to a combination of rotation, ejection, and infall motions. For all sources, we noticed a systematic dispersion of the direction of the velocity gradient from inner to outer scales in the envelope (see Fig. 19). We discuss in Sect. 5.4 whether this shift in direction of the velocity gradient is due to a transition from rotation-dominated inner envelope to collapse-dominated outer envelope at r > 1500 au, or is due to the different molecular tracers used for this analysis. In most sources, the gradient moves away from the equatorial axis as the scale increases. Only three sources (IRAM04191, L1521F, and L1527) show a gradient close to the equatorial axis with ΔΘ < 30° at 2000 au in N₂H⁺ emission from the combined dataset while four sources show a complex gradient and five sources have a ΔΘ > 60°.

4.3 High dynamic range position-velocity diagrams to probe rotational motions

To investigate rotational motions and characterize the angular momentum properties in our sample of Class 0 protostellar envelopes, we build the position-velocity (PV_rot) diagrams along the equatorial axis. We assumed the position angle of the equatorial axis as orthogonal to the jet axis reported in Table 1. The choice of this equatorial axis allows us to maximize sensitivity to rotational motions and minimize potential contamination on the line of sight due to collapsing or outflowing gas (Yen et al. 2013). The velocities reported in the PV_rot diagram are corrected for the inclination i of the equatorial plane with respect to the line of sight (see Table 1). We note that the correction for inclination is a multiplicative factor, thus if this inclination angle is not correctly estimated, the global observed shape is not distorted.

The analysis described in detail in Appendix C allows us to build a PV_rot diagram with a high dynamic range from 50 au up to 5000 au for each source as follows (see the example of L1527 in Fig. 15):

Firstly, to constrain the PV_rot diagram at the smallest scales resolved by our dataset (~0.5′′), we use the PdBI C¹⁸O datasets that we analyze in the (u,v) plane to avoid imaging and deconvolution processes (see label “C¹⁸O PdBI” in Fig. 15). We only kept central emission positions in the channel maps at a position angle < |45°| with respect to the equatorial axis (see Appendix C).

Secondly, since the C¹⁸O extended emission is filtered out by the interferometer, we used the combined C¹⁸O emission to populate the PV_rot diagram at the intermediate scales of the protostellar envelopes (see label “C¹⁸O combined” in Fig. 15). The C¹⁸O molecule remains the most precise tracer when the temperature is higher than ~20 K because below, the C¹⁸O molecule freezes onto dust ice mantles. To determine the transition radius R_trans between the two tracers, we calculate the C¹⁸O and N₂H⁺ column densities along the equatorial axis from the combined integrated intensity maps (see Appendix D and green points in Fig. 15).

Then, at radii r > R_trans, the N₂H⁺ emission traces better the envelope dense gas. We use the combined N₂H⁺ emission maps to analyze the envelope kinematics at larger intermediate scales. When the N₂H⁺ column density profile reaches a minimum value due to the sensitivity of the combined dataset, this dataset is no longer the better dataset to provide a robust information on the velocity (see label “N₂H⁺ combined” in Fig. 15).

Finally, we use the 30m N₂H⁺ emission map to populate the PV_rot diagram at the largest scales of the envelope (see label “N₂H⁺ 30m” in Fig. 15).

The CALYPSO datasets allow us to continuously estimate the velocity variations along the equatorial axis in the envelope over scales from 50 au up to 5000 au homogeneously for each protostar. Figure 16 shows the PV_rot diagrams built for all sources of the sample. The systemic velocity used in the PV_rot diagrams are determined in Appendix E.

The method of building PV_rot diagrams described above and in Appendix C corresponds to an ideal case with a detection of a continuous blue-red velocity gradient along the direction perpendicular to the outflow axis in the velocity maps. In practice, the direction of velocity gradients is not always continuous at all scales probed by our observations(see Table 2 and Fig. 19). For some sources, to constrain the PV_rot diagrams, we did not take the kinematic information at all scales of the envelope into account. Velocity gradients can be considered as probing rotational motions if the following criteria are met:

Firstly, we only consider the significant velocity gradients reported in Table 2 with a blue- and a red-shifted velocity components observed on each side of the protostellar embryo, itself at the systemic velocity v₀. For example, we only take the C¹⁸O emission from the 30m map into account for L1521F (see Fig. 11).

Secondly, we only take the velocity gradients aligned with the equatorial axis (ΔΘ < 60°) into accountin order to minimize contamination by infall and ejection motions. For example, we do not report in the PV_rot diagrams the N₂H⁺ velocity gradients from the combined maps for L1448-NB, IRAS2A, SVS13-B, and GF9-2 (see Table 2).

Finally, we do not consider the discontinuous velocity gradients which show an inversion of the blue- and red-shifted velocity components along the equatorial axis from inner to outer envelope scales. For example, we do not report in the PV_rot diagrams the N₂H⁺ velocity gradients from the 30m maps for IRAM04191 and L1157 (see Figs. 10 and 13).

When velocity gradients with a blue- and a red-shifted velocity components observed on each side of the protostellar embryo are continuous from inner to outer envelope scales but shifted from the equatorial axis (ΔΘ ≥ 60°), we only report upper limits on rotational velocities in the PV_rot diagrams. The sources in our sample show specific individual behaviors, therefore we followed as closely as possible the method of building the PV_rot diagram adapting it on a case-by-case basis.

5 Discussion

In this section, we discuss the presence of rotation in the protostellar envelopes from the PV_rot diagrams (see Sect. 5.1 and Fig. 16). We build the distribution of specific angular momentum associated with the PV_rot diagrams (see Sect. 5.2) and explore the possible solutions to explain the j(r) profiles observed in the inner (r < 1600 au, see Sect. 5.3) and outer (r > 1600 au, see Sect. 5.4) parts of the envelopes.

5.1 Characterization of rotational motions

We assume that the protostellar envelopes are axisymmetric around their rotation axes, and thus, the velocity gradients observed along the equatorial axis, and reported in the PV_rot diagrams, are mostly due to the rotational motions of the envelopes. We model the rotational velocity variations by a simple power-law model v ∝ r^α without taking the upper limits into account. This method has been tested with an axisymmetric model of collapsing-rotating envelopes by Yen et al. (2013). As long as rotation dominates the velocity field on the line of sight, which depends on the inclination and flattening of the envelope, Yen et al. (2013) obtain robust estimates of the rotation motions at work in the envelopes. First, we fix the power-law index at α = −1 to compare to what is theoretically expected for an infalling and rotating envelope from a progenitor core in solid-body rotation (Ulrich 1976; Cassen & Moosman 1981; Terebey et al. 1984; Basu 1998). The reduced χ² values of fits by an orthogonal least-square model are reported in the second column of Table 3. Then, we let the power-law index vary as a free parameter: the best power-law index and the reduced χ² found for each protostellar envelope in our sample are reported in the third column of Table 3. Figure 16 shows the PV_rot diagrams adjusted by a power-law for the sources of the CALYPSO sample.

The power-law indices of the PV_rot diagrams from our sample are between −1.1 and 0.8. Five sources (L1448-2A, IRAS2A, SVS13-B, L1527, and GF9-2) show rotational velocity variations in the envelope scaling as a power law with an index close to −1. This is consistent with the expected index for collapsing and rotating protostellar envelopes. The reduced χ² are ~1.5 for these sources except for IRAS2A and SVS13-B for which it is better (~0.2). L1521F and L1157 show a power-law index close to 0 with a very low reduced χ² (≤0.3, see Table 3). These flat PV_rot diagrams (v_rot ~ constant) suggest differential rotation of the envelope with an angular velocity of $\Omega =\frac{v_{\text{rot}}}{r}\propto r^{-1}$. For two other sources (IRAS4B and IRAM04191), the best indices are compatible with −0.5, which could suggest Keplerian rotation at scales of r < 1000 au. However, the reduced χ² are also satisfactory (~1) when we fix the power-law index at α = −1 (see Table 3). Thus, for these two sources, our CALYPSO datasets only allow us to estimate a range of the power-law indices between −1 and −0.5 (see Table 3).

Rotational velocity variations along the equatorial axis between 50 and 5000 au in L1448-NB cannot be reproduced satisfactorily by any single power-law model (χ² > 2, see Table 3). However, considering only the points at r < 400 au for L1448-NB, we obtain a power-law index of −0.9±0.2 with a good reduced χ² of 0.4, as expected for a collapsing and rotating envelope (see Table 3).

We found a positive index α for IRAS4A of 0.8 (see Table 3). It could be an indication of solid-body rotation of $\Omega =\frac{v_{\text{rot}}}{r}~$constant. However, we observe that the velocity in the PV_rot diagram decreases from 2000 to ~600 au and re-increases at small scales (see panel f of Fig. 16). Thus, the velocity gradient is not uniform on the scales traced by the PV_rot diagram as would be expected for a solid-body rotation (v_rot ∝ r). Moreover, points at radii < 600 au are consistent with an infalling and rotating envelope (see panel f of Fig. 16): considering only these points, we obtain a power-law index of −1.3 ± 0.6 with a good reduced χ² value of 0.6 (see Table 3). There is a dip in the C¹⁸O emission at r < 350 au that could be due to the opacity (see Figs. H.20 and H.19), thus, below this radius the information on velocities could be altered. To date, no observations have identified any solid-body rotating protostellar envelope. Numerical models also favor differential rotation of the envelope (Basu 1998). The interpretation of the velocity field as tracing solid-body rotation in the envelope of IRAS4A is therefore unlikely to be correct.

For the sources IRAS2A, IRAM04191, and L1157, the reduced χ² is also good (~1) when the PV_rot diagrams of these sources are ajusted by a model with a fixed index of α = −1 (see Table 3). We determine position and velocity from four different and independent methods and we did not consider the uncertainties of the connection between the different tracers and datasets. The uncertainty on the indices reported in Table 3 may thus be underestimated. On the other hand, although we determined the systemic velocity by maximizing the overlap of the blue and red points, this method does not allow a more accurate determination than 0.05 km s⁻¹. The systematic error of 0.05 km s⁻¹, added to previous velocity errors of the points in the PV_rot diagrams to take this uncertainty on the systemic velocity into account (see Appendix E), can be overestimated and thus lead us to underestimate the χ². For these three sources, the CALYPSO dataset only allow us to estimate a range of power-law indices between −1 and the α value reported in the fifth column of Table 3. The uncertainties on the indices reported in Table E are statistical errors and a systematic uncertainty of ±0.1 has to be added to account for the uncertainties in the outflow directions and thus the equatorial axis directions (see Table 1). Moreover, despite the choice of the equatorial axis, the rotational velocities could be contaminated by infall at the small scales along this axis due to the envelope geometry.

To conclude, the organized motions reported in the PV_rot diagrams and modeled by a power-law function with an index α ranging from −2 to 0 are consistent with differential rotational motions (Ω ∝ r^ɛ, with −3 < ɛ < −1 here). We identified rotational motions in all protostellar envelopes in our sample except in IRAS4A.

5.2 Distribution of specific angular momentum in the CALYPSO Class 0 envelopes

5.2.1 Specific angular momentum due to rotation motions

Assuming the motions detected along the equatorial axis are dominated by differential rotation for 11 of the 12 sources in our sample, we use the measurements reported in the PV_rot diagrams toderive the radial distribution of specific angular momentum in the protostellar envelopes due to rotation. In this part of the study, IRAS4A is excluded. The specific angular momentum is $j=\frac{J}{M}=\frac{I\Omega }{M}$ with the moment of inertia I defined as I ∝ Mr² (Belloche 2013). Thus, the specific angular momentum is calculated from the rotational velocity: j = v_rot(r) × r. We plot all the specific angular momentum profiles obtained for the CALYPSO subsample in panel b of Fig. 17. The individual distribution of specific angular momentum j(r) for each source is given in Appendix H. This is the first time that the specific angular momentum distribution as a function of radius within a protostellar envelope is determined homogeneously for a large sample of 11 Class 0 protostars. We performed a least-square fit of the j(r) profiles for each source individually, using a model of a simple power-law and a broken power-law model to identify the change of regimes. The broken power-law model function is defined with a break radius r_break as follows: $$\begin{array}{l}j(r)\propto {\left(\frac{r}{r_{\text{break}}}\right)}^{{\beta }_1} \text{when} r<r_{\text{break}},\\ j(r)\propto {\left(\frac{r}{r_{\text{break}}}\right)}^{{\beta }_2} \text{when} r>r_{\text{break}}.\end{array}$$

We report in Table F.1 the power-law indices fitting the best individual profiles and the associated reduced χ². For the broken power-law fits, only results with a reduced χ² better than the one obtained with a simple power-law model and with a break radius value r_break to which the j(r) profile is really sensitive, have been retained.

Two sources (L1448-NB and L1448-C) are better reproduced by a broken power-law model than a simple power-law model where the χ² are ~2: this allows us to identify a change of slope from a relatively flat profile to an increasing profile at larger radius in the envelope (β ~1), with break radii between 500 and 700 au. For the other sources, we also identified at scales of r < 1300 au a flat profile of specific angular momentum with β < 0.5 (L1448-2A, IRAS2A, SVS13-B, IRAS4B, L1527, and GF9-2) while the specific angular momentum profile at scales of r > 1000 au shows a steeper slope with β ~1 (L1521-F). However, two sources of the sample (IRAM04191 and L1157) stand out as the sources showing a steep increase in their specific angular momentum profile at scales of r < 1000 au (β ≥ 0.7), similar to the indices found at large radii in the sources showing a break in their j(r) profiles. We note that for the flat profiles (β < 0.5; L1448-2A, SVS13-B, IRAS4B, and GF9-2), and IRAM04191 and L1157, the specific angular momentum distribution is only constrained at scales < 1300 au. Most of the sources in our sample are better reproduced by a broken power-law model with a break radius (1000 ± 500) au and an increasing profile at larger radius in the envelope (β ~1.4) than a simple power-law model.

In his review, Belloche (2013) plotted the observed specific angular momentum as a function of rotation radius for several objects along the star-forming sequence. In this plot (panel a of Fig. 17), he identifies three regimes in the distribution of specific angular momentum, that can be broadly associated with different evolutionary stages. During the prestellar regime, on large scales, the apparent angular momentum of molecular clouds (Goldsmith & Arquilla 1985) and dense cores (Goodman et al. 1993; Caselli et al. 2002) appears to follow the power-law relation j ∝ r^1.6. During protostellar regime, between 100 and ~6000 au (0.03 pc), a few points in different protostellar envelopes suggest the specific angular momentum is relatively constant (j ~ 10⁻³ km s⁻¹ pc, Ohashi et al. 1997a; Belloche et al. 2002; Chen et al. 2007). During disk and binary regime, below 100 au, measurements in disks and Class II binaries (Chen et al. 2007) show a decrease of j following a trend characteristic of Keplerian rotation (j ∝ r^0.5).

Thus, from previous observational studies on rotational motions, finding a break at r ~ 1000 au between two trends of specific angular momentum within Class 0 protostars was unexpected. Although the velocity gradients observed in the outer part of the protostellar envelopes (r > 1000 au) are consistent with rotational motions, the observed j regime at these scales is not expected from pure rotational motions.

5.2.2 Apparent specific angular momentum

The radius range of j(r) distribution due to rotation motions is not homogeneous between sources (see Table F.1). To identify whether the radius of ~1000 au is a critical radius between two trends of specific angular momentum in each source, we derive the radial distribution of the apparent specific angular momentum |j_app| at all scales in the envelopes. To build |j_app|(r) distribution,we consider the gradients observed at all envelope scales, including also the reversed gradients and the shifted ones at scales of r ≳1000 au (see Fig. 19 and Sect. 5.4) which were excluded in the construction of the PV_rot diagrams in Fig. 16 because they are not consistent with rotational motions. By considering these velocity gradients, we add points in the outer envelopes but the trend observed in the inner envelopes do not change (see Tables F.1 and F.2). Thus, the |j_app|(r) distribution helps us to understand the origin of the trend and the velocity gradients observed at r > 1000 au. We plot all the apparent specific angular momentum profiles obtained for the CALYPSO subsample in Fig. 18. We also report the apparent specific angular momentum of IRAS4A which was identified as the only source that did not show any rotational motions in our sample (see Sect. 5.1). As for j profiles, we performed a least-square fit of the |j_app|(r) profiles for each source individually and we report the indices of the power-law models in Table F.2.

We create the median |j_app|(r) profile of the CALYPSO subsample. We first resampled the individual profile of each source in steps of 100 au and normalized it by the value at 600 au, then we took the median value of individual profiles at each radius step. The median profile is shown in gray on Fig. 18. From a broken power-law fit, we obtain a relatively flat profile (j_app ∝ r^0.3±0.3) at radii smaller than 1570±300 au and an increasing profile (j_app ∝ r^1.6±0.2) in the outer envelope. The radius of ~1600 au therefore appears to be a critical radius which delimits two regimes of angular momentum in protostellar envelopes: the specific angular momentum decreases down to ~1600 au and then tends to become constant.

The change of behavior of j_app above the break radius could be due to a change of tracer to study the kinematics in the outer envelope. However, we do not find any systematic consistency between r_app,break and the transition radius R_trans between the two tracers C¹⁸O and N₂H⁺. Even if for SVS13-B, R_trans is in the error bars of r_app,break, for three sources (L1448-NB, L1448-C, and IRAS4A) it is not consistent, and for IRAS4B, we do not observe a change of regime for j_app at r ~1600 au (see Tables F.1 and F.2). Moreover, for L1521F, only the C¹⁸O emission shows a velocity gradient allowing us to constrain the kinematics at scales of r > 1600 au (see Fig. 11) and we find the same trend of j_app (β_app ~1.2) than in all other sources where we used N₂H⁺ to constrain the outer part of the envelopes. The other sources (L1448-2A, IRAS2A, IRAM04191, L1527, L1157, and GF9-2) show a negative value of the apparent angular momentum at outer envelope scales due to a renversal of the velocity gradients (see Fig. 18, Table F.2, and Sect. 5.4.1). For two of these sources (L1448-2A and GF9-2) the radius where the gradient reverses along the equatorial axis, resulting in a negative j_app with respect to the inner envelope scales, is consistent with R_trans and r_app,break. For two sources (IRAS2A and IRAM04191), R_trans is consistent with the radius where the gradient reverses along the equatorial axis but not with r_app,break. For the last two sources (L1527 and L1157), the three radii are all different from each other. The different individual behaviors in the CALYPSO sample allow us to conclude that our finding that protostellar envelopes are characterized by two regimes of angular momentum does not result from our use of two different tracers.

From the median |j_app|(r) profile without normalization of the individual profiles at 600 au, we find a mean value of specific angular momentum in the inner parts of the envelopes (r < 1600 au) of ~6 ×10⁻⁴ km s⁻¹ pc. This value isslightly lower but compatible with the estimates made by Ohashi et al. (1997a) and Chen et al. (2007) in four Class 0 or I sources (j ~10⁻³ km s⁻¹ pc at r< 5000 au). It is also consistent with the studies by Yen et al. (2015a,b) which find values between 5 × 10⁻³ and 5 × 10⁻⁵ km s⁻¹ pc in the inner envelope (r < 1500 au). Yen et al. (2015a) estimate a specific angular momentum of ~ 5 × 10⁻⁴ km s⁻¹ pc at r~100 au for L1448-Cand L1527. Moreover, our values for L1157 are consistent with their upper limit estimate of 5 × 10⁻⁵ km s⁻¹ pc in the inner envelope (r < 100 au) of L1157.

The high angular resolution and the high dynamic range of the CALYPSO observations allow us to identify the first two regimes within individual protostellar envelopes: values at radii ≳1600 au (j_app ∝ r^1.6 on average, see Table F.2) seem to correspond to the trend found in dense cores at scales >6000 au while the values stabilize around ~6 ×10⁻⁴ km s⁻¹ pc on average at radii <1600 au. This study resolves for the first time the break radius between these two regimes deeper within the protostellar envelopes at around ~1600 au instead of ~6000 au. In a study of ammonia emission in the outer envelopes of two Class 0 objects, Pineda et al. (2019) find an increasing angular momentum profile scaling as r^1.8 from 1000 to 10 000 au, with values ~3 × 10⁻⁴ km s⁻¹ pc at radii ~1000 au. They do not detect the break around ~1600 au found in the CALYPSO sample. This break radius from which the profiles are found to be flat in the inner envelope may depend on the evolutionary stage of the accretion process during the Class 0 phase as suggested by Yen et al. (2015b). It could be due to the propagation of the inside-out expansion wave during the collapse (Shu 1977): assuming a median lifetime or half life of ~5 ×10⁴ yr for Class 0 protostar envelopes (Maury et al. 2011; see also Evans et al. 2009) at sound velocity (~ 0.2 km s⁻¹), one obtains a radius ~2000 au. This radius is on the same order of magnitude as the observed break radius between the two regimes observed in the distribution of specific angular momentum of sources in our sample. In this case, the break radius could be an indication of the age of the protostars, except for four sources (L1448-NB, IRAM04191, L1521F, and L1157) in our sample where we do not observe this break radius. Beyond this radius, the outer envelope may not have collapsed yet, and could therefore retain the initial conditions in angular momentum of the progenitor prestellar core.

This could be an explanation for the increase in angular momentum observed at the scales of r > 1600 au (j_app ∝ r^1.6 on average, see Table F.2), consistent with the prestellar stage (j ∝ r^1.6). We discuss the properties and physical origin of these two regimes in more details in the next sections.

5.3 Conservation of angular momentum in Class 0 inner envelopes

In this section, we focus on the relatively constant values of specific angular momentum observed in the inner envelopes at scales of r ≤ 1600 au in the j(r) profiles due to rotation motions (see Fig. 17). From these flat profiles, we find that the matter directly involved in the formation of the stellar embryo has a specific angular momentum ~3 orders of magnitude higher than the onein T-Tauri stars (j ~ 2 × 10⁻⁷ km s⁻¹ pc, Bouvier et al. 1993). We discuss constant values of specific angular momentum as conservation of angular momentum to test disk formation as a possible solution to the angular momentum problem.

It is difficult to constrain the time evolution of specific angular momentum for a given particle from angular momentum distributions which are snapshots of the angular momentum distribution of all particles at a given time during the collapse phase. During the collapse of a core initially in either solid-body rotation or differential rotation, particles conserve their specific angular momentum during the accretion on the stellar embryo (Cassen & Moosman 1981; Terebey et al. 1984; Goodwin et al. 2004). In the case of a protostellar envelope with a density profile ρ ∝ r⁻², an observed flat profile j(r) = constant requires, since each particle at different radii has the same specific angular momentum, an initially uniform distribution of angular momentum. This does not agree with the steep increase in specific angular momentum we observe at scales of r > 1600 au in the j(r) profiles. The break in the specific angular momentum profile could be due either to a faster collapse of the inner envelope caused by an initial inner density plateau (Takahashi et al. 2016) or to a change of dominant mechanisms responsible for the observed velocity gradients from inner to outer scales of the envelope.

In our sample, we distinguish eight sources with a relatively flat j(r) profile in the inner envelope (β < 0.5, see Table F.1): L1448-2A, L1448-NB, L1448-C, IRAS2A, SVS13-B, IRAS4B, L1527, and GF9-2. We estimate a centrifugal radius that would be obtained when the mass currently observed at ~100 au collapses and based on the mean value of specific angular momentum observed today ⟨j_100 au⟩ as follows: $$R_{\text{cent}}=\frac{{\langle j_{\text{100} \text{au}}\rangle }^2}{G M_{\text{100} \text{au}}}.$$(2)

The lower limit of the mass enclosed within 100 au, M_100 au, is the mass of the envelope $M_{\text{100} \text{au}}^{\text{dust}}$ estimated from the PdBI 1.3 mm dust continuum flux (Maury et al. 2019), assuming optically thin emission, a dust temperature at 100 au computed with Eq. (D.2) and corrected by the assumed distance (see Table 1). This mass estimate does not include the mass of the central stellar object, M_⋆ : since the embryo mass is unknown for most sources in our sample, we consider an upper limit of $M_{\text{100} \text{au}}=M_{\star }+M_{\text{100} \text{au}}^{\text{dust}}$ assuming M_⋆ = 0.2 M_⊙ for each source in our sample. This value of 0.2 M_⊙ corresponds to the stellar mass in the Class 0/I protostar L1527 from kinematic models of the Keplerian pattern in the disk (Tobin et al. 2012; Ohashi et al. 2014; Aso et al. 2017). The range of values for M_100 au are reported in the third column in Table 4. The calculated range of centrifugal radii associated with M_100 au is listed for each source in the fourth column in Table 4. We note that if M_⋆ of a source is smaller than that of L1527, then the centrifugal radius value we calculated is underestimated.

Since the embryo mass is uncertain and M_100 au may be underestimated if the dust emission is not optically thin, we compute the mass enclosed within r < 100 au, including the stellar embryo mass, needed to form a disk the size of $R_{\text{disk}}^{\text{dust}}$ with the ⟨j_100 au⟩ observed. The values are reported in the last column of Table 4.

For all the sources in our sample, the upper limits of the R_cent range are larger than 150 au and systematically larger than the continuum disk candidate radii $R_{\text{disk}}^{\text{dust}}$ from Maury et al. (2019) reported in the fifth column of Table 4. Moreover, Maret et al. (2020) only detect possible Keplerian rotation in two protostars in our sample (L1527 at radii ~90 au and L1448-C at r ~ 200 au) from the CALYPSO data. Thus, most R_cent values is expected to be less than 100 au. The upper R_cent values are probably overestimated because the contribution of the embryo mass to M_100 au is excluded.

Comparing the lower limits of the R_cent range with the candidate disk radius, we find a good agreement for most sources in our sample except for L1448-NB. We find a larger centrifugal radius (~500 au) than the observed disk size (<50 au) calculated considering only the main protostar L1448-NB1 of the binary system. Since in this study, we are interested in the kinematics of the whole system, we must consider all the continuum structure and not only that of the main protostar. Considering NB1 and NB2, Maury et al. (2019) resolve a circumbinary structure with a radius of (320 ± 90) au centered on the middle of the two components. Given the uncertainties, the latter value is consistent with the lower centrifugal radius estimated here. At these scales, Tobin et al. (2016a) observe a spiral structure surrounding the multiple system and interpreted it as a gravitationally unstable circumbinary disk. On the other hand, Maury et al. (2019) suggest that this component is due to orbital motions and tidal arms between the companions and Maret et al. (2020) do not detect any Keplerian rotation at radii < 170 au. Thus, the nature of this additional structure surrounding the multiple system is still unclear. As a consequence, the increase in specific angular momentum we measured at small scales could not only trace the rotation of the disk or the envelope but may be contaminated by gravitational instabilities due to orbital motions or a fragmented disk surrounding the system. Given the large uncertainties on the dust disk radii, we found a good agreement between centrifugal radii and $R_{\text{disk}}^{\text{dust}}$ for L1448-C and L1527. Moreover, the dust radius (50 au in L1527, Maury et al. 2019) does not necessarily exactly correspond to the centrifugal radius which was first detected in L1527 from observations of SO emission at 100 ± 20 au (Sakai et al. 2014). For this source, our estimate of R_cent (~70 au) is consistent with previous kinematic studies which detect a proto-planetary disk candidate with a radius of 50−90 au (Ohashi et al. 2014; Aso et al. 2017; Maret et al. 2020). Moreover, we observe a slight increase in the specific angular momentum we measured at r < 80 au. It could be due to the transition from the envelope to the disk.

The hypothesis of collapsing material with conservation of angular momentum, resulting in disk formation, at r < 100 au is therefore plausible for most sources in our sample. We notice that L1448-NB, in which Tobin et al. (2016a) claim the detection of a large candidate disk, shows the highest value of specific angular momentum at r < 1600 au of the CALYPSO sample, consistent with the angular momentum observed in the proto-planetary disks surrounding the T-Tauri stars which are estimated to be 1−6×10⁻³ km s⁻¹ pc (Simon et al. 2000; Kurtovic et al. 2018; Pérez et al. 2018). It could suggest an increase in the angular momentum of the disk during its evolution. In this case, the mean value of j(r) in the inner envelope would be lower in the less evolved than in the more evolved Class 0 protostars, and it would increase with time until reaching the value contained in the T-Tauri disks. In this scenario, L1448-NB would be one of the most evolved objects in the sample. However, the borderline Class 0/I protostar L1527, which is the most evolved object of the CALYPSO sample, has a specific angular momentum of ~6 × 10⁻⁴ km s⁻¹ pc at the inner envelope scales (see Table 4). In the same way, L1448-C has a specific angular momentum less than one order of magnitude lower than the values observed in the Class II disks while Maret et al. (2020) suggest the presence of a Keplerian disk in the inner envelope. As most of the CALYPSO inner protostellar envelopes have an order of magnitude less specific angular momentum than in Class II disks, we discuss below several possible explanations:

• (i)

  a part of the angular momentum inherited by the T-Tauri disks may not come from the rotating matter contained in the inner envelope accreted during the Class 0 phase. During the Class I phase, the mass accreted could come from regions further away from the envelope (r ≫1600 au) with a possibly higher specific angular momentum;

• (ii)

  disks may expand with time due to the transfer of angular momentum from their inner regions to their outer ones. Unfortunately, the specific angular momentum does not contain information about the mass. Large values of j(r) may be carried by low masses at the outer disk radius but may remain difficult to quantify. To this day, the mechanisms at work in disk evolution remain an open question. Some studies, for example, suggest that viscous friction may be responsible for the disk expansion (Najita & Bergin 2018);

• (iii)

  the specific angular momentum of the proto-planetary disks may be biased toward high values from historical, large and massive disks. A new population of small T-Tauri disks with radii between 10 and 30 au has been observed thanks to PdBI and ALMA (Piétu et al. 2014; Cieza et al. 2019). Assuming a small rotationally-supported disks around a stellar object including a total mass of 0.1−1 M_⊙ (Piétu et al. 2014), one expects a specific angular momentum between 10⁻⁵ and 10⁻⁴ km s⁻¹ pc, values which are similar to those we obtained in the inner Class 0 protostellar envelopes with the CALYPSO sample. However, to this day, no resolved observations of gas kinemactics of these small Class II disks allow us to estimate observationally their specific angular momentum.

5.4 Origin of the velocity gradients at r > 1600 au

At outer envelope scales, we detect velocity gradients (~ 2 km s⁻¹ pc⁻¹ at ~10 000 au, see Table 2) in the CALYPSO single-dish maps. They may not be directly related to rotational motions of the envelopes but rather to other mechanisms. Indeed, we observe in the CALYPSO dataset a systematic evolution of the orientation of the gradients between the inner and outer scales in the envelope (see Table 2). Figure 19 shows the orientation of the mean velocity gradient observed at different scales of the envelope with respect to the position angle of the gradient observed at scales ~100 au. The clear dispersion (~100° on average, see Fig. 19) of gradient position angle across scales within individual objects may be due to a change of dominant mechanisms responsible for the observed gradients from inner to outer scales of the envelope. From the literature, velocity gradients are often measured in the outer protostellar envelopes along the equatorial axis and they are interpreted as due to rotational motions or infall from a filamentary structure at scales of 1500−10 000 au (Ohashi et al. 1997a; Belloche et al. 2002; Tobin et al. 2011). In this section, we explore the possible origins of the velocity gradients found at scales of r > 1600 au and used to build the |j_app|(r) profiles (see Fig. 18).

5.4.1 Questioning the interpretation of counter-rotation

Six sources in the sample show a clear reversal of the orientation of the mean velocity gradient (|Θ − Θ_small| > 130°) from the inner to the outer envelope scales: IRAS4A, IRAS4B, L1527, IRAM04191, L1157, and GF9-2. We note that the kinematics at scales where we observed reversed velocity gradients (r > 1600 au) with respect to the small scales were not taken into account to build the PV_rot diagrams in Fig. 16, or the specific angular momentum profiles shown for the full sample in panel b of Fig. 17. Indeed, these profiles were aimed at characterizing the rotational motions in the envelopes and the angular momentum due to this rotation: a reversal of the rotation, if real, would require a more complex model than the power-law (v ∝ r^α) model we adopted in Sects. 4.3–5.3. In this section, we discuss these complex patterns in more detail.

In IRAM04191, we observed velocity gradients at outer envelope scales of r > 1600 au consistent with those observed previously by Belloche et al. (2002) and Lee et al. (2005) (Θ ~ 100°, see Table 2). However, in the inner envelope, we noticed a velocity gradient with a direction of Θ = −83° (see bottom middle panel in Fig. 10 and Table 2). In L1527, we found small-scale velocity gradients (Θ ~ 0° at r ~ 1000 au) consistent with those previously observed by Tobin et al. (2011) which are in the opposite direction compared to the large-scale one (r ~ 8000 au, Goodman et al. 1993). Tobin et al. (2011) interpret this reversal of velocity gradients as counter-rotation but it also could be due to infalling motions that dominated the velocity field at the outer envelope scales (Harsono et al. 2014).

Our study suggests that reversals of velocity gradients are common in Class 0 protostellar envelopes. However, the asymmetrical velocity gradients (for IRAS4B, GF9-2), the filamentary structures traced by the integrated intensity at scales of r > 2000 au (for IRAS4A, IRAS4B, L1527, and GF9-2), and a strong external compression of the cloud hosting IRAS4A and IRAS4B (Belloche et al. 2006) lead us to exclude the observed reversed gradients as counter-rotation of the envelope. Moreover, only MHD models with Hall effect succeed to form envelopes in counter-rotation. These models form a thin layer of counter-rotating envelopes at the outer radius of the disk (r ~ 50–200 au; Tsukamoto et al. 2017). This envelope layer is in counter-rotation compared to the formed disk and the protostellar envelope at r > 200 au as a consequence of the Hall effect generated by the rotation of the disk which changes the angular momentum of the gas at the disk outer radius. Therefore, these models cannot explain the inversions of rotation in the different layers of the envelope at scales of r > 3000 au as observed in our sample. Historically, the gradients observed from single-dish mapping at r > 3000 au have been used to quantify the amplitude for the angular momentum problem. However, incorrectly interpreted as pure rotational motions in the envelope, the resulting angular momentum measurements and the expected disk radii would be significantly overestimated.

Recent studies on the angular momentum of the protostellar cores from hydrodynamical simulations of star formation are questioning the standard model of star formation from a collapsing core initially in solid-body rotation (Kuznetsova et al. 2019; Verliat et al. 2020). They show that the angular momentum of synthetic protostellar cores is not directly related to the initial rotation of the synthetic cloud, and Keplerian disks can be formed from a simple non-uniform perturbation in the initial density distribution. In this scenario, the angular momentum observed in inner protostellar envelopes and disks may not been inherited from larger-scale initial conditions but generated during the collapse itself.

5.4.2 Contribution of infalling motions and core-forming motions

The misalignments between the gradients observed in the envelopes at inner and outer envelope scales suggest a change of dominant mechanisms at r > 1600 au. At large scales, infalling motions of the envelope can dominate rotational motions. In the hypothesis of a flattened infalling envelope, infall motions are expected to produce a velocity gradient projected in the plane of the sky that is oriented along the minor axis of the envelope, namely at the same position angle as the outflow. In L1448-NB, SVS13-B and L1527, we detect velocity gradients aligned with the outflow axis at r > 3000 au while at small scales the gradients are consistent with the equatorial axis (see Table 2). These three sources could be good candidates of the transition from collapse to rotation between large and small scales. This scenario is also suggested in the study of Ohashi et al. (1997b). They suggested that at outer envelope scales of r ~ 2000 au, the protostellar envelope L1527 is not rotationally supported (v_rot ~ 0.05 km s⁻¹) but is dominated by the collapse (v_inf ~ 0.3 km s⁻¹).

Currently, there are very few constraints on the infall velocities at scales of r > 1600 au in the CALYPSO protostellar envelopes. Belloche et al. (2002) estimate an infall velocity of v_inf ~ 0.15 km s⁻¹ at r ~ 1000 au from radiative transfer modeling of CS and C³⁴S emission in IRAM04191. In the dense core L1544, Tafalla et al. (1998) suggest also an infall velocity of ~0.1 km s⁻¹ at scales >3000 au. The velocity offset, with respect to the systemic velocity assumed for each source, found along the equatorial axis at > 1000 au with CALYPSO is reported in Table 5. For most sources, we find typical velocity offsets ≲0.3 km s⁻¹ at scales of 1600 au (see Table 5), consistent with infall velocities found in IRAM04191 and L1544, except for IRAS4A. IRAS4A harbors a velocity of ~ 0.5 km s⁻¹ at r ~ 1000 au. This result is consistent with those of Belloche et al. (2006) at ~2000 au. They suggest that a fast collapse is triggered by an external compression from the cloud in which the source is embedded. Thus, for all sources in our sample, the velocity gradient misalignment could be due to a change of mechanism dominating the velocities projected on the line of sight. This suggests either rotational velocities much smaller than infall velocities or a non axisymmetric geometry of the kinematics at outer envelope scales.

Moreover, Herschel observations have shown that most solar-type prestellar cores and protostars form in filaments (André et al. 2014). Indeed, the column density maps of the Herschel Gould Belt Survey program⁴ (André et al. 2010)reveal that the CALYPSO protostars are embedded in or lie in the immediate vicinity of filamentary structures with $N_{{\text{H}}_2}>{10}^{21}$ cm⁻². Thus, the large-scale kinematics in protostellar envelopes could be contaminated or dominated by the kinematics of the filaments. Kirk et al. (2013) studied the velocity field of Serpens-South in the Aquila molecular cloud and showed a complex kinematics with longitudinal collapse along the main filament, radial contractions, and accretion streams from the cloud tothe main filament. The longitudinal collapse of the filament could be responsible for the large-scale gradients observed in our protostellar envelopes as observed in the Serpens-Main region by Dhabal et al. (2018). Several studies also highlighted transverse velocity gradients perpendicular to the main filament that suggested the material may be accreting along perpendicular striations (Palmeirim et al. 2013; Dhabal et al. 2018; Arzoumanian et al. 2018; Shimajiri et al. 2019). Palmeirim et al. (2013) estimate the velocity of the infalling material to be 0.5−1 km s⁻¹ at r ~ 0.4 pc in the B211/L1495 region. In our velocity maps at 10 000 au along the equatorial axis, we measure typical velocities < 0.3 km s⁻¹ in most of the sources (see Table 5) except in L1448-2A, IRAS4A, and IRAS4B. These three sources exhibit velocities of 0.5−1 km s⁻¹ consistent with infall velocities estimated at filamentary scales. In these cases, host-filament motions could dominate the kinematics in outer protostellar envelopes (r > 1600 au).

5.4.3 Contamination by turbulent motions from cloud scales

All sources of the CALYPSO subsample (except L1448-NB) reveal a steep increase in apparent specific angular momentum with the radius at ~1600−10000 au scales, with an average trend of j_app ∝ r^1.6±0.2 (see Table F.2 and Fig. 18). This trend is similar to that observed in prestellar cores and clumps at scales > 10 000 au (see Figs. 17 and 18). Indeed, Goodman et al. (1993), Caselli et al. (2002) and Tatematsu et al. (2016) show a trend between the size of prestellar cores and their observed mean angular momenta at scales on the order of 0.1 pc: j(r) distribution scaling as r^1.2−1.7. From this dependency of j with core radius and the linewidth-size relation, Tatematsu et al. (2016) suggest that non-thermal motions (turbulence) are related to the origin of angular momentum observed in these 0.1 pc cores.

Burkert & Bodenheimer (2000) studied numerical models of turbulent molecular clouds with a symmetric density profile and a Gaussian or random velocity field. They showed that 0.1 pc cores with random motions exhibit most of the time velocity gradients that, interpreted as rotation, would have specific angular momentum values of j ~ 3 × 10⁻³ km s⁻¹ pc (10²¹ cm² s⁻¹) and would scale as j ∝ r^1.5. This is in good agreement with observed values at 0.1 pc from the literature (Goodman et al. 1993; Caselli et al. 2002; Tatematsu et al. 2016). Our observations showing a trend of j ∝ r^1.6 at scales ~1600–5000 au could be either a signature of the turbulent cascade from the large-scale ISM propagating with subsonic properties to 1600 au envelope scales, or gravitationally-driven turbulence due to large-scale collapse motions at the interface between filaments and cores (Kirk et al. 2013).

From the analysis of the gas velocity dispersion in molecular line observations, Goodman et al. (1998) and Pineda et al. (2010) identify the dense cores at a typical scale of 0.1 pc as the first velocity-coherent structures decoupled from the turbulent cloud. In this case, we would expect a quiescent structure with subsonic motions at radii < 0.1 pc and the interpretation of ISM turbulent cascade with supersonic motions as a consequence of the steep increase of j_app at scales < 10 000 au would no longer be valid. Except L1448-2A, IRAS4A, and IRAS4B which exhibit velocities 0.5−1 km s⁻¹ consistent with supersonic turbulent motions, all sources in our sample show typical velocities ≲0.3 km s⁻¹ consistent with subsonic-transonic turbulent motions. This could suggest that the power-law behavior of j_app ∝ r^1.6 observed in the outer envelopes (r > 1600 au) could be a scaling law due to the tail of a low velocity subsonic-transonic turbulent cascade.

At scales of r > 1600 au, we observe typical velocity linewidths ≲1 km s⁻¹ (see Fig. 20). We note that the linewidths tend to decrease from ~1600 au to larger scales in the outer envelopes and they do not show scaling laws with the radius as expected from turbulent motions in the ISM (Larson 1981; see Fig. 20 and Table G.1), but the velocity structure at these scales seems to show multiple components in velocity for several sources (L1527, L1448-C, IRAS2A, SVS13-B, IRAS4A, IRAS4B). As we can not disentangle them and identify exactly which component comes from the outer envelope or the host cloud for example, we need either a more elaborate model than a Gaussian or a HFS model to analyze the spectra or a more suitable tracer to determine more robustly the linewidths of the outer envelopes.

6 Conclusions

In the framework of the CALYPSO survey, we analyzed the kinematics of Class 0 protostellar envelopes. The main results of our study are listed below.

• 1.

  We identify differential rotation motions in 11 sources in a sample of 12 Class 0 protostellar envelopes. The only exception is IRAS4A: the motions reported in the PV_rot and modeled by a power-law function are consistent with a solid-body rotation, but the velocity gradient is not uniform in the inner envelope at r < 2000 au as wouldbe expected.

• 2.

  This is the first time that the specific angular momentum distribution as a function of envelope radius is determined homogeneously for a large sample of 11 Class 0 protostars. The high angular resolution and the high dynamic range of the CALYPSO observations allow us to identify two distinct regimes: the apparent specific angular momentum decreases as j_app ∝ r^1.6±0.2 down to ~1600 au and then tends to become relatively constant around ~6 × 10⁻⁴ km s⁻¹ pc down to ~50 au.

• 3.

  The values of specific angular momentum measured in the inner Class 0 envelopes suggest that material directly involved in the star formation process (<1600 au) typically encloses the same order of magnitude in specific angular momentum as what is inferred in small T-Tauri disks (r ~ 10 au). The constant values of j at 50−1600 au allow us to determine good estimates of the centrifugal radius in the Class 0 protostars of the CALYPSO sample, which compare well with the disk radii estimated from the dust continuum (Maury et al. 2019). This suggests that the specific angular momentum is conserved during the accretion on the stellar embryo, resulting in disk formation.

• 4.

  At scales of r > 1600 au, we conclude that the velocity gradients observed in the outer envelope with respect to small scales are not due to pure rotational motions or counter-rotation motions but related to other mechanisms. Historically, the gradients observed from single-dish mapping at r > 3000 au have been interpreted as rotation and used to quantify the amplitude of the angular momentum problem for star formation. Thus, if the gradients are incorrectly interpreted as rotation, the angular momentum problem for star formation and the expected disk radii may have been significantly over-estimated. Moreover, we find no robust hints that envelopes are rotating with typical velocities higher than the sound speed at scales of r > 1600 au. This suggests that the origin of angular momentum in the outer protostellar envelopes could be the gravitationally-driven turbulence due to large-scale collapse motions at the interface between filaments and cores, or the dissipation of the large-scale ISM turbulent cascade propagating with subsonic velocities to 1600 au envelope scales.

Acknowledgements

We thank the IRAM staff for their support carrying out the CALYPSO observations. 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). M.G. thanks the Max-Planck Institute for Radio Astronomy for its support toward the end of this work. We would like to thank Cecilia Cecarelli for comments and suggestions on the estimation of column density, and Nagayoshi Ohashi and Jaime E. Pineda for valuable discussions on the interpretation. We also thank the referee for the useful comments which helps to improve this paper.

Appendix A Details of IRAM 30m observations

Appendix B Properties of molecular line emission maps

Appendix C Details of construction method of PV_rot diagrams

In this appendix, we present in detail our method of quantifying centroid velocity variations. We have built large dynamic range PV_rot diagrams (see Sect. 4.3) by determining the position and the velocity in the three datasets. We analyzed the PdBI datasets in the (u,v) plane while the combined and 30m datasets are analyzed in the image plane.

C.1 Analysis of PdBI datasets in the (u,v) plane

The PdBI datasets provide the best spatial resolution available in the CALYPSO sample, allowing us to constrainthe velocity variations at the smallest scales of the protostellar envelopes. Because the (u,v)-coverage of our PdBI dataset is limited, inside the C¹⁸O PdBI emitting size of the source (see R_i in Tables C.1 and B.3), we work directly on the visibilities to avoid the complex imaging and deconvolution processes involved in image plane analysis and to characterize the peak of emission at a given velocity with a higher astrometric precision. We fit the PdBI visibilities of each channel with an elliptical Gaussian source in order to determine the centroid position of the emission in each velocity bin. In the cases where the emission shows irregular intensity distributions with multiple peaks, this analysis recovers the most intense peak at a given velocity. Our results are consistent with those of Maret et al. (2020) who are doing a similar analysis in the (u,v) plane. We note that they are not strictly identical because Maret et al. (2020) consider only pixels at scales < 2′′ from the central position to study the disk kinematics whereas we are interested in the kinematics of the envelope scales between 50 and 5000 au. Maret et al. (2020) also did some tests from synthetic visibilities with an interferometer and find that working in the uv-plane provides a good estimate of the rotation profile at small scales.

We report in the PV_rot diagrams the fitting results of sources showing a velocity gradient in at least two channels in the PdBI channel maps (see Appendix H). We identify a gradient as a variation in the central position of the C¹⁸O emission from one side of the continuum peak in the redshifted velocities to the other side in the blueshifted velocities with respect to the systemic velocity. The fitting results in channel maps of sources which do not exhibit such an organized spatial distribution of velocities are considered as upper limits. We only consider the fitting results in channels where the C¹⁸O emission is detected with a S/N higher than 5 and with a fit central position consistent with the equatorial axis, namely at a position angle < |45°| with respect to the equatorial axis. These criteria on the position angle prevent the selection of PdBI C¹⁸O points from being affected by the typical error of ±10° on the direction of the equatorial axis. When these three criteria are satisfied, we project the fit centroid position onto the equatorial axis to constrain the PV_rot diagram at small scales. The position errors are derived from our modeled elliptical Gaussian emitting source. The errors on the velocity are related to the channel width (0.2 km s⁻¹ in C¹⁸O emission from the PdBI datasets). We note that most of the CALYPSO sources show an optically thin emission at small scales (except IRAS4A, see Appendix H.19), thus, the intensity peak may be located in the equatorial plane where the density is typically higher. In the case that the centroid position in a channel map projected onto the plane of the sky is not along the equatorial axis while the centroid position in 3D actually belongs to the equatorial plane, the apparent distance in the channel map would underestimate the true distance of the centroid emission to the rotation axis because of its location in the third dimension. We did not apply any further correction on the velocities due to the position angle with respect to the equatorial axis of the fit central position. The associated rotational velocity of a fit central position close to 45° with respect to the equatorial axis could be over-estimated because more contaminated by the ejection or collapse motions at small scales than a fit central position closer to the equatorial axis. Restraining the analysis to the fit central positions with a position angle < |20°| with respect to the equatorial axis results in an additional uncertainty of ±0.1 on the indices reported in Table E. This uncertainty is consistent with the systematic uncertainty of ±0.1 which has to be added to account for the uncertainties in the equatorial axis directions (see Sect. 5.1).

Seven sources in our sample show a small-scale velocity gradient aligned along the equatorial axis in C¹⁸O emission: L1448-2A, L1448-NB, L1448-C, IRAS2A, SVS13-B, L1527, and GF9-2 (see Figs. H.1, H.4, H.8, H.12, H.16, H.35, and H.43). For IRAM04191and L1521F (see Figs. H.27 and H.31, respectively), the weak detection in C¹⁸O emission from the PdBI datasets does not allow the diagram to be constrained to the smallest scales of the envelope (r < 60 au). For IRAS4A, IRAS4B, and L1157 (see Figs. H.21, H.23, and H.39, respectively), the channel maps do not exhibit an organized spatial distribution of velocities along the equatorial axis at r ≲ 350 au: the central emission fit show a position angle >|45°| with respectto the equatorial axis, suggesting a contamination by the outflow kinematics. However, PdBI velocity maps exhibit weak gradients along the equatorial axis but working with the visibilities does not allow us to disentangle it from the emission from the outflows. Thus, for these three sources, we analyzed the PdBI observations in the image plane as described in the next section to constrain the PV_rot diagrams at scales of r < 600 au.

C.2 Analysis of combined and 30m datasets in the image plane

We analyze the combined and 30m datasets in the image plane to probe the intermediate and outer scales of the envelope. For the pixel (i.e., a position) along the equatorial axis, we report the centroid velocity determined in the velocity maps in Sect. 4.2. The errors on the velocity given by the Gaussian or HFS functions are reported on the PV_rot diagrams (see Fig. 16 and Appendix H). The position errors are related to the pixel size, which is itself relatedto the spatial resolution of the datasets.

To build PV_rot diagrams weonly consider the velocity gradients with blue and red-shifted components with respect to the systemic velocity, centered on the central protostar and along the equatorial axis axis as expected from rotational motions. The values from velocity maps that do not exhibit such an organized spatial distribution of velocities are considered as upper limits on the rotational velocity. Another requirement is the robustness of the centroid velocity offset with respect to the systemic velocity. We only considered spectra with a S/N higher than 5. From a Gaussian measurement, we assume the centroid velocity cannot be robustly determined with an accuracy better than one third of the spectral resolution: we only consider the C¹⁸O points where the relative centroid velocity, v − v_lsr, is $\ge \frac{\Delta v}{3}$, with Δv the spectral resolution. From a HFSmeasurement, the centroid velocity is determined much more precisely thanks to the larger number of components: we only consider the N₂H⁺ points where $v-v_{\text{lsr}}\ge \frac{\Delta v}{3\times \sqrt{6}}$. We note that the weakest component which is 7 times weaker than the strongest one is never detected with a S/N higher than 2 in the spectra.

C.3 Construction method of PV_rot diagrams

The analysis described above allows us to determine the centroid position of the emission at a given velocity and the centroid velocity at a given position, respectively. By putting the results end-to-end, we build a PV_rot diagram with a high dynamic range from 50 up to 5000 au for each source as follows (see the example of L1448-C in Fig. 15):

Firstly, at the smallest scales resolved by our dataset (~0.5′′) and up to the PdBI C¹⁸O emission size radius of the sources R_i (see Tables B.3 and C.1, and Fig. 15), we use the PdBI C¹⁸O datasets to constrain the PV_rot diagram (see label “C¹⁸O PdBI” in Fig. 15).

Secondly, since the C¹⁸O extended emission is filtered out by the interferometer, we used the combined C¹⁸O emission to populate the PV_rot diagram at angular radii between R_i and the 30m half-power beam width (HPBW ~ 6′′ at the C¹⁸O frequency, see Table B.1). However, because the C¹⁸O molecule freezes onto dust ice mantles at temperatures below ~20K, the maximum radius up to which it remains the best tracer of the inner envelope can be smaller than the 30m HPBW. To determine the maximal radius R_trans where the C¹⁸O molecule remains the best tracer, we calculate the C¹⁸O and N₂H⁺ column densities along the equatorial axis from the combined integrated intensity maps (see Appendix D and green points in Fig. 15). R_trans is defined one of the following criteria, depending on the source (see Table C.1):

• (i)

  the radius from which the C¹⁸O column density changes from a smooth decreasing profile to a noisy dispersion,

• (ii)

  the maximum radius where the C¹⁸O emission is detected with a S/N higher than 5 along the equatorial axis,

• (iii)

  the radius below which the N₂H⁺ column density no longer traces the inner region and the profile flattens,

• (iv)

  the radius where the N₂H⁺ profile reaches its maximum in column density.

For most sources, R_trans is smaller than the 30m HPBW. We used the combined C¹⁸O emission map to constrain the gas kinematics in the PV_rot diagram up to R_trans (see label “C¹⁸O combined” in Fig. 15).

Then, at radii r > R_trans, the N₂H⁺ emission traces better the envelope dense gas. We use the combined N₂H⁺ emission maps to analyze the envelope kinematics at intermediate scales, from R_trans to the 30m HPBW (~14′′ at the N₂H⁺ frequency, see Table B.2). However, the N₂H⁺ column density profile can reach a minimum value at radii R_trans < r < 30m HPBW due to the sensitivity of the combined datasets. In this case, the combined dataset is no longer the better dataset to provide a robust information on the velocity. We defined R_int the maximum radius up to which we use N₂H⁺ emission from the combined dataset (see Table C.1 and see label “N₂H⁺ combined” in Fig. 15).

Finally, beyond R_int, we use the 30m N₂H⁺ emission map to populate the PV_rot diagram up to the largest scales of the envelope (see label “N₂H⁺ 30m″ in Fig. 15). We note R_out the maximum radius of the protostellar envelopes (see Table C.1).

The sources in our sample show specific individual behaviors, therefore we adapted the method of building the PV_rot diagram described above on a case-by-case basis. Nevertheless we respected the transition radii as closely as possible (see Table C.1).

Appendix D Details of column density estimate

To interpret the gas kinematics using both C¹⁸O and N₂H⁺ emission in a robust way, it is important to constrain from which scales in the envelope their emission comes from. This is especially important because we merged single-dish and interferometric datasets which are sensitive to different scales of the envelope (see Tables B.3 and B.4). Therefore, we built column density maps to quantify the scales probed by C¹⁸O and N₂H⁺ molecular lines. We assumed that the observed molecular transitions are optically thin (τ_ν << 1) and we used Eq. (D.1) defined in Mangum & Shirley (2015): $$\begin{array}{ll}{N}_{\text{thin}}=\hfill & \frac{3 k}{8 {\pi }^3 {\mu }^2 B_{\text{0}} R_{\text{i}}} \left(T_{\text{ex}}+\frac{hB}{3k}\right) \frac{\exp \left(\frac{E_{\text{u}}}{T_{\text{ex}}}\right)}{\exp \left(\frac{h\nu }{k{T}_{\text{ex}}}\right)-1}\hfill \\ \hfill & \times \frac{{\displaystyle \int T_{\text{R}}}\text{\hspace{0.17em}d}v}{J_{\nu }(T_{\text{ex}})-J_{\nu }(T_{\text{bg}})},\hfill \end{array}$$(D.1)

where E_u is the energy of the upper level, μ is the dipole moment of the molecule, B₀ the rotational constant, h Planck’s constant, k Boltzmann’s constant, R_i the relative intensity of the component if the transition has a hyperfine structure, T_ex the excitation temperature, and T_bg the cosmic background temperature (Mather et al. 1994). J_ν is the effective radiation temperature defined by $J_{\nu }(T)=\frac{h\nu }{k} \frac{1}{\exp (\frac{h\nu }{kT})-1}$. We report the values of each parameter for the two emission lines in Table D.1.

Assuming local thermodynamical equilibrium (LTE), we can determine the excitation temperature from the gas temperaturewithin the emission size of each tracer observed in the combined dataset (see Table C.1). This hypothesis is valid when the density is higher than the critical density of the transition of C¹⁸O (2−1) and N₂H⁺ (1₀₁−0₁₂), estimated at ~8.4 × 10³ cm⁻³ (Flower 2001) and ~2.6 × 10⁵ cm⁻³ (Daniel et al. 2005), respectively at T ~ 10 K. Belloche et al. (2002) estimate for IRAM04191, one of the lowest envelope masses of our sample (see Table 1), a density higher than 10⁵ cm⁻³ up to r ~ 4000 au. Thus, LTE is a good assumption for both transitions in our sample. As the dust temperature in the envelope is a good approximation for the gas kinetic temperature (Ceccarelli et al. 1996), we used dust temperature profiles from CALYPSO PdBI observations assuming that the temperature distribution depends on the radius r from the central stellar object and the internal luminosity of the source L_int (see Table 1) as follows (Terebey et al. 1993): $$T_{\text{dust}}(r)=60 \times {\left(\frac{r}{13400 \text{au}}\right)}^{-q} {\left(\frac{L_{\text{int}}}{{10}^5{L}_{\odot }}\right)}^{q/2},$$(D.2)

with q = 0.4 for the CALYPSO sample (Terebey et al. 1993; Maury et al. 2019). We weighted the radial temperature distribution with the dust column density as a function of r, that is the amount of material at each radius, to estimate a single robust mean value of gas temperature (see Table D.2). For C¹⁸O, we used the PdBI flux at 1.3 mm calculated at different radii by Maury et al. (2019). For N₂H⁺, we assumed a density profile of ρ ∝ r^−p, with the index p equal to the values determined at 3 mm by Maury et al. (2019). The excitation temperature values adopted for each source to derive the column density are given in Table D.2. There are two exceptions in the CALYPSO sample: IRAM04191 and L1521F, which have the lowest luminosities of the CALYPSO sample (see Table 1). The dust temperature profile of these sources from the PdBI dataset cannot be determined because the temperature profile is dominated by external heating. Belloche et al. (2002) and Tokuda et al. (2014, 2016) estimate the dust temperature at ~10 K at scales ~2000 au. From Eq. (D.2) and from values of L_int estimated for these two sources (0.05 L_⊙ for IRAM04191, André et al. 2000, and 0.035 L_⊙ for L1521F, Tokuda et al. 2016) we determined mean excitation temperatures of 20 K for the C¹⁸O molecule and 10 K for the N₂H⁺ molecule (see Table D.2).

We note that we did not use the excitation temperature determined by fitting the HFS line profile for the N₂H⁺ molecule. For most sources in the CALYPSO sample, the average spectrum from the combined dataset within the central 14′′ ×14′′ does not follow the relative intensities of the hyperfine components expected under LTE conditions (3, 3, 7, 5, 3, 5, 1 from the isolated hyperfine component 1₀₁−0₁₂; Endres et al. 2016). Thus, either the emission does not satisfy the LTE conditions or this is an effect of thepartial opacity of emission. Moreover, some adjacent hyperfine components are too close to be spectroscopically separated despite the good spectral resolution of the combined dataset (0.13 km s⁻¹ at 3 mm). To minimize the error propagation, we only considered the isolated hyperfine component (1₀₁−0₁₂) to determine the N₂H⁺ column density (Caselli et al. 1995, 2002).

We determined the C¹⁸O (2−1) and N₂H⁺ (1₀₁−0₁₂) opacities forour source sample assuming optically thin emission and using the maximum temperature of the emission spectrum T_peak: $${\tau }_{\nu }=\frac{T_{\text{peak}}}{J_{\nu }(T_{\text{ex}})-J_{\nu }(T_{\text{bg}})}.$$(D.3)

We calculated the opacity on the spectra of each pixel from the combined dataset for both molecular transition. We used the excitation temperature values mentioned above and a standard value of the cosmic background temperature T_bg = 2.7 K (Mather et al. 1994). We report the average τ_ν on the size of the combined maps of C¹⁸O and N₂H⁺ emissions (5′′ × 5′′ and 20′′ × 20′′, respectively) in Table D.2. We noticed that the sources have average τ_ν < 0.4 but some pixels have opacity values >0.5. We determined the column density from each pixel using Eq. (D.1) and the parameter values listed in Table D.1 for the C¹⁸O and N₂H⁺ molecules. For the calculation of column density, we distinguished two opacity regimes according to the value of τ_ν of each spectrum. We consider an optically thin regime when τ_ν < 0.4 and, in thiscase, the density column is determined by Eq. (D.1). When τ_ν ≥ 0.4, the optically thin hypothesis is no longer a reasonable hypothesis and a correction factor $\frac{{\tau }_{\nu }}{1-e^{-{\tau }_{\nu }}}$ (Goldsmith & Langer 1999) must be applied to the Eq. (D.1) to calculate the column density. Indeed, the column density determined in the optically thin hypothesis by Eq. (D.1) is underestimated by more than 20% for a value of τ_ν = 0.4.

Figure D.1 shows as an example the column density maps from the combined dataset obtained for L1448-2A. The column density maps of the other sources from our sample are provided in Appendix H (see Figs. H.5, H.9, H.13, H.17, H.19, H.24, H.28, H.32, H.36, H.40, and H.44).

Appendix E Systemicvelocity estimate

The systemic velocities of the sources are the reference velocities to quantify centroid velocity variations in the PV_rot diagrams. The systemic velocity is defined as the average velocity of an object. The linear velocity gradient modeling allows us to determine a systemic velocity, v₀, in the different datasets and tracers. For a given source, we noticed that v₀ is different depending on the region of the envelope considered (see Table 2).

In the first place, we determine these values by fitting the HFS line profile of the optically thin transition N₂H⁺ (1−0) from the 30m datasets. Table E.1 reports these first estimates and compares them to the values found in the literature. The differences of a few tenths of a km s⁻¹ between ourvalues and the literature can be explained by the different tracers used in our study compared to the literature. Indeed, different tracers do not trace exactly the same material in the protostellar envelopes, as we observed with the C¹⁸O and N₂H⁺ in our study(see Sect. 4.2 and Table 2). There could also be systematic uncertainty on the rest frequencies or on the calibration of the spectral axis for instance. This first estimate in the quieter parts of the envelopes, could be more robustthan a value determined in the inner parts which can be dominated by rotation, collapse, or ejection. On the other hand, the turbulence at large scales could affect the systemic velocity if the outer envelope is turbulent or if there is a contamination of the turbulence on the line of sight from cloud scales. In the case of an axisymmetric envelope with a symmetric kinematics (infall and rotation), the systemic velocity could be more accurate at small scales.

If the reference velocity value is well chosen, blue- and red-shifted velocity points should overlap if the envelope is axisymmetric. In practice, the gradients may not be symmetric with respect to the central position of the source. This asymmetry can cause shifts in position between the blue-shifted and the red-shifted emission. We built the PV_rot diagrams by maximizing the overlap between blue- and red-shifted velocity points to produce the best symmetric pattern. We applied this method independently for the PV_rot diagram points from the C¹⁸O and N₂H⁺ emission to determine the best couple (v_sys, r_orig) adapted to each tracer, with v_sys the systemic velocity and r_orig the central position of the gradient along the equatorial axis. We fit a power-law function (see Sect. 5.1) exploring a range of ±0.7 km s⁻¹ around the first estimate of the systemic velocity determined at outer envelope scales and a range of ±0.2′′ along the equatorial axis around the dust continuum peak (see Table 1). Table E.2 gives the reduced χ² associated with the best couple of parameters corresponding to the best overlap of the blue- and red-shifted points. This method does not allow for a more accurate determination of the systemic velocity value than 0.05 km s⁻¹ given the errors on the velocity of the points populating the PV_rot diagrams. We therefore added a systematic error of 0.05 km s⁻¹ to previous velocity errors determined above for both tracers. The central positions and systemic velocities obtained with this method are on average within < 0.2 km s⁻¹ and < 0.1′′ of the first estimates of these values using the average N₂H⁺ spectrum (see Appendix E and Table E.1) and the dust continuum peak (see Table 1).

Appendix F Details on j distributions

Appendix G Details on linewidth distributions

Appendix H Comments on individual sources

H.1 L1448-2A

L1448-2A (also known as L1448-IRS2 or Per-emb-22) is a Class 0 protostar in the L1448N complex in the Perseus cloud at a distance previously estimated to be 235 pc (Hirota et al. 2011) but determined at (293 ± 20) pc by recent Gaia parallax measurements (Ortiz-León et al. 2018). CO emission maps showed the presence of a bipolar flow in the northwest-southeast direction (O’Linger et al. 1999). Another study by Wolf-Chase et al. (2000) suggests that the source shows distinct bipolar outflows, namely a signature of the presence of a binary system. Observations with VLA confirmed that L1448-2A is a binary system separated by ~170 au (Tobin et al. 2016b). 1.3 mm dust continuum emission also resolves the binary system with a separation of ~180 au (Maury et al. 2019, see Table 1). Podio & CALYPSO (in prep.) report two angle values (+63° (blue), +140° (red), see Table 1) because the observed red and blue flow cavities are not aligned.

Figures 1 and 3 show the integrated intensity and centroid velocity maps obtained for this source from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission respectively. The C¹⁸O emission from the PdBI datasets does not peak on the main protostar L1448-2A1 but in the middle of the binary system (see Table 1 and Fig. 1). Thus, the origin of the coordinate offsets is chosen to be the middle of the binary system (RA: 03^h25^m22^s.380, Dec.: 30°45′13′′.21; see Table 1) to study the kinematics in the protostellar envelope of L1448-2A.

The gradients observed in the PdBI and combined velocity maps of the C¹⁸O emission (see bottom left and middle panels on Fig. 3) have an orientation of 99° ≤ Θ ≤107°, they show a difference >60° with respect to the equatorial axis (see Table 2). Therefore, these gradients at r < 700 au are a merging of the kinematics between several mechanisms: ejection by outflows, orbital motions of the binary system, and rotational motions of the envelope along the equatorial axis. The N₂H⁺ emission from the PdBI datasets (see top left panel on Fig. 3) shows reversed gradients at scales of r ~ 2000 au with respect to those observed at smaller scales (Θ ~ −87°, see Table 2). In the outer envelope (see top middle and right panels on Fig. 3), the velocity gradients are not continuous in the same direction and not source-centered. Therefore, these gradients seem to be due to an external contamination. We notice that the gradients from the 30m datasets at scales of r > 4500 au (see right top and bottom panels on Fig. 3) have a PA of −8 − 14° and are in the opposite to the direction of the bipolar outflows (see Table 2). Moreover, the integrated intensity at these scales seems to trace an elongated structure from east to west that could be compatible with part of the filament.

The panel a of Fig. 16 shows the PV_rot diagram of L1448-2A built from the velocity gradients observed at scales of r < 1300 au. The index of the fitting by a power-law (α ~ −0.9, see Table 3) is consistent with an infalling and rotating protostellar envelope. Therefore, we constrain the radial distribution of the specific angular momentum of L1448-2A at radii of 60−1300 au (see Fig. H.2).

H.2 L1448-NB

L1448-NB (or L1448-IRS3) is a Class 0 protostar in the L1448N complex in the Perseus cloud at (293 ± 20) pc (Ortiz-León et al. 2018). This source is part of a multiple system: L1448-NA, L1448-NB, and L1448-NC. NA and NBare separated by 1700 au, and NB and NC are separated by 4900 au (Maury et al. 2019). L1448-NB is itself also a multiple system: Maury et al. (2019) find a binary system with a separation of ~140 au. ALMA 1.3 mm observations with an angular resolution of ~0.2′′ showed L1448-NB as a triple protostar system with binaries separated by ~60 au and a third core separated by ~180 au from others (Tobin et al. 2016a).

Figures H.3 and 4 show the integrated intensity and centroid velocity maps obtained for L1448-NB from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission respectively. The C¹⁸O emission from the PdBI dataset does not peak on the main protostar L1448-NB1 but on the secondary NB2 resolved by 1.3 mm continuum emission of the binary system (see Table 1 and Fig. H.3). Thus, the origin of the coordinate offsets is chosen to be the secondary protostar NB2 of the binary system (RA: 03^h25^m36^s.315, Dec.: 30°45′15′′.15; see Table 1) to study the kinematics in the protostellar envelope of L1448-NB.

The gradients observed in the PdBI and combined velocity maps of the C¹⁸O emission (see bottom left and middle panels on Fig. 4) are along a northeast-southwest axis and have a PA of Θ ~ 50°, namely a difference >50° with respect to the equatorial axis (see Table 2). Therefore, these gradients at r < 1000 au are a merging of the kinematics between several mechanisms: ejection by outflows, orbital motions of the binary system NB1 and NB2, and rotational motions of the envelope along the equatorial axis. They are consistent with those observed from SMA observations of C¹⁸O emission by Lee et al. (2015) and Yen et al. (2015a). In the outer envelope at r > 1000 au, velocity gradients are in the direction of the bipolar outflows. They have an angle of 64° ≤ ΔΘ ≤ 90° with respect to the equatorial axis (see Table 2). Therefore, these gradients are not dominated by rotational motions along the equatorial axis but by ejection motions from bipolar outflows or due to an external contamination.

The panel b of Fig. 16 shows the PV_rot diagram of L1448-NB built from the velocity gradients observed at scales of r < 1700 au. The index of the fitting by a power-law (α ~ −0.9, see Table 3) is consistent with an infalling and rotating protostellar envelope. Therefore, we constrain the radial distribution of the specific angular momentum of L1448-NB at radii of 150−1700 au (see Fig. H.6).

H.3 L1448-C

L1448-C (or L1448-mm) is located in the Perseus molecular cloud at a distance of (293 ± 20) pc (Ortiz-León et al. 2018). This source was firstly detected as a radio source at 2 cm (Curiel et al. 1990) associated with a strong millimetric continuum emission (Bachiller et al. 1991). This Class 0 protostar (Barsony et al. 1998) has a powerful and very collimated outflow (~70 km s⁻¹; Bachiller et al. 1990, 1995; Guilloteau et al. 1992; Hirano et al. 2010). Podio & CALYPSO (in prep.) estimate the PA of this outflow at −17° (see Table 1). Maret et al. (2020) detect for the first time hints of Keplerian rotation at scales of r ~ 200 au.

Figures H.7 and 5 show the integrated intensity and centroid velocity maps obtained for L1448-C from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission respectively. The gradients observed in the PdBI and combined velocity maps of the C¹⁸O emission (see bottom left and middle panels on Fig. 5) have a PA of −140° < Θ < −121°, namely an angle difference < 30° with respectto the equatorial axis (see Table 2). This gradient is consistent with the one observed in C¹⁸O (2−1) emission with SMA (Yen et al. 2013, 2015b). The 30m datasets show no velocity gradient (see bottom right panel on Fig. 5). From the N₂H⁺ emission (see top panels on Fig. 5), we observed a clear velocity gradient with an orientation consistent with the outflow axis (Θ ≥ 152 °, see Table 2). However, we detect a weak velocity gradient along the equatorial axis. The kinematics seems to be dominated by bipolar jets and outflows or by external contamination from scales >1500 au. Curiel et al. (1999) interpreted the gradient along the equatorial axis and the gradient along the outflows at larger scales as a suggestion that the envelope is rotating and contracting at similar velocities.

The panel c of Fig. 16 shows the PV_rot diagram of L1448-C built from the velocity gradients observed at scales of r < 4000 au. The index of the fitting by a power-law (α ~ −1, see Table 3) is consistent with an infalling and rotating protostellar envelope. From SMA observations, Yen et al. (2013) obtain an index α ~ −1 on scales of 100−3000 au, consistent with our results. Figure H.10 shows the radial distribution of the specific angular momentum of L1448-C at radii of 100−4000 au.

H.4 IRAS2A

IRAS2A (also known as NGC 1333-IRAS2A) is located in the molecular cloud NGC 1333 in the Perseus complex at (293 ± 20) pc (Ortiz-León et al. 2018). This Class 0 protostar was firstly detected via observations of dust continuum emission at 450 μm, 850 μm, and at 3 mm (Sandell et al. 1994; Sandell & Knee 2001; Jørgensen et al. 2004). The CO emission showed the presence of two bipolar outflows in the north-south east-west direction (Knee & Sandell 2000), suggesting the presence of a binary system. Dust continuum observations with VLA resolve IRAS2A into a protobinary system separated by ~140 au (Tobin et al. 2015), which is not resolved by the CALYPSO observations (Maury et al. 2019). Codella et al. (2014a) and Podio & CALYPSO (in prep.) detect jet associated with the principal component with a PA of +205° and a tentative SiO and SO jet emission associated with a possible secondary with a PA = −65°. In this study of the kinematics, we only used the first value (see Table 1).

H¹³CN observations from SMA with an angular resolution < 1′′ (~200 au) do not showorganized velocity gradients or hints of Keplerian motions (Brinch et al. 2009). CALYPSO methanol observations suggest the presence of a weak velocity gradient oriented in the direction perpendicular to the outflow axis, consistent with low rotational motions of the inner envelope (Maret et al. 2014).

Figures H.11 and 6 show the integrated intensity and centroid velocity maps obtained for IRAS2A from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission respectively. The gradients observed in the PdBI and combined velocity maps of the C¹⁸O emission (see bottom left and middle panels on Fig. 6) have an angle difference of ΔΘ ≤ 20° with respectto the equatorial axis (see Table 2). These gradients seem to be due to rotational motions of the envelope slightly contaminated by the outflows or the orbital motions of the multiple system. The other panels on the Fig. 6 show velocity gradient with a PA along a northwest southeast axis and a angle difference of −20° < Θ < −20°. These gradients are not due to rotational motions of the envelope but motions from larger scales. Moreover, the integrated intensity trace a filamentary structure which are consistent with the filament in which the source is embedded.

The panel d of Fig. 16 shows the PV_rot diagram of IRAS2A built from the velocity gradients observed in C¹⁸O emission at scales of r < 1500 au. The index of the fitting by a power-law (α ~ −0.7, see Table 3) is consistent with a Keplerian rotation. However, we also obtain a good reduced χ² when we fit the PV_rot diagram by an infalling and rotating envelope by fixing the index of the power-law at −1 (~0.9, see Table 3). Thus, for this source, the CALYPSO dataset only allow us to estimate a range of the power-law index between −1 and −0.7 (see Table 3) at the scales of 80−1500 au probed in our analysis. Moreover, Keplerian rotation is not detected at the smaller envelope radii investigated by Maret et al. (2020). Thus, the Keplerian rotation due to a large disk is not a robust interpretation for this source. Figure H.14 shows the radial distribution of the specific angular momentum of IRAS2A at radii of 80−1500 au.

H.5 SVS13-B

SVS13 is a multiple system located in the molecular cloud NGC 1333 in the Perseus complex at a distance of (293 ± 20) pc (Ortiz-León et al. 2018). This system is composed of three main sources which are called SVS13-A, B, and C, aligned in the northeast southwest direction (Looney et al. 2003). SVS13-A is a Class I protostar and SVS13-B is a Class 0 protostar. These sources are gathered in a big filamentary structure observed at 450 μm and 1.3 mm (Chandler & Richer 2000; Hull et al. 2014). Observations at 70 μm detected SVS13-A and C but not SVS13-B, which suggests that the latter is deeply embedded in the filamentary structure (Chen et al. 2009). SVS13-A has a powerful outflow observed in the northwest southeast direction by Bachiller et al. (2000); Plunkett et al. (2013). SVS13-B also harbors a very collimated outflows (Bachiller et al. 1998). Podio & CALYPSO (in prep.) estimate the PA of the outflows at 167° (see Table 1).

A velocity gradients with a magnitude of 28 km s⁻¹ pc⁻¹ and symmetric with respect to SVS13-A and B is detected in N₂H⁺ emission from PdBI observations. This gradient suggests that the binary system is physically linked (Chen et al. 2009). A study of dust continuum emission at 8 mm by Segura-Cox et al. (2016) suggests a Keplerian disk with a radius ≲25 au. A similar study from CALYPSO observations at 1.3 mm also suggests an unresolved Keplerian disk with a radius < 60 au (Maury et al. 2019).

Figures H.15 and 7 show the integrated intensity and centroid velocity maps obtained for SVS13-B from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission respectively. The C¹⁸O emission is weakly detected (~5σ) from PdBI observations for this source: the emission is dominated by its companion SVS13-A. The channel maps on Fig. H.16 allow us to constrain the PV_rot diagram with one point at r ~100 au in the envelope (see panel e of Fig. 16). We observed a velocity gradient along the equatorial axis (v − v_sys < 0.4 km s⁻¹, see panel e of Fig. 16) in the combined velocity maps of the C¹⁸O emission (r < 500 au). From the PdBI and combined velocity maps in N₂H⁺ emission (see top left and middle panels on Fig. 7), the velocity gradient is dominated by the intrinsic velocity of the protostar SVS13-A. At larger scales (r > 4500 au), the velocity gradients have an angle difference ΔΘ ≥ 70° with respect to the equatorial axis (see Table 2). These gradients are not due to rotational motions of the envelope but motions from larger scales. Moreover, the integrated intensity trace a filamentary structure which are consistent with the filament in which the source is embedded (see right panels on Fig. 7).

The panel e of Fig. 16 shows the PV_rot diagram of SVS13-B built from the velocity gradients observed in C¹⁸O emission at scales of r < 500 au. The index of the fitting by a power-law (α ~ −0.9, see Table 3) is consistent with an infalling and rotating envelope. Figure H.18 shows the radial distribution of the specific angular momentum of SVS13-B at radii of 100−500 au.

H.6 IRAS4A

The class 0 protostar IRAS4A is located in the molecular cloud NGC 1333 in the Perseus complex at (293 ± 20) pc (Ortiz-León et al. 2018), in the vicinity of another young multiple system, IRAS4B (Lay et al. 1995). It harbors a ~ 2′′ binary system (Looney et al. 2000; Girart et al. 2006; Jørgensen et al. 2007; López-Sepulcre et al. 2017; Maury et al. 2019) and the dust continuum emission at (sub)millimeter wavelengths is dominated by the main protostar IRAS4A1. Santangelo et al. (2015) and Podio & CALYPSO (in prep.) estimate the PA of the outflows at 180° (see Table 1).

Figure H.20 shows the integrated intensity obtained for IRAS4A from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission respectively. C¹⁸O emissionfrom the PdBI dataset is centered on the main protostar IRAS4A1 (RA: 03^h29^m10^s.537, Dec.: 31°13′30′′.98; see Table 1 and Fig. H.20). This is why, the origin of the coordinate offsets is chosen to be the main protostar to study the kinematics in the protostellar envelope of IRAS4A. We also noticed a hole in the integrated intensity and column density (see Fig. H.19) that suggests the C¹⁸O emission may be optically thick at scales of r ≲ 100 au.

The gradients observed in the PdBI and combined velocity maps of the C¹⁸O emission (see bottom left and middle panels on Fig. 8) have an angle difference of ΔΘ ≤ 11° with respectto the equatorial axis (see Table 2). These gradients are consistent with those observed in C¹⁷O emission from SMA observations (Ching et al. 2016). We notice that the gradient is aligned along and symmetric with respect to the binary system IRAS4A1 and IRAS4A2. This could suggest that the binary system is physically linked by orbital motions.

Belloche et al. (2006) find two velocity components at r ~ 6000 au by analyzing CS, C₃₄S and, N₂H⁺ lines obtained with the 30m. This suggests that IRAS4A is collapsing with a velocity of v < 7.30 km s⁻¹ triggered by a fast external compression with a velocity of v > 7.30 km s⁻¹. Figure 7 shows the centroid velocity maps obtained for SVS13-B from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission respectively. To minimize contamination by external cloud compression, the velocity maps were constructed by fitting and removing a second gaussian component fixed at 7.7 km s⁻¹.

The channel maps do not exhibit an organized spatial distribution of velocities along the equatorial axis : the central emission fit show a position angle >|45°| with respect to the equatorial axis, suggesting a possible contamination by the outflow kinematics (see Fig. H.21). To minimize this contamination and probe rotational motions in the equatorial axis, we used the method in the image plane in the PdBI dataset instead of working in the visibilities to constrain in the inner envelope (200−800 au).

The velocity gradients observed in the PdBI velocity maps of the N₂H⁺ emission (see top left panels on Fig. 8) have a direction of Θ ~ −69° consistent with the equatorial axis. However, at larger scales (r > 1000 au), the N₂H⁺ emission is dominated by external cloud compression. The gradient have a direction of Θ ≥ 32°. Moreover, the integrated intensity of N₂H⁺ emission from the 30m datasets (see right top panel on Fig. H.20) trace the filamentary structure which are consistent with the filament in which the source is embedded.

The panel l of Fig. 16 shows the PV_rot diagram of IRAS4A built from the velocity gradients observed at scales of r < 1600 au. This is the only source with increasing velocities from small to large scales. The index of the fitting by a power-law (α ~ 0.8, see Table 3) could be interpreted by solid-body rotation, but as the velocity gradient is not uniform as expected, this is not a robust interpretation for this source (see Sect. 5.1). Moreover, by considering only point at r < 550 au, the index of the fitting by a power-law is consistent with an infalling and rotating envelope (α ~ −1.3, see Table 3). We did not take this source into account to build the radial distribution of the specific angular momentum due to rotation motions.

H.7 IRAS4B

IRAS4B (orNGC 1333-IRAS4B) is a Class 0 protobinary system in the Perseus molecular cloud at a distance of (293 ± 20) pc (Ortiz-León et al. 2018). The binary companions 4B1 and 4B2 are separated by ~11′′ (Looney et al. 2000; Jørgensen et al. 2007). The source harbors a young bipolar outflow almost aligned along a north-south direction (Choi 2001). Podio & CALYPSO (in prep.) distinguish two outflows: a first associated with the main protostar 4B1 at 167° and a second with a PA of −99° (see Table 1). IRAS4B is in the vicinity of IRAS4A which undergo a collapse triggered by a fast external compression with a velocity of v > 7.30 km s⁻¹ (Belloche et al. 2006). Thus, collapse of IRAS4B could have triggered in the same way. It is why, to minimize contamination by external cloud compression, the velocity maps were constructed by fitting and removing a second Gaussian component fixed at 7.7 km s⁻¹ (see Fig. H.22). Previous studies of C¹⁸O (2−1) emission from SMA observations are detected no organized velocity gradients consistent with rotational motions (Yen et al. 2013, 2015b).

Figure 9 shows the centroid velocity maps obtained for IRAS4B from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission. The gradients observed in the PdBI and combined velocity maps of the C¹⁸O emission (see bottom left and middle panels on Fig. 9) have an angle difference of ΔΘ ≤ 20° with respect to the equatorial axis (see Table 2). The velocity gradients observed in the PdBI velocity maps of the N₂H⁺ emission (see top left panels on Fig. 9) is also consistent with the equatorial axis.

However, the channel maps do not exhibit an organized spatial distribution of velocities along the equatorial axis: the central emission fit show a position angle >|45°| with respect to the equatorial axis, suggesting a possible contamination by the outflow kinematics (see Fig. H.23). To minimize this contamination and probe rotational motions in the equatorial axis, we used the method in the image plane in the PdBI dataset instead of working in the visibilities to constrain in the inner envelope (200−500 au).

At outer envelope scales (r > 1000 au), the N₂H⁺ emission shows a reversed gradient dominated by external cloud compression. Moreover, the integrated intensity of N₂H⁺ emission from the 30m datasets (see right top panel on Fig. H.22) trace the a filamentary structure which are consistent with the filament in which the source is embedded.

The panel g of Fig. 16 shows the PV_rot diagram of IRAS4B built from the velocity gradients observed at scales of r < 1100 au. The index of the fitting by a power-law (α ~ −0.6, see Table 3) is consistent with a Keplerian rotation. However, we also obtain a good reduced χ² when we fit the PV_rot diagram by an infalling and rotating envelope by fixing the index of the power-law at −1 (~0.6, see Table 3). Thus, for this source, the CALYPSO dataset only allow us to estimate a range of the power-law index between −1 and −0.6 (see Table 3) at the scales of 170−1100 au probed in our analysis. Moreover, Keplerian rotation is not detected at the smaller envelope radii investigated by Maret et al. (2020). Thus, the Keplerian rotation due to a large disk is not a robust interpretation for this source. Figure H.25 shows the radial distribution of the specific angular momentum of IRAS4B at radii of 170−1100 au.

H.8 IRAM04191

At a distance of 140 pc (Torres et al. 2009), IRAM04191 + 1522 (hereafter IRAM04191) is located in the southern part of the Taurus cloud in the viciny of the Class I protostar IRAS 04191. Its envelope mass (0.5 M_⊙, André et al. 2000), its low luminosity (0.1 L_⊙, Dunham et al. 2006), and its temperature suggest that is one of the youngest accreting and isolated Class 0 sources known in Taurus (t ~1−3 × 10⁴ yr, André et al. 1999). This source harbors a strongly collimated outflow with a PA = 28° (Belloche et al. 2002) and an inclination angle to the line of sight estimated at ~40° (André et al. 1999; Maret et al. 2014). A new study by Podio & CALYPSO (in prep.) estimates the PA of the outflows at 20° (see Table 1). We used this later value for our study. From CALYPSO observations, Maury et al. (2019) estimate an envelope mass of 0.45 M_⊙ and did not resolve candidate disk down to <50 au.

The 30m observations of several molecular lines suggests signatures of the envelope infalling motions with an estimated velocity of ~0.2 km s⁻¹ to r ~1000 au and ~0.1 km s⁻¹ to r ~ 10000 au (Belloche et al. 2002). These observations associated with radiative transfer model allowed kinematics of the outer envelope to be constrained. A strong velocity gradient of ~40 km s⁻¹ pc⁻¹ was detected at scales of ~3500 au along the northwest-southeast direction, namely in the equatorial axis. At ~1500 au, Lee et al. (2005) detect a velocity gradient of ~100 km s⁻¹ pc⁻¹. These gradient are interpreted as rotational motions of the protostellar envelope (Belloche et al. 2002; Takakuwa et al. 2003).

Figures H.26 and 10 show the integrated intensity and centroid velocity maps obtained for IRAM04191 from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission, respectively. The C¹⁸O emission is weakly detected (~5σ) from PdBI observations for this source which is one of the lowest luminosity sources in our sample (see bottom left panel on Fig. 10). However, the combined map allows the kinematics of the protostellar envelope between 50 and 1000 au to be constrained. At those scales (see bottom middle panel on Fig. 10), we noticed a velocity gradients with a direction of Θ = −83°, consistent with the equatorial axis. At r > 1000 au, the velocity gradients observed in the velocity maps of the N₂H⁺ emission (see top panels on Fig. 10) have a direction of Θ ~ 100°, consistent with those detected by Belloche et al. (2002) at similar scales. Therefore, there is a reversal of the kinematics between inner (~500 au) and outer envelope scales (>1000 au).

The panel h of Fig. 16 shows the PV_rot diagram of IRAM04191 built from the velocity gradients observed at scales of 50−800 au. The index of the fitting by a power-law (α ~ −0.3, see Table 3) is consistent with a Keplerian rotation. However, we also obtain a good reduced χ² when we fit the PV_rot diagram by an infalling and rotating envelope by fixing the index of the power-law at −1 (~1.1, see Table 3). Thus, for this source, the CALYPSO dataset only allow us to estimate a range of the power-law index between −1 and −0.3 (see Table 3) at the scales of 50−800 au probed in our analysis. Moreover, Keplerian rotation is not detected at the smaller envelope radii investigated by Maret et al. (2020). Thus, the Keplerian rotation due to a large disk is not a robust interpretation for this source. Figure H.29 shows the radial distribution of the specific angular momentum of IRAM04191 at radii of 50−800 au.

H.9 L1521F

At a distance of 140 pc (Torres et al. 2009), L1521F (or MC27) is one of the densest cores in Taurus cloud. This source was firstly identified as a starless core (Codella et al. 1997; Onishi et al. 1999; Crapsi et al. 2004). Its very low luminosity (L < 0.004; Maury et al. 2019) indicates that it is among the youngest known Class 0 protostars (André et al. 2000) and still may preserve the initial conditions of star formation (Terebey et al. 2009). The emission of the molecular line HCO⁺ showing a high central density and an asymmetric signature typical of the collapse suggests that the source is in the early stages of gravitational collapse (Onishi et al. 1999). Its low luminosity of L_bol < 0.04 (Maury et al. 2019; Ladjelate et al. in prep.) indicating that is one of the known Class 0 youngest prostars (André et al. 2000), L1521F could still preserve its initial conditions (Bourke et al. 2006; Terebey et al. 2009). N₂H⁺ and NH₃ emissions from VLA observations show a flattened structure perpendicular to the outflows (Tobin et al. 2010). Podio & CALYPSO (in prep.) do not detect clear jet signature in SiO and SO but only blue and red-shifted clumps in CO. In our study of the kinematics, we used the PA of the outflows at 240° observed in CO and HCO⁺ by Tokuda et al. (2014, 2016) (see Table 1). Recent ALMA observations of ¹²CO emission showed a clear velocity gradient in the equatorial axis consistent with the presence of a disk of radius of ~10 au surrounded a stellar embryo of ~0.2 M_⊙ (Tokuda et al. 2017).

Figure H.30 show the centroid intensity maps obtained for L1521F from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission. The C¹⁸O emission is weakly detected (~5σ) from the PdBI observations for this source which is one of the lowest luminosity sources in our sample (see bottom left panel on Fig. 11). The emission peak does not correspond to the continuum emission peak at 1.3 mm determined from the PdBI dataset by Maury et al. (2019) (see Table 1). This is consistent with the existence of a starless high-density core, MMS-2 (RA:4^h28^m38^s.89, Dec.:+26°51′33′′.9) detected by ALMA Cycle 0 observations in dust continuum and H¹³CO(3–2) emissions (Tokuda et al. 2014). The protostar MMS-1 (RA:4^h28^m38^s.96, Dec.: +26°51′35′′) is also not spatially resolved in H¹³CO(3–2) emission. The L1521F protostar is not spatially resolved by the emission H¹³CO (3−2) but a clear velocity gradient is observed along the northwest southeast axis at outer envelope scales consistent with the equatorial axis (Tokuda et al. 2014). This velocity gradient is in the same direction as the elongated structure we observed in combined N₂H⁺ emission (see top middle panel on Fig. 11) and also observed by Tobin et al. (2010). The absence of velocity gradient at small scales (r < 1500 au) is interpreted as due to low dust temperature of T ~ 10K (Tokuda et al. 2014, 2016) which is not sufficient to excite high density molecules (E_upper = 15.8 K for the molecular transition C¹⁸O (2−1); CDMS).

Figures 11 show the centroid velocity maps obtained for L1521F from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission. Only the 30m velocity map of the C¹⁸O emission (see bottom right panel on Fig. 11) shows a gradient with a blue- and a red-shifted velocity components along the equatorial axis at radii ~ 3000 au. The direction of this gradient (Θ ~ −8°) is consistent with the equatorial axis axis and the gradients previously observed by Tokuda et al. (2014). However, we noted that the integrated intensity of the C¹⁸O emission from the 30m datasets does not trace the envelope but a more complex and diffuse structure.

The panel i of Fig. 16 shows the PV_rot diagram of L1521F built from the velocity gradients observed at scales of 1500−4200 au. The index of the fitting by a power-law is close to 0 (see Table 3). However, we also obtain a quite good reduced χ² when we fit the PV_rot diagram by an infalling and rotating envelope by fixing the index of the power-law at −1 (~1, see Table 3). Thus, for this source, the CALYPSO dataset only allow us to estimate a range of the power-law index between -1 and 0 (see Table 3) at the scales of 1500−4200 au probed in our analysis. Figure H.33 shows the radial distribution of the specific angular momentum of L1521F at radii of 1500−4200 au.

H.10 L1527

At a distance of 140 pc (Torres et al. 2009) in the Taurus molecular cloud, L1527 is the closest Class 0/I source which exhibits a candidate proto-planetary disk with an outer radius of 50−90 au around a 0.2−0.3 M_⊙ central object (Tobin et al. 2012; Ohashi et al. 2014). A recent study by Aso et al. (2017) find an outer radius of 74 au and a central mass of 0.45 M_⊙. The embryo mass is close of the envelope mass (~1 M_⊙), which sets the source at the boundary between the phases Class 0 and Class I. This source is an ideal target for studying the kinematics of the disk and the envelope at high angular resolution because it is seen almost edge-on, with an angle ~90° with respect to the sky plane (Ohashi et al. 1997a). L1527 would also be embedded in a filament with an angle of PA ~ 135° (Marsch & HGBS, in prep).

Figures H.34 and 12 show the integrated intensity and centroid velocity maps obtained for L1527 from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission, respectively. At scales of r < 2000 au (see left and middle panels on Fig. 12), show a clear velocity gradient along the north-south axis. These velocity gradients are consistent with those observed in C¹⁸O (2−1) emission with CARMA, SMA and ALMA (Tobin et al. 2012; Yen et al. 2013, 2015b; Ohashi et al. 2014). The linear gradient fits give orientations of −9° < Θ < 26° which is consistent with the equatorial axis direction assumed perpendicular to the outflows (see Table 2). Velocity maps from the 30m datasets at scales ~3000−5000 au (see right panels on Fig. 12) show gradients along the east-west axis opposite to the outflow direction, with an orientation ~113−123° (see Table 2).

The panel j of Fig. 16 shows the PV_rot diagram of L1527 built from the velocity gradients observed at scales of r < 2000 au. The index of the fitting by a power-law is consistent with a rotating and infalling envelope (α ~ −1.1, see Table 3). This result is consistent with previous studies of Yen et al. (2013) and Ohashi et al. (2014) which are found indices of approximately −1 and approximately −1.2 on scales of 40−500 au and 100−1500 au, respectively.Ohashi et al. (2014) constrain the specific angular momentum at j ~ 6 ×10⁻⁴ km s⁻¹ pc at radii of 100−1600 au. This value is consistent with the value we found on the Fig. H.37 which shows the radial distribution of the specific angular momentum of L1527 at radii of 40−2000 au.

H.11 L1157

In the Cepheus molecular cloud at a distance previously estimated to be 250 pc (Looney et al. 2007) but determined at (352 ± 18) pc by recent Gaia parallax measurements (Zucker et al. 2019), L1157 is a Class 0 protostar deeply embedded within a large circumstellar envelope (Gueth et al. 2003; Beltrán et al. 2004). The source harbors a powerful asymmetrical bipolar outflows at large scales (Tafalla & Bachiller 1995; Zhang et al. 1995; Gueth et al. 1996). Podio et al. (2016) and Podio & CALYPSO (in prep.) estimate the PA of the outflows at 163° (see Table 1). Previous study of C¹⁸O (2−1) observations with SMA showed no clear sign of rotational motions but kinematic models predict L1157 exhibit a possible disk with an outer radius < 5 au (Yen et al. 2015b).

Figures H.38 and 13 show the integrated intensity and centroid velocity maps obtained for L1157 from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission, respectively. Velocity maps of C¹⁸O emission from the PdBI and combined datasets (see bottom left and middle panels on Fig. 13) show gradients with a direction of Θ ~ 13−35° (see Table 2). They are dominated by the kinematics of the outflows but a weak velocity gradient is observed along the equatorial axis at radii r < 2000 au. Velocity maps of N₂H⁺ emission from PdBI and combined datasets (see top left and middle panels on Fig. 13) trace cavities along the outflow axis. These gradients are consistent with those detected in CARMA observations by Chiang et al. (2010) and Tobin et al. (2011). At scales of r > 2500 au, velocity maps from the 30m datasets (right panels on Fig. 13) show a gradient with a direction of Θ ~ −130°, in the opposite direction with respect to small scales observed along the equatorial axis.

The panel k of Fig. 16 shows the PV_rot diagram of L1157 built from the velocity gradients observed at scales of r < 700 au. The index of the fitting by a power-law is close to 0 (see Table 3), consistent with differential rotation of the envelope. At each PdBI emission channel, the central emission fit show a position angle >|45°| with respect to the equatorial axis, suggesting a contamination by the outflow kinematics (see Fig. H.39). To minimize this contamination and probe rotational motions in the equatorial axis, we used the method in the image plane in the PdBI dataset instead of working in the visibilities to constrain in the inner envelope (80−200 au). Figure H.41 shows the radial distribution of the specific angular momentum of L1157 at radii of 80−700 au.

H.12 GF9-2

GF9-2 (L1082C Benson & Myers 1989 or GF9-Core Ciardi et al. 2000) is located in the filamentary cloud GF9-2 (or LDN 1082) at a distance of 200 pc (Wiesemeyer 1997; Wiesemeyer et al. 1998). In this study, we adopt this distance but it is very uncertain and some studies estimated a higher distance between 440−470 pc (Viotti 1969, C. Zucker, priv. comm.) and 900 pc (Reid et al. 2016). Due to its low luminosity (0.3 L_⊙, Wiesemeyer 1997), GF9-2 is has been cataloged as a transitional object between the prestellar and the Class 0 phases (Wiesemeyer et al. 1999). GF9-2 would be the youngest source in our sample and would therefore retain the initial conditions of its gravitational collapse because the central embryo does not have a powerful bipolar flow (Furuya et al. 2006). Podio & CALYPSO (in prep.) detect for the first time the outflows with a PA of 0° (see Table 1).

Figures H.42 and 14 show the integrated intensity and centroid velocity maps obtained for GF9-2 from the PdBI, combined, and 30m CALYPSO datasets for the C¹⁸O and N₂H⁺ emission, respectively. The CALYPSO datasets allow us to constrain the kinematics of this source at the envelope scales for the first time. Only the PdBI and combined maps of the C¹⁸O emission (see bottom left and middle panels on Fig. 14) show a velocity gradient associated with the source. The gradient is consistent with the equatorial axis, with a direction of −160° < Θ < −130° (see Table 2). Moreover, the N₂H⁺ integrated intensity from the 30m datasets at scales of r > 3000 au traces th e filament structure in which the source is embedded.

The panel l of Fig. 16 shows the PV_rot diagram of GF9-2 built from the velocity gradients observed at scales of r < 900 au. The index of the fitting by a power-law is consistent with an infalling and rotating envelope (approximately −0.8, see Table 3). Figure H.45 shows the radial distribution of the specific angular momentum of GF9-2 at radii of 70−900 au.

References

---

All Tables

All Figures