Low dust emissivities and radial variations in the envelopes of Class 0 protostars: possible signature of early grain growth

Context. Analyzing the properties of dust and its evolution in the early phases of star formation is crucial to put constraints on the collapse and accretion processes as well as on the pristine properties of planet-forming seeds. Aims. In this paper, we aim to investigate the variations of the dust grain size in the envelopes of the youngest protostars. Methods. We analyzed Plateau de Bure interferometric observations at 1.3 and 3.2 mm for 12 Class 0 protostars obtained as part of the CALYPSO survey. We performed our analysis in the visibility domain and derived dust emissivity index ( β 1 − 3mm ) proﬁles as a function of the envelope radius at 200–2000 au scales. Results. Most of the protostellar envelopes show low dust emissivity indices decreasing toward the central regions. The decreasing trend remains after correction of the (potentially optically thick) central region emission, with surprisingly low β 1 − 3mm < 1 values across most of the envelope radii of NGC 1333-IRAS 4A, NGC 1333-IRAS 4B, SVS13B, and Serpens-SMM4. Conclusions. We discuss the various processes that could explain such low and varying dust emissivity indices at envelope radii 200–2000 au. Our observations of extremely low dust emissivity indices could trace the presence of large (millimeter-size) grains in Class 0 envelopes, in which case our results would point to a radial increase of the dust grain size toward the inner envelope regions. While it is expected that large grains in young protostellar envelopes could be built via grain growth and coagulation, we stress that the typical timescales required to build millimeter grains in current coagulation models are at odds with the youth of our Class 0 protostars. Additional variations in the dust composition could also partly contribute to the low β 1 − 3mm we observe. We ﬁnd that the steepness of the β 1 − 3mm radial gradient depends strongly on the envelope mass, which might favor a scenario in which large grains are built in high-density protostellar disks and transported to the intermediate envelope radii, for example with the help of outﬂows and winds.


Introduction
Investigating the properties of dust in the early phases of star formation is essential to understand the physico-chemistry taking place in the collapsing core as well as to trace the pristine properties of planet-forming material building up the protoplanetary disks. From star-forming cores down to disk scales, the collisions and interactions between grain particles increase and the physical conditions and chemistry evolve, affecting the grain population, their progressive growth, and reprocessing. Deriving the dust properties from observations is not trivial. Submillimeter and millimeter wavelength observations of the dust thermal emission have been extensively used to put constraints on the dust opacity. In this regime indeed, the opacity has a strong frequency dependence κ ν = κ 0 (ν/ν 0 ) β , with β the so-called dust emissivity index. Analyzing this dust emissivity index and its variations can help us to probe the grain size distribution (Draine 2006;Natta et al. 2007;Testi et al. 2014). Works based on PRONAOS, Herschel, the Balloon-borne Large Aperture Sub-millimetre Telescope (BLAST), or James Clerk Maxwell Telescope (JCMT) measurements have revealed that star-forming clumps present significant variations in their grain emissivity (Chini et al. 1997;Dupac et al. 2003;Martin et al. 2012). No evidence of major grain growth was discovered toward star-forming filaments such as OMC 2/3 (Sadavoy et al. 2016) or the center of pre-stellar cores such as L1544 (Chacón-Tanarro et al. 2017, at least at the angular resolution at which these studies were performed. In some Perseus star-forming clumps, however, Chen et al. (2016) found dust emissivity indices as low as 1.0, while these clumps have values of about 1.6 in the diffuse interstellar medium (ISM; Planck Collaboration XI 2014; Juvela et al. 2015). These results are in line with models predicting that grain coagulation and the formation of aggregates take place in dense molecular clouds (Ossenkopf & Henning 1994) and that the submillimeter and millimeter emissivity increasing with fluffiness (Bazell & Dwek 1990;Ossenkopf 1993). Grain growth processes inside molecular clouds are also required to explain the "coreshine" phenomenon, a scattering process at 3.6 and 4.5 µm that dominates absorption Steinacker et al. 2010Steinacker et al. , 2015Andersen et al. 2013;Lefèvre et al. 2014) and that was first discovered using Spitzer observations.
At the other end of the star formation timeline, observations of pre-main-sequence objects (in particular T Tauri stars) with circumstellar disks reveal β values ranging from 2 down to −1 (Beckwith & Sargent 1991); these values are compatible with the presence of millimeter (Pérez et al. 2012(Pérez et al. , 2015Tazzari et al. 2016) and even centimeter grain sizes (Testi et al. 2003; see also models from D' Alessio et al. 2001). How and when dust grains evolve during the building of the star and disk system, i.e., during the protostellar phase, is less clear. Several analyses, mostly on individual objects, have tried to address these questions. While Agurto-Gangas et al. (2019) have found no clear sign of grains larger than 100 µm in the inner envelope of the Class I protostar Per-emb-50, Jørgensen et al. (2007) showed that the inner envelopes of deeply embedded protostars could present dust emissivity indices of ∼1 that may be indicative of grain growth in these young objects. No evidence seems to be found yet of an evolution of the grain growth from the inner regions of the envelopes to the planet-forming disks (e.g., Ricci et al. 2010;Testi et al. 2014).
Several studies have suggested that grain growth could occur even earlier than the Class I phase, i.e., at the Class 0 phase. Class 0 objects are the youngest protostars: most of their mass resides in an extended envelope (with outer radii as large as 10 4 au) that is being accreted onto a protostellar embryo (André et al. , 2000 during very short timescales (t ≤ 5 × 10 4 yr; Maury et al. 2011). Kwon et al. (2009), Chiang et al. (2012, and Li et al. (2017) showed that a dozen of Class 0 protostellar envelopes present β at millimeter wavelengths lower than 1.7, and for some objects even lower than 1. β values of 1 have been obtained at 1.3 mm from models using grain growth alone (Miyake & Nakagawa 1993;Ossenkopf & Henning 1994;Draine 2006). Many questions remain, however, about whether β values lower than 1 are realistic or about how large (millimeter-size) grains can form in Class 0 envelopes despite the long timescales and high densities required to build them. Current radiative transfer models of young protostellar envelopes and disks already have tried to incorporate variations in the grain composition (e.g., ices) and size distribution (Whitney et al. 2003;Robitaille et al. 2006;Facchini et al. 2017). More observations are, however, needed to put constraints on the dust properties within the protostellar envelopes during the disk-forming stage, that is the Class 0 stage, and to identify the environmental parameters (protostar mass, temperature, and hosting cloud conditions) leading to these variations to probe the pristine material building up the planet-forming disks and build a coherent end-to-end scenario of the grain evolution.
Class 0 protostars are ideal candidates to investigate the potential grain growth processes taking place during the main accretion phase and establish the pristine dust properties from which the future disk originates. In this work, we aim to perform a resolved and multifrequency millimeter analysis of the variations of the dust properties in a sample of 12 Class 0 protostellar envelopes observed as part of the Continuum And Lines in Young ProtoStellar Objects survey (CALYPSO 1 ). Interferometric observations at millimeter wavelengths are ideal to perform this analysis: most of the envelope/disk emission is optically thin in this regime, limiting the biases in accessing the dust opacity and studying its spatial variations. We provide details on the sample and data reduction in Sect. 2 and analyze the 1-3 mm β radial profiles in Sect. 3. We correct our β 1−3mm profiles for the central regions in Sect. 4 and describe the implications of these corrections for our results. We analyze the dust emissivity values and radial gradients and their dependence on the global source properties in Sect. 5. Our results are summarized in Sect. 6. 1 PI: Ph. André; http://irfu.cea.fr/Projets/Calypso/

The sample
The CALYPSO sources were selected from various nearby (d = 140−436 pc) star-forming clouds and observed at 1.3 and 3.2 mm using the Institut de Radioastronomie Millimétrique/Plateau de Bure Interferometer (IRAM/PdBI; Maury et al. 2019). Both the most extended configuration (A array) and the intermediate antenna configuration (C array) were used to probe the dust continuum emission over a wide range of baselines. For this analysis, we selected the 12 brightest protostars of the CALYPSO sample, namely sources with a 1.3 mm PdBI integrated flux higher than 100 mJy (according to Table 4 in Maury et al. 2019). The sample is diverse in terms of morphologies (including single objects as well as close and wide binaries), bolometric luminosities (from 1 up to 30 L ), or bolometric temperatures (from 25 to 60 K). The source characteristics are summarized in Table 1. For further details on the morphology of the dust continuum emission in this sample of sources, we refer to the maps shown in Appendix B of Maury et al. (2019). Maps of the gas kinematics are shown in Gaudel et al. (2019).

Data processing
We worked on the common range of baselines sampled at 1.3 and 3.2 mm, thus ranging from 20 to 220 kλ. The data calibration (and self-calibration for some of the sources) was performed with GILDAS 2 using the Continuum and Line Interferometer Calibration (CLIC) package. Maury et al. (2019) have provided details on the observations, data reduction, and Plummer models fitted to the dust continuum visibility profiles. Figure 1 shows the real part of the visibilities at 1.3 and 3.2 mm (circles and triangles, respectively); the baselines are averaged every 20 kλ. Errors are estimated by adding in quadrature the calibration uncertainties (∼10% at 94 GHz and ∼15% at 231 GHz) and the rms dispersion around the weighted mean in each bin. Our analysis is restricted to the properties of the dust continuum emission around the primary protostars, as defined in Table 2 of Maury et al. (2019). The two sources L1448-N and L1448-2A host a secondary source that greatly contributes to the visibility profiles. For these sources, we chose to show only the amplitudes of the shortest baselines derived using the center of the binary system as the phase center and thus tracing the general dust properties in the outer envelope of the binary system. We additionally indicated the fluxes at the ∼200 kλ baseline obtained using the primary source as the phase center to show how the dust properties connect between these two scales.

Free-free + synchrotron contribution
In protostars, thermal (free-free) emission produced by a jet or a photo-evaporative wind (see Anglada et al. 1996;Lugo et al. 2004;Pascucci et al. 2012;Tychoniec et al. 2018a, for observations and theoretical models) and to a lesser extent, from nonthermal (synchrotron) emission (André 1996) originating from the presence of magnetic fields, can contribute to the 1-3 mm emission within the central 50 au around the protostar. We estimated the contribution to the 1.3 and 3.2 mm fluxes not coming from thermal dust emission from radio measurements beyond 2 cm. Details for each source are provided in Appendix A. The contribution is systematically lower than 10% at 1.3 mm, but can be non-negligible at 3.2 mm. Correcting for  (André et al. 2010;Ladjelate, in prep.) and rescaled to the distances used in this analysis. ( f ) Envelope mass. The M env were rescaled to the distances used in this analysis. (g) A mean distance of 436 pc for the Serpens Cloud has been determined by Ortiz-León et al. (2018b) using the Gaia-DR2 results and seems to correspond to the average distance determined from the Main cloud or the star-forming region W40. Ongoing work reanalyzing the Gaia data toward Serpens South seems, however, to suggest a first layer of extinction around 300-350 pc while its associated young stellar objects could be distributed in a larger extinction layer up to distances of 350 pc (Palmeirim et al., in prep.). We use this new distance for the paper. As we are working on flux ratios, this does not affect the slope we present for SerpS-MM18. This, however, has an impact on the envelope mass discussed in Sect. 6, which would be higher by a factor of 1.6 if the distance is 436 pc. this compact contamination is one motivation for the "central region correction" (c.r.c) we apply and discuss in Sect. 4.

Optically thick emission at 1.3 and 3.2 mm
If regions of our Class 0 envelopes are optically thick, the intensity should only depend on the temperature of the emitting material and the brightness temperature should be equal to the actual temperature in the region. Following Motte et al. (1998), we have where S is the flux density (here taken at 1.3 mm), Ω the solid angle associated with the physical scale probed, τ the optical depth, and B the Planck function derived from the dust temperature. The dust temperature is estimated with Eq. (6) (cf. Sect. 3). If the medium is optically thick and if tau > 1, (1-exp(−τ))*B tends toward B. If the medium is optically thin, we should observe S/Ω < B(T ). We thus use the radially determined fluxes at 231 GHz maps divided by Ω and compare them to the B(T ) profiles. Figure 2 shows how the two quantities vary with the scales in astronomical units. The absence of overlap between the two profiles indicates that most of the continuum emission of the envelope scales we probe is optically thin, except perhaps for the 1.3 mm emission coming from the <200 au regions of IRAS 4A, IRAS 4B, or SerpM-SMM4. In this case again the central region correction applied in Sect. 4 allows us to correct for the potential central optical thickness effects affecting the visibility profiles.

Observational results
An interferometer only samples the uv plane at discrete spatial frequencies. Voids in this under-sampled Fourier space are a challenge for interferometric image reconstruction. We thus decided to perform our analysis of the grain size distribution variations directly in the visibility domain. As described in many analyses (Harvey et al. 2003;Maury et al. 2019, among others), if the density and temperature distributions in the envelope follow a power law, ∝ r −p and ∝ r −q (Adams 1991;Motte & André 2001), respectively, then in the Rayleigh-Jeans approximation, the intensity will follow a power law I(r) ∝ r −(p+q−1) . The resulting visibility distribution will be a power law as well, expressed as V(b) ∝ b (p+q−3) . The ratio of the amplitude visibility measured at two wavelengths at a common baseline thus traces the spectral index of the dust composing the emitting material within a given area corresponding to the scale probed by the specific baseline. The variations of the millimeter spectral index β mm can then be used as a signature of variation of the grain size distribution with lower β values associated with the presence of more emissive grains, as expected for larger grains (Kruegel & Siebenmorgen 1994).

Observed spectral index profiles with uv distance
The dust continuum emission fluxes at 1.3 and 3.2 mm, as measured in bins of baselines, are shown in Fig. 1 (2)) are overlaid (blue line). Their values are shown on the right y-axis. The error bars indicate the flux uncertainties and include the calibration errors. These errors are used to derive the uncertainties on α 1−3mm shown as the shaded blue area. As L1448-N and L1448-2A host a secondary source, we only show the continuum flux densities at the shortest baselines but also indicate the flux at ∼200 kλ obtained using the primary source as the phase center (see Sect. 2.2 for more details).
also shows the variations of the observed spectral index α 1−3mm defined as where ν 0 and ν 1 are equal to 231 and 94 GHz, respectively. Figure 1 indicates, overlaid on the visibilities, the variations of the observed spectral index α 1−3mm (blue lines, scale provided on the right axis) for each object.
In order to express the uv distances in terms of angular separations, we use the relation linking the brightness distribution I of an object to its complex visibility given by the Fourier transform where (x, y) represent the angular coordinates on the sky (in radians) and (u, v) the coordinates describing the baseline.
For an assumed circularly symmetric object, we can use the zeroth-order Bessel function J 0 to write where x 0 is the first zero of J 0 . The visibility reaches its first zero when b V = 0 = 1.22(λ/θ), resulting in a baseline-angular scale relation given by These angular scales θ, over which the dust continuum emission is integrated to compute the spectral index are indicated on the top axis of Fig. 1. We observe that most of the objects show a clear decrease of the observed spectral index α 1−3mm toward larger uv distances, thus smaller envelope radii. Several objects of the sample have already been studied using other millimeter interferometric observations with similar resolutions, allowing us to confront our observations with previous works. L1157 has been studied using CARMA by Kwon et al. (2009) andChiang et al. (2012) who derived the profile of the dust emissivity index as a function of uv distance using the Rayleigh-Jeans approximation and the optically thin assumption. Our results are consistent with the relatively flat profile they observe on 10-40 au and 25-2500 au scales, respectively: L1157 has one of the weakest gradient of our sample. For L1527, we find a very similar flat profile of the observed spectral index as a function of uv distance than that found by Tobin et al. (2013).
The observed spectral indices derived in baseline bins as a function of their associated 3.2 mm flux densities are shown in

Dust emissivity index radial profiles
In the Rayleigh-Jeans approximation (and assuming optically thin emission at millimeter wavelengths), α values translate to a dust emissivity index β = α-2. This is, however, only a proxy for β as the Rayleigh-Jeans approximation does not hold across protostellar envelopes because of their low dust temperatures. To offer a more robust approximation, we recalculate β 1−3mm assuming that our objects are centrally illuminated spherical envelopes in which the temperature varies with radius (Terebey et al. 1993;Motte & André 2001) as The protostellar internal luminosities L int (in L ) are tabulated in Table 1. We then use the temperature profiles to calculate the binned β values using where ν 0 and ν 1 are equal to 231 and 94 GHz, respectively, and Emissivity index β . Relation between the dust emissivity index β 1−3mm and the observed spectral index α 1−3mm -2. The dashed line represents the 1:1 relation, while the plain line represents the polynomial regression to the data (β 1−3mm = −0.07 (α 1−3mm −2) 2 + 1.12 (α 1−3mm −2) + 0.10). Each symbol represents a different protostellar envelope and symbols are color coded with respect to our modeled dust temperature.
in Fig. 4. We immediately observe that the relation deviates from the 1:1 relation as one points toward colder regions (with α 1−3mm -2 underestimating β at colder temperatures). This has a direct impact on the β 1−3mm radial gradients, steepening most of them. For a visual illustration, maps of β 1−3mm have been generated from the 1.3 and 3.2 mm dust continuum cleaned maps (see Table C.1 and Appendix C) but we performed the following analysis on the visibility datasets rather than on the maps. Figure 5 (black lines) shows how the dust emissivity index β 1−3mm varies within the envelopes for each object of our sample. Our β 1−3mm values can be reconnected with those presented in Bracco et al. (2017) on larger scales derived with Herschel and the Neel-IRAM-KID-Array (NIKA). These authors found that the emissivity β decreases radially in the two protostellar cores studied in Taurus while staying flat in the pre-stellar core they analyzed; values of β = 1.0-1.5 nicely match the range of β 1−3mm we observe in most of our sources beyond 2000 au. We provide the β 1−3mm values we derive at 500 au in Table 2. These values range from 0.38 to 0.92. We note that the typical value found in the diffuse ISM is β ∼ 1.6 (Planck Collaboration XI 2014; Juvela et al. 2015). The β 1−3mm values observed at 500 au are thus extremely low for the whole sample. If absolute flux calibration affects the β 1−3mm measured, the fact that all sources present low β 1−3mm values would imply a systematically underestimated 1.3 mm flux (or overestimated 3.2 mm flux) that is not expected since the data comes from >30 independent tracks observed over a period of less than years. Phase noise could produce a loss of 1.3 mm flux at long baselines, but this does not affect our analysis since we limit ourselves to baselines <200 kλ, thus recovering the dust continuum emission at envelope scales.
Clear radial variations appear in the envelopes of several objects. To quantify these variations, we fit the β 1−3mm radius trend in lin-log scale and call "β 1−3mm gradient" the slope of the relation linking beta with the logarithm of the radius (in au) as follows:  Table 2). The colored lines indicate the trend in various angular bins with respect to the outflow direction. If PA is the position angle with respect to the outflow axis, the red lines indicate the trends for regions located close to the outflow direction (  where A is the β 1−3mm gradient. With this definition, A = 2 indicates that β 1−3mm doubles its value over an order of magnitude in radius. Gradients are calculated using baselines below 100 kλ to limit the noise present at longer baselines and only probe the dust emissivity variations at envelope scales (i.e., 300-4000 au, depending on the source distance). We provide the β 1−3mm gradients in Table 2. Some objects show flat β 1−3mm profiles, in particular IRAS 4B, SVS13B, and L1157. In L1157, Kwon et al. (2015) found that the dust emissivity index changes with radius and has a mean β 1−3mm value of 0.76, but values as low as 0.1 at 80 kλ. We do not detect such low values at these scales. Our results are closer to those of Chiang et al. (2012) in that respect, namely that β 1−3mm does not vary much depending on the envelope radius in L1157.

Investigating dust emissivity variations with respect to source morphology
The source physical properties (high densities in the equatorial plane, outflows) could have an impact on the grain populations and thus affect the structure and evolution of the central regions. To analyze potential local variations of the dust emissivity index, we chose to analyze the visibilities corresponding to different angular bins separately. The outflow opening angles are still poorly constrained and can vary significantly from one source to another. In Hsieh et al. (2017), the median opening angle for the outflow of their low-mass objects is 42 • (their Table 1). We thus chose to study the visibilities in 45 • angular bins separately, namely (i) close to the outflow direction, (ii) close to the equatorial plane, or (iii) close to the ±45 • PA.
The separation of the visibilities with respect to their uv position angle is performed using the GILDAS/MAPPING 3 software. The outflow position angles are provided in Table 1. We plot the radial trends of the dust emissivity located in these different regions in Fig. 5. We did not perform the analysis for L1448-2A and L1448NB1 because they are close binaries driving nonaligned protostellar jets, which confuse the assumption of a simple geometry on the envelope for these sources. We observe that for most of the objects, the trends are similar in all directions. Our data therefore does not point to significantly different dust properties in the outflow direction or the equatorial plane, at least at the 80-2500 au scales we trace in this analysis. While this is surprising at first sight, we stress that such analysis should be repeated on sensitive data probing the dust continuum emission both at the very small and very large envelopes scales, so the outflow cavity walls can be detected and explored individually for example. A better knowledge of the outflow opening angles and orientation would help us refine our crude separation in outflow or equatorial plane regions. We finally note that in L1157, the 3.2 mm emission is more extended than the 1.3 mm emission in the outflow direction , an elongation also observed at 2.7 mm with PdBI (Beltrán et al. 2004) and 3 mm with CARMA (Chiang et al. 2010) 4 . Gueth et al. (1997 suggested that the asymmetry could be linked with an extended component arising from the compressed outflow cavity edges (their Fig. 10). In the visibility plots, we observe that the dust emissivity index seems to be slightly lower in the outflow direction. The trend is also observable in the β 1−3mm map in

Applying a sub-arcsecond component correction
The central sub-arcsecond emission can affect the observed spectral index estimated at larger scales when the analysis is performed in the uv domain, as shown by the work of Miotello et al. (2014) in which unresolved components were shown to affect the emission at all uv distances. Their results are reinforced by the tests performed by our team using predictions from synthetic observations of MHD simulations that show that the α 1−3mm values of the simulation are recovered for most of the envelope once our central region correction is applied (see Appendix B). If the brightness temperature profiles derived for our sources suggest that the Class 0 envelopes studied in this work are mostly optically thin at millimeter wavelengths (see Sect. 2.4 and Fig. 2), the low β 1−3mm values we observe could be due to unresolved optically thick emission in the central regions of some of our sources, but also partly driven by the presence of a central compact component with a different grain population and thus a different emissivity. Moreover, we saw in Sect. 2.3 that on small scales (beyond 200 kλ), the emission could be non-negligibly contaminated by nonthermal dust emission, especially at 3.2 mm, which would also bias our β 1−3mm values. To correct for these contaminations, we subtracted the >200 kλ average amplitude from the shorter visibilities at both wavelengths to recalculate corrected β 1−3mm values. These average amplitudes are tabulated in Table B.1. Because of the asymmetry of the systems at small scales, we again did not perform the analysis for the two binary sources L1448-N and L1448-2A.
We also did not include L1448-C whose 3 mm flux is dominated at most scales studied in this work by the compact component and for which our correction would lead to extremely uncertain β 1−3mm values. To take the uncertainties of the correction into account, the corrected β 1−3mm values at 500 au and gradients are estimated by generating 2000 sets of radial β 1−3mm values varying randomly within their individual error bars following a normal distribution around their nominal value. The corrected β 1−3mm values and gradients are provided in Table 2. The orange dashed lines in Fig. 5 indicate the β 1−3mm corrected radial trends when the central region correction is applied.
Robustness of the correction -We performed a series of tests before choosing our final 200 kλ threshold to ensure that our findings were robust. When we subtracted 150 kλ average amplitude for instance, we find that most of the sources have the same β 1−3mm radial trends. Two sources have a modified slope compared to the 200 kλ correction. SerpS-MM18 has a steeper slope for the 150 kλ correction but the gradient stays within the 40% error bar quoted in Table 2. The second source is L1527 for which a 150 kλ correction leads to a flat beta profile. We note however that the error bar for this source is one of the largest of our sample: the flat β gradient for the 200 kλ correction quoted in Table 2 is also consistent with no gradient at all in this source, within the error bars.

Impact on the β 1-3mm values and radial gradients
After correction, the β 1−3mm values at 500 au range from 0.4 to 1.2. The β 1−3mm values and radial gradients are not affected in the same way for all the sampled sources. Since a significant correction could be the sign of an optically thick compact dust emission, we first explored the link with the presence of a disk. Maury et al. (2019) have detected candidate disks in all but four CALYPSO sources (L1448-2A, L1448-NB1, IRAS 4A1, and SerpM-S68N). The disk is resolved for IRAS 4B, L1448-C, Serp-MM4, and L1527. In the two sources presenting the largest disk-like structures in the dust emission, namely Serp-MM4 and IRAS 4B (radii of 290 and 125 au, respectively), we observe that the sub-arcsecond correction does not significantly affect our previous estimates. This suggests that there is not a significant difference between the dust properties in the envelope and those obtained integrating over the inner 1 region. On the contrary, in L1527 (estimated disk radius of 54 au in Maury et al. 2019), the correction is relatively large, suggesting that the small disk could be partially optically thick or host more emissive (thus lower β 1−3mm ) grains than the envelope. The interpretation is very similar for IRAS 2A1: the presence of an unresolved disk (suggested by Maury et al. 2019) with different grain properties could explain the lower β 1−3mm values derived before the correction is applied. The 1.3 mm and 3.2 mm profiles of L1157 and the sub-arcsecond correction are very similar to those of L1527. Our results even suggest that the β 1−3mm gradient could be reversed once the correction is applied, i.e., where β 1−3mm values decrease toward larger scales, reconciling the gradient estimated in the uv domain with the β 1−3mm map derived in the L1157 envelope (see Appendix C and the previously mentioned 3.2 mm extension observed in the outflow direction). For sources with no disk detected in Maury et al. (2019), namely SerpM-S68N and IRAS 4A1, we find that those are only weakly affected by our central region correction.

Emissivity index dependence on protostar properties
The observed β 1−3mm values stay low in most of our sources despite the central region correction. If the (corrected) β 1−3mm values at 80-2500 au scales mostly stay above 1 in the envelopes of L1157, L1527, and SerpM-S68N, the sources IRAS 4A, IRAS 4B, SVS13, or SerpM-SMM4 have, on average, β 1−3mm values below 1 up to 2000 au scales. No correlation was found between the average β 1−3mm value and age in a sample of Class I and II objects studied by Ricci et al. (2010). We investigate other correlations in this work. In Fig. 6, we compare the β 1−3mm value at 500 au with the bolometric temperatures T bol for each object: the correlation is weak (Pearson correlation coefficient R = 0.5) between the two parameters. We note that this result is consistent with the absence of correlation found between the two parameters in the analysis of Froebrich (2005). Like these authors, we also find an anticorrelation between β 1−3mm and M env (R = −0.79), which partly suggests smaller grains on average in protostars with smaller envelope mass (Fig. 6). We note that M env estimates strongly depend on the assumed distance and are impacted by the difficulty to separate the envelope from the parent cloud. The M env estimates are also usually calculated assuming a fixed dust opacity; for instance with κ 1.3 , the dust opacity at 1.3 mm is fixed to a commonly used value of 0.01 cm 2 g −1 . This could contribute to the anticorrelation we observed between M env (or N H2 ) and β 500au .

Possible steeper emissivity radial evolution in more massive envelopes
As far as the β 1−3mm variations with envelope radius are concerned, IRAS 2A1, IRAS 4A1, SerpM-SMM4, SerpM-S68N, and SerpS-MM18 present a significant evolution of β 1−3mm in their envelopes. All these sources have a relatively large envelope mass (>6 M ). We note in particular that the two Serpens sources SerpS-MM18 and SerpM-SMM4 show the same β 1−3mm gradient. The sources IRAS 4B, SVS13B, L1157, and L1527 have more moderate gradients (<0.5). Those sources have lower envelope masses than the previous group. In Fig. 7, we show that the steepness of the gradient linking β 1−3mm to the envelope scales probed by our observations seems to be strongly correlated to the envelope mass of the object (R = 0.93) by over an order of magnitude. The relation suggests more homogeneous dust properties in the lowest mass sources of our sample. We note that the more evolved object L1527 presents an homogeneous β 1−3mm in its envelope. One explanation could be that the grain populations have been completely mixed over longer timescales. Investigating more potential relations with environmental properties, we do not find correlations between the β 1−3mm gradients with the internal luminosity (Fig. 7); neither do we find steeper β 1−3mm gradients in objects presenting larger snowlines nor a correlation between β 1−3mm values and the position of the CO snowline itself (Fig. 5; Anderl et al. 2016;Gaudel et al. 2019).
Complex organic molecules (COMs) are supposedly formed on the surfaces of dust grains and are released in the gas phase when the dust icy mantles sublimate (Rawlings et al. 2013). Belloche et al (in prep.) performed a COM analysis for the CALYPSO sample, studying the abundance of molecules such as CH 3 OH, C 2 H 5 OH, CH 3 OCH 3 , CH 3 OCHO, CH 3 CHO, NH 2 CHO, CH 3 CN, C 2 H 5 CN, a-(CH 2 OH) 2 , CH 2 (OH)CHO, HNCO or NH 2 CN. IRAS 2A1, IRAS 4B, and SerpS-MM18  Table 2. Colored/empty symbols represent the values after/before the central region correction (c.r.c).
have clear detections of most of the cited COMs but do not particularly share a common trend in terms of β 1−3mm values or gradients. Both SerpM-S68N and L1157 only have CH 3 OH, CH 3 OCHO, and CH 3 CN detected but the former has the steepest β 1−3mm gradient while the latter the lowest. The absence of correlation between COM abundances and the radial variations of the dust emissivity index we observe thus suggests that the decrease of the emissivity index at small envelope radii is not driven by the loss of the icy mantles.
Finally, we plot in Fig. 8 the variation of β 1−3mm as a function of the dust temperature at the scale it is measured (Eq. (6)), allowing a source-to-source local comparison. The large range of β 1−3mm values obtained at similar dust temperatures across the sample envelopes suggests that the local dust temperature does not have a major influence on the dust emissivity itself.

Millimeter-size grains: grown in situ or transported
Several sources have β 1−3mm values below 1. Miyake & Nakagawa (1993), Ossenkopf & Henning (1994), or Draine (2006 have shown that grain growth only could already produce A5, page 10 of 17 β 1−3mm ≤ 1 at millimeter wavelengths. A recent study by Ysard et al. (2019) has quantified the size distribution effects on the shape (peak, β) of the spectral energy distribution (SED) of dense interstellar regions including young stellar objects. Their results (see their Fig. 10) clearly show that the millimeter spectral indices decrease with an increasing size distribution. For most of the material mixtures probed, β values <1 require grain sizes larger than 100 µm. Many questions remain regarding how these large grains can form in the youngest protostellar envelopes despite the long timescales and high densities required to build them. Ormel et al. (2009) showed that long support mechanisms are needed for coagulation to happen and impact the β 1−3mm observed in molecular clouds. These authors found that millimeter-size particles require ice-coated grains and coagulation times of 10 7 yr to be formed from ISM-like grains. This timescale shrinks when considering hydrogen number densities >10 6 cm −3 (densities that are reached in our protostellar envelopes) or if the coagulation processes start on already grown grains, as expected to explain the coreshine effect observed in cold cores. The denser the medium, the stickier the grains: this enhances the coagulation processes in the inner envelopes during the protostellar collapse.
More recently however, Wong et al. (2016) proposed new coagulation models and found that higher densities (n H = 10 10 cm −3 ) would be required to grow grains above 300 µm in reasonable (a few free-fall time) timescales. To explain the presence of millimeter-size grains in the Class 0 envelopes, these authors suggested that large grains could be built up in the highdensity central region (in the protostellar disk), then could be "launched" into the inner envelope via the gas dragged by the outflow. The transportation conditions would depend on many criteria: for instance, the central mass of the stellar embryo, i.e., the larger the mass, the less efficient the transportation; the gas velocity, i.e., low v gas is needed to transport millimeter grains; the opening angle of the outflow; its mass loss rate; and the sputtering and shattering rate, namely, the ability of a grain to survive gas friction and grain-grain collisions.
Since the outflow momentum flux correlates well with the circumstellar envelope mass (Bontemps et al. 1996), the correlation we observe between M env and the β 1−3mm gradient could be linked with stronger outflow energetics in high-mass protostellar envelopes, leading to a more efficient transport of large grains from the inner disk to the inner/intermediate envelope scales. Interestingly, the detection of strongly polarized (from a few percent to up 10%) dust emission at millimeter wavelengths following precisely protostellar outflow cavity walls (Maury et al. 2018;LeGouellec et al. 2019) suggests a population of large grains (Valdivia et al. 2019) in these specific locations, which strengthens this scenario. The scales on which dust grains could be reinjected into the envelopes using this mechanism and outflow dynamics (gas velocity and grain motions) need to be constrained observationally, especially at the small scales at which the outflow is launched and the cavity develops. Current efforts in the theoretical aspects of these investigations, in particular the ongoing implementation of the dust grain dynamics into dusty collapse simulations (e.g., Lebreuilly et al. 2019) along with further tests on the decoupling of the dust and gas and potential grains segregation sorting in the equatorial plane (Bate & Lorén-Aguilar 2017), might provide valuable predictions to carry out additional observational tests to constrain the full dust evolution from the core to disk scales.

Grain composition effects
The grain chemical composition is expected to vary during the collapse of protostellar envelopes as a result of, for example, grain-grain collisions and thermal desorption. Composition effects could thus also participate in the low β values we observed at millimeter wavelengths. So far, the very low values of β 1−3mm ( 1) have been difficult to explain from optical properties derived in laboratory studies of interstellar dust analogs. Köhler et al. (2015) showed that additional carbonaceous mantles or ices would not lead to such low beta. Studies on magnesiumrich or iron-rich amorphous silicates (Demyk et al. 2017a,b) have shown that the mass absorption coefficient strongly varies with wavelength but the derived millimeter β values were still higher than 1. Recent work using the optical constants from The Heterogeneous dust Evolution Model for Interstellar Solids (THEMIS; Jones et al. 2013) dust model has highlighted the strong influence of the chemical composition on the emission of various mixtures of amorphous hydrogenated carbon grains and amorphous silicates (Ysard et al. 2019). Their study, however, shows that large grains seem to be a necessary condition to explain β mm < 1. Coupeaud et al. (2011) suggested that astrosilicate grains with a stoichiometry close to that of pyroxene could have β < 1 in the millimeter range, at which pyroxenes are more abundant, for a given pressure, at lower temperatures (Kessler-Silacci et al. 2005). Several studies of massive protostars and protostars from the Orion molecular cloud complex seem to confirm that their amorphous silicate population is more dominated by amorphous pyroxene than the typical ISM dust (Demyk et al. 1999;Poteet 2012). The object-to-object variations of β 1−3mm could finally be partly linked with variations in the grain porosity itself (see for instance Kataoka et al. 2014). In conclusions, while a varying grain size distribution remains our favored hypothesis to explain the radial trends we observe, a larger spectrum of laboratory experiments would help predict the observational signatures expected for a variety of interstellar dust analogs and from the collisions of these grains. More observations of the dust polarized emission (and their use as a constraint for polarized dust models) could also provide important additional information on the dust composition (Guillet et al. 2018).

Other potential caveats
Studies of 3 mm dust continuum measurements in nearby filaments (Schnee et al. 2014;Mason et al. 2019) recently suggested that the emission at 3 mm deviates from dust emission extrapolated from a simple modified blackbody SED fitted to shorter wavelengths. The emission is inconsistent with free-free emission or spinning dust models and the enhanced emission at 3 mm has been attributed to the presence of amorphous dust grains. The question regarding whether 3 mm dust continuum emission can be safely used to study grain size variations in protostellar environments remains, however, open. Observations at additional (sub)mm wavelengths but also at higher resolution would enable us to sample better the local submillimeter SEDs to disentangle different grain models. Future James Webb Space Telescope (JWST) instruments will also provide the opportunity to observe gases, silicate features (e.g., analyzing the 10-18 µm ratio as done in Demyk et al. 1999), or ices in absorption, probing more deeply into the warm interior of the protostellar envelopes. Galván-Madrid et al. (2018) have recently explored the effects of self-obscuration in protostellar disks with a radially decreasing temperature gradient in producing extremely A5, page 11 of 17 A&A 632, A5 (2019) low spectral indices. Given that most of the emission arising from the Class 0 protostellar envelopes studied in this work is optically thin, our results are probably not affected by these self-obscuration effects.
Finally, our team is currently producing a series of synthetic observations to quantify the impact of various parameters (inclination, size distribution, luminosity, and envelope mass) on the observed millimeter continuum and polarized emission (Valdivia et al., in prep.). This work has also been initiated to investigate the instrumental effects linked with interferometric observations and get a better handle on whether they could partly bias the β 1−3mm trends we observe.

Conclusions
We study the dust properties in a sample of 12 Class 0 protostellar envelopes observed at 94 and 231 GHz with PdBI as part of the CALYPSO survey. We perform our analysis in the visibility domain. We derive dust emissivity index profiles as a function of the envelope scales probed. Since the central arcsecond emission can affect the observed spectral index estimated at larger scales, we correct the visibility profiles for this potential contamination. Our analysis leads to the following conclusions: (i) Most of the Class 0 envelopes present low dust emissivity indices β 1−3mm : the observed β 1−3mm often reach values lower than 1, in particular in the envelopes of IRAS 4A, IRAS 4B, SVS13B, and SerpM-SMM4. Many sources also present a strong decrease of β 1−3mm toward small envelope radii. Low values of β 1−3mm are found preferentially in sources presenting these large radial variations. Both results suggest that grain growth could already be at work in Class 0 inner envelopes.
(ii) The steepness of the β 1−3mm gradient depends strongly on the envelope mass, suggesting less grain evolution in low-mass objects and a more diverse grain mixture from the outer to the inner regions in more massive envelopes.
(iii) If accretion and coagulation models predict the presence of large grains in protostellar envelopes, the timescales required seem to be too large for millimeter-size grains to be grown directly in the envelopes from ISM-like grains. Potential transportation mechanisms of large grains formed in the high-density inner envelopes and protostellar disks back to the intermediate scales of the envelopes might explain the observations and should be observationally investigated in the future.
(iv) Variations in the grain composition, stoichiometry, or porosity also strongly influence the dust emission at millimeter wavelengths and could partially contribute to the radial trends observed at envelope scales. Recent predictions from a range of interstellar grain mixture analogs suggest however that a certain level of grain growth would still be necessary to explain β 1−3mm values <1 (Ysard et al. 2019).
Higher resolution observations, with JWST and ALMA, will help probe deeper into the core and connect the dust properties down to the disk/envelope transition to analyze locally if observations fit with the different accretion, coagulation, and transport model predictions. A&A 632, A5 (2019) Appendix A: Free-free and synchrotron contributions to the 1.3 and 3.2 mm emission Notes. (a) As in Tobin et al. (2015), the spectral slopes α ff+s are defined by F ν = F 0 (λ/λ 0 ) −α ff+s , with λ 0 = 1 mm. In order to quantify the contribution to the PdBI observations not coming from thermal dust emission, we assume that the ≥2 cm emission can be considered as free from thermal dust emission (as observed in IRAS 2A by Tobin et al. 2015) and considered as compact. We gathered fluxes beyond 2 cm published in the literature. For sources with more than two observational constraints, we estimated the free-free + synchrotron spectral slope α ff+s as where λ 0 = 1 mm. We then extrapolated the fitted spectrum to the millimeter regime and compare it with the 1.3 and 3.2 mm PdBI fluxes at 200 kλ (so on our smallest scales) to estimate the contamination at both wavelengths. The references of the used radio fluxes, the derived α ff+s slopes, and the contamination to the 1.3 mm and 3.2 mm emission are reported in Table A.1. For SerpM-SMM4, where only the 3.6 cm flux was found, we fixed α ff+s to 1.1, which is the α ff+s value of IRAS 4A1, i.e., the source in the sample that is the closest in terms of envelope mass, temperature, and internal luminosity. For SerpS-MM18, the α ff+s value was derived by Kern et al. (2016) from observations at 4.75 and 7.25 GHz. We note that several millimeter sources reside in the 4 beam, as observed in Plunkett et al. (2018). The high α ff+s found by Kern et al. (2016) (∼2-2.3) suggests that thermal dust emission might still contribute to the centimeter emission. As in SerpM-SMM4, we decided to fix α ff+s to 1.1 for this source. The contribution is systematically lower than 10% at 1.3 mm, but it can become non-negligible at 3.2 mm. Calculated using the radio sources of Kern et al. (2016) for SerpS-MM18, the contamination accounts for all the 200 kλ 1.3 mm and 3.2mm fluxes; these fluxes have a large α ff+s (>2) usually associated with optically thick emission. The observations have shown that the analysis of the spectral index performed in the uv domain can be biased by unresolved optically thick emission or the presence of central compact components. We subtracted the >200 kλ average amplitude to correct for these contaminations. These average amplitudes are tabulated in Table B.1. In order to further test the effects of this correction on the 1.3 and 3.2 mm visibility profiles and deduced β 1−3mm , we generated synthetic observations starting from a nonideal MHD simulation of the gravitational collapse of a prestellar core. The simulation is generated using the RAMSES code (Teyssier 2002;Fromang et al. 2006). It was carried out by P. Hennebelle in a fashion similar to the non-ideal MHD models of protostar formation shown in Masson et al. (2016). The dust radiative transfer output from this MHD model that we use here was already presented in Valdivia et al. (2019). The core has a mass of 2.5 M , a mass-to-flux ratio µ = 6, and a rotational to gravitational energy ratio of 10 −2 . The effective resolution is  Maury et al. 2018). We selected a snapshot in which the central embryo has accreted 0.2 M and we post-processed a zoomed-in region of 8000 au using the radiative transfer code POLArized RadIation Simulator (Reissl et al. 2016, POLARIS). We assumed a gas-to-dust mass ratio of 100 and a MRN-like (Mathis et al. 1977) dust grain size distribution (n(a) ∝ a −3.5 ; 5 nm ≤ a ≤ 1 µm). The composition is 62.5% astronomical silicates and 37.5% carbonaceous grains: the simulation includes neither dust evolution nor dust dynamics. The chosen accretion luminosity and distance are 1 L and 120 pc, respectively. Figure B.1 shows the integrated dust emission of the protostellar envelope at 1.3 and 3.2 mm, with a density distribution is close to a simple power law, a pseudo-disk, and bipolar outflow bubbles extending to up to 800 au. We show how the predicted spectral index (Eq. (2)) varies locally in the synthetic envelope in the bottom panel, where α ∼ 3 in the outer part of the envelopes, increasing up to 3.3 in the outflow bubble, then decreasing drastically below 2 in the central 50 au where an optically thick disk component resides. The 1.3 and 3.2 mm emission maps are used as inputs for the ALMA simulator of the CASA 5 software. We ran the simobserve task to produce four simulated measurement sets for compact and extended ALMA antenna configurations, C43-1, C43-3, C43-5, and C43-8 6 . We used integration times of 7200 s per configuration and assumed a PWV of 0.7 mm (i.e., under good atmospheric conditions). Thermal noise (default ATM model) was added to the simulated data. We then extracted the concatenated visibilities using GILDAS/MAPPING. Figure B.2 shows the predicted observations in the uv domain at 1.3 and 3.2 mm and their associated spectral indices α. The α values plunging at small scales can be explained by the central (∼50 au) region that, with its weak (α < 3) values, affects the whole visibility profile. To correct for this contribution, we masked the inner region: we calculated the fluxes at 1.3 m and 3.2 mm at 200 kλ (to be sure to encompass the 50 au central region) and removed them from the visibility amplitudes at all shortest baselines. The corrected trend for α is overlaid in Fig. B.2 (dashed blue line). The α ≥ 3 values observed in the α map (Fig. B.1 bottom) are recovered for most of the envelope; there is a nearly constant α after the correction. We also reran the tests using a simulation in which the temperature is fixed to 10 K during the post-processing, but this does not affect our results: we conclude that temperature gradients along the line of sight should only have a minor effect on our results.

Appendix C: Maps of β 1-3mm
For a qualitative assessment of the spatial variations of the dust emissivity index β 1−3mm , we performed an imaging of the dust continuum emission (including the secondary sources) using GILDAS/MAPPING. We applied a robust weighting scheme, compromise between the natural weighting and the uniform weighting that minimizes the synthesized beam sides lobes in particular. Table C.1 provides the synthesized beam and the rms of the obtained 94 and 231 GHz maps. In order to build maps of β 1−3mm , we need to have both 94 and 231 GHz maps tracing the same physical scales at the same resolution. We thus apply an additional uv tapering on the visibilities at 231 GHz to obtain beams at 231 GHz similar (in size and orientation) to those at 94 GHz. Maps are produced with the same pixel size. The "tapered" sizes and rms of these new maps are also given in Table C.1 (231 GHz-tapered columns). We finally used the temperature profiles estimated for each source to apply Eq. (7) on a pixel-by-pixel basis. Figure C.1 shows the β 1−3mm maps of the sample. β 1−3mm is shown at locations where the emission at 1.3 and 3.2 mm is detected with a signal-to-noise ratio (S/N) higher than two. We observe that the dust emissivity index changes in the envelopes of most of the sources as observed in the uv domain. However, when β 1−3mm values are averaged over various aperture sizes, we observe that β 1−3mm gradients are not as clearly observed as we probe smaller and smaller scales compared to the results obtained from the uv domain analysis. This demonstrates the difficulty of analyzing β 1−3mm spatial variations directly from interferometric maps affected by filtering. Only pixels with a 2σ detection at both wavelengths are shown. The contours indicate the 94 GHz emission detected at 3, 5, 10, 30, and 80σ. The blue and red arrows indicate the direction of the bipolar jets. The insets show how β 1−3mm evolves as values get integrated to larger and larger scales. The x-axis of the inset plots is labeled in astronomical units.