Revealing the dust grain size in the inner envelope of the Class I protostar Per-emb-50

A good constraint of when the growth of dust grains from sub-micrometer to millimeter sizes occurs, is crucial for planet formation models. This provides the first step towards the production of pebbles and planetesimals in protoplanetary disks. Currently, it is well established that Class II objects have large dust grains. However, it is not clear when in the star formation process this grain growth occurs. We use multi-wavelength millimeter observations of a Class I protostar to obtain the spectral index of the observed flux densities $\alpha_\mathrm{mm}$ of the unresolved disk and the surrounding envelope. Our goal is to compare our observational results with visibility modeling at both wavelengths simultaneously. We present data from NOEMA at 2.7 mm and SMA at 1.3 mm of the Class I protostar, Per-emb-50. We model the dust emission with a variety of parametric and radiative transfer models to deduce the grain size from the observed emission spectral index. We find a spectral index in the envelope of Per-emb-50 of $\alpha_{\rm env}$=$3.3\pm0.3$, similar to the typical ISM values. The radiative transfer modeling of the source confirms this value of $\alpha_{\rm env}$ with the presence of dust with a $a_\mathrm{max}$$\leq$100 $\mu$m. Additionally, we explore the backwarming effect, where we find that the envelope structure affects the millimeter emission of the disk. Our results reveal grains with a maximum size no larger than $100$ $\mu$m in the inner envelope of the Class I protostar Per-emb-50, providing an interesting case to test the universality of millimeter grain growth expected in these sources.


Introduction
Disks and envelopes around protostars play a fundamental role in the process of planet formation since they contain the ingredients out of which planets are formed . Thanks to detailed studies of protoplanetary disks at several submm and mm wavelengths such as HL Tau (Carrasco-González et al. 2016), CY Tau, DoAr 25, and FT Tau (Pérez et al. 2015;Tazzari et al. 2016) , it is now well established that the radial profiles of their grain size distributions are compatible with millimeter size grains. However, it is not yet clear at which stage of the star and planet formation process dust grains start to efficiently coagulate and evolve from µm size particles to macroscopic dimensions. Ormel et al. (2009) studied in detail the possibility of grain growth in pre-stellar cores and found that while it is easy to grow to micron size particles, the growth to millimeter or centimeter size pebbles requires high densities and relatively long Based on observations carried out under project number S16AT with the IRAM NOEMA Interferometer. IRAM is supported by IN-SU/CNRS (France), MPG (Germany) and IGN (Spain). timescales of ∼10 7 yr, much longer than the lifetimes of dense cores. This is also explored recently in Chacón-Tanarro et al. (2017), where they calculate the grain size in the center of the pre-stellar core L1544, finding that only in the central 300 AU, grain size can grow to about 200µm. In the earliest protostellar phases, e.g. during the Class 0 stage, the protostar is fully embedded in the parent envelope, while in the Class I phase, the envelope is partially dissipated and the disk emission can be better separated from the envelope. Therefore, Class I protostars can more easily address the start of planetesimal formation and constrain the initial conditions of the evolution of protoplanetary disks. The possibility for the first large solids to assemble during the early phases of disk evolution would have important implications. If the process starts already in the Class I stage it would imply a much more effective and rapid planetesimal formation phase in the disks. In fact, if large (mm to cm-size) dust particles from the inner envelope (Chiang et al. 2012a;Tobin et al. 2013;Miotello et al. 2014) are deposited in the disk at large radii during the disk formation stage, they would be much less affected by the radial transport and fragmentation processes, which ad-

The source
Per-emb-50 is a protostar located in the active cluster forming region NGC1333 in the Perseus cloud (see Figure 1), at a recently revised distance of 293 pc (Ortiz-León et al. 2018;Zucker et al. 2018). It is classified as a Class I protostar from the slope of its SED in the near-, mid-infrared ("Cores to Disks" or c2d Spitzer Legacy project from, Evans et al. 2003). Based on Bolocam 1.1mm data, the bolometric temperature is T bol =254±23 K. The rescaled bolometric luminosity is L bol =13.7±3 L , making it one of the brightest Class I sources in Perseus. High angular resolution observations conducted at 8mm in the VLA Nascent Disk and Multiplicity (VANDAM) survey, provide a lower limit for the disk mass and outer disk radius. The rescaled values from Segura-Cox et al. (2016) for mass and radius are: M disk =0.28-0.58 M and r out =27-32 AU, respectively. Literature values for envelope mass, disk mass, disk inclination, and other parameters are presented in Table 1. We note that some of these physical parameters were calculated using the 230 pc from Hirota et al. (2008) or 250 pc in the case of Bolocam observations, therefore, we rescale the limits taking into account the different distance adopted. Even though Per-emb-50 presents a small disk at 8 mm, it is the perfect candidate for studying the growth in the inner envelope and their dust properties. References.

SMA observations
The Submillimeter Array (SMA) data shown in this paper are from the MASSES legacy program (Mass Assembly of Stellar Systems and their Evolution with the SMA, PI: I.W. Stephens, M. Dunham; e.g., Stephens et al. (2018)). Per-emb-50 was observed at 1.3 mm with the receiver centered at 220.69 GHz, in the Extended (eight antennas) and Subcompact (seven antennas) configuration with ASIC (Application Specific Integrated Circuit) correlator during September 2015 and November 2014, respectively. Additionally, Per-emb-50 was observed during October 2015 with SWARM (SMA Wideband Astronomical ROACH2 Machine) correlator at 1.3mm in extended configuration (see Table 2 for more details). Weather conditions were good, with zenith optical depths at 220 GHz of τ 220 = 0.07 -0.15. Calibration was done in MIR while imaging was done in MIRIAD (Sault et al. 1995), using the standard calibration procedure. We inspected the amplitudes and phases of the calibrators on each baseline in order to look for variations or noisy data, which were manually flagged. Corrections for system temperatures were applied in order to calibrate the atmosphere attenuation in the visibility amplitudes. Detailed information of the calibration can be found in Stephens et al. (8C) in operation during the first track, and 6 antennas (6C) in operation for the second. Antennas were based on stations E10, W20, W10, N20*, E18, N11* 2 , N17 and E04. The projected baselines were from 7.8 kλ and 102 kλ. Per-emb-50 was observed for hour angles from -5.8 to 1.5 h for 8C, and from -5.3 to 1.4 h for 6C. In total we spent 9 hours on source. 0333+321 was used as phase/amplitude calibrator. The sources LkHa101 and MWC349 were used for the flux calibration, while the quasars 3C84 and 3C454.3 were used for the bandpass calibration. We consider an absolute flux uncertainty of 10%. The total bandpass for the 110 GHz continuum measurement was 2 GHz. Data reduction and image synthesis were carried out using the GILDAS software (Guilloteau & Lucas 2000) with the procedure of MAPPING> Selfcal. The continuum map ( Fig. 2) was produced using natural weighting and the resulting beam size is 2. 1 × 1. 6 at P.A. 34.84 • . The clean map has a rms noise level of 2.1 mJy beam −1 .

Observational Analysis
Since we are working with interferometric data, the best way to analyze our source is working on the visibility domain. This is to avoid biases in the model-data comparison that are introduced by the CLEAN algorithm, u-v sampling, and the imaging process.
In Fig.2 we plot the real visibility as a function of the deprojected baseline length (uv-distance). The deprojected uvdistances are given by R= Lay et al. 1997). The values for inclination i, and position angle, PA, are presented in Table1. In Fig.3, we show images of the SMA and NOEMA observations. Per-emb-50 appears as a point source in these two images, so we do not resolve the embedded disk. Consequently, since the disk is unresolved, then it contributes as a constant component at all baselines. At long baselines, we expect that the amplitude of the visibility is dominated by the disk component, while in shorter baselines the resolved envelope dominates. For Per-emb-50, the value of the amplitudes start becoming constant above 47 kλ (see Fig.2  Using the fluxes between the u-v ranges 47-80 kλ at 1.3 and 2.7 mm, we obtain therefore the average value α mm in the unresolved disk, which is α disk =1.71±0.3. As shown in the Fig.2, an excess of emission is present at short baselines (<47 kλ) at 2.7 mm and 1.3 mm, which correspond to physical scales of 1500-3000 AU. The excess values at these baselines, after subtracting the disk visibilities, are:  =10.1 ± 5 mJy. The excess, even if not very pronounced at 2.7 mm, is detected at 1.3 mm, therefore this indicates the presence of extended emission related to the inner envelope. Considering the average fluxes at these very short baselines, and using the same u-v distances ranges at both wavelengths, we recover an average value for α env bigger than the typical ISM values, α env =4.0 ± 0.8. The uncertainty in α disk and α env are estimated following the procedure shown in Appendix A, with the absolute flux uncertainty of 10% for 2.7 mm data and 20% for 1.3 mm data added in quadrature. We present α env and its change as a function of deprojected baseline in Fig.4. If we translate this value to the spectral index of the dust opacity, we obtain a β env ∼2.0, which is similar to the values found in the ISM.
These preliminary results are showing a discrepancy with previous studies on spectral indexes in Class I protostar or even younger sources with dust opacity indexes α mm <3 (Miotello et al. 2014;Kwon et al. 2009;Chiang et al. 2012a). To investigate possible explanations, we will perform a partial and full radiative transfer modeling on envelope and disk to take into account possible deviations from the optically thin and Rayleigh Jeans regimes, which can affect the values of α mm .

Modeling
In order to model the Class I protostar and compare with the observations, we consider appropriate physical structure and conditions of the source, including the envelope structure, density, properties of the dust grains, and we predict the 1.3 and 2.7 mm emission with a u-v modeling described below. In the first step, we fit the Per-emb-50 data with a parametric modeling in uv space (Section 4.1) in order to address the α values, as well as visibility comparisons. Afterward, we use the radiative transfer tool RADMC-3D (Dullemond et al. 2012) in two ways: (a) to apply the modeling approach of Miotello et al. (2014), where the disk and envelope are modeled separately (Section 4.2), (b) to compute the emission for the new modeling presented in this work, that include a self-consistent radiative transfer model for the disk and the envelope (Section 4.3). In the following sections, we discuss and compare the details of the results of each modeling case.

Parametric Model
We implement a model that consists of an extended envelope described by a Gaussian, and an unresolved disk (point source) that has constant flux at all baselines. Therefore, the combined amplitude profile, that depends on the uv distance, defined as √ u 2 + v 2 , and frequency, ν, is described by: where F e and F d are the flux density from the Gaussian emission (extended envelope) and point source emission (unresolved disk) respectively, α e and α d are the spectral indexes of the two components, and σ is the width of the Gaussian given by σ ≈ FWHM/2.355. In this simple model, we first set the flux from the disk at 1.3 mm based on the average value reported in Sec. 3, F d =63 mJy. Then, four parameters are explored: F e , α e , α d and σ. The Markov chain Monte Carlo (MCMC) method, implemented as a python package emcee (Foreman-Mackey et al. 2013), is utilized to calculate the posterior probability distributions of each of these parameters. For each model we used the 750 steps after the burn-in and 400 walkers (see Appendix B for more details). The results from this simple model will be discussed in the following section.

Parametric model results
In Fig. 5, model visibilities are compared with observational data at each u−v sample and wavelength. In Table 3 we present the best-fit parameters found for this parametric model. The values of the flux spectral index in the disk and envelope are consistent with the observational analysis (see Section 3), but their errors are highly dominated by the systematic error of absolute fluxes and the statistical error of the data. We estimate that the uncertainty on α mm for the envelope and disk using a simplistic approximation for non correlated errors is ±0.3. Additionally, from this simple model we can constrain the size of the region where the envelope emission arises, which the 1sigma width (from Table 3) 43 kλ (1405 AU). From the model we can derive the flux from the disk at 2.7 mm and the prediction of the total flux at baseline=0 kλ or zero spacing flux, F 2.7mm zero and F 1.3mm zero . The results are shown in Table 4. We use the derived parameters from our parametric model (Table 4) to estimate the zero spacing flux at 1.1 mm, F 1.1mm zero . We find the flux is only 127.4 mJy beam −1 which is much lower than the single dish flux of 612±18 mJy reported by Enoch et al. (2006). This discrepancy is related to the resolution of the observations. While our interferometric data is sensitive to the inner envelope of this source, the 31 beam size of Bolocam is recovering the extended emission, which is affected by blending effects, especially in a crowded region such as NGC 1333 (see Fig. 1). Since this simple model is not taking into account properties of the dust grains and density profiles for both the envelope and disk, we also analyze Per-emb-50 with more detailed dust radiative transfer models.
Article number, page 5 of 17 A&A proofs: manuscript no. main For this model, we adopted the procedure described by Miotello et al. (2014), where they analyzed two Class I protostars with a 2−step model. The disk is modeled adopting the two-layer model by Dullemond et al. (2001), whose output spectrum is taken as central source of illumination in the envelope model. The envelope, on the other hand, is modeled using RADMC-3D (Dullemond et al. 2012).

Modeling Protostar and Disk
We adopt a simple disk model heated by protostellar radiation.
We calculate the properties of the central protostar, assuming that it emits black body radiation, characterized by a radius R , effective temperature T eff , and mass M . To obtain T eff we assume that Per-emb-50 lies along the birthline for intermediate mass stars by Palla & Stahler (1990). Given the rescaled bolometric luminosity L bol reported in Table 1, we estimate T eff =5011 K. With L bol and T eff , we can estimate R eff using the Stefan-Boltzmann law: Then, with R eff =5.01 R , we use the mass vs. radius relation for a spherical protostar accreting at a rate of 10 −5 M yr −1 from Palla & Stahler (1991), to deduce an effective mass of M eff =2.9 M (Table. 5). Additionally, we add a disk structure defined by an inner and outer radius, r in and r out , an inclination angle i, and a dust surface density profile that follows a simple power law, where Σ 0 is the surface dust density fixed at r Σ 0 = 1 AU from the central protostar, and where p=1 since the quality of the data is not sensitive enough to discriminate between different values of p. The disk inclination i is fixed to 67 • as found by Segura-Cox et al. (2016). Since the mm-SED is not sensitive to r in , we set r in =0.1 AU. r out and M disk can be constrained by our observations assuming a dust opacity (see Section 4.2.3) and gas-to-dust mass ratio of 100.

Modeling the Envelope
We adopted the rotating and collapsing spheroid structure by Ulrich (1976) to model the envelope. The density of this envelope structure is given by, where ρ 0 is the density in the equatorial plane at the centrifugal radius R rot of the envelope, and θ 0 is the solution of the parabolic motion of an infalling particle given by: The outer radius of the envelope is fixed at 8 800 AU, which is equivalent to the 30 aperture of Bolocam. In this case we will use the envelope mass derived by Bolocam to compare with the models. We computed ρ 0 by imposing a total envelope mass M env , and R rot , which can have a significant influence on the amplitude as a function of baseline, it was left free to vary. Outflow cavities are not included in this model. RADMC-3D is used to compute the temperature of the envelope, with the implementation of Eq. (5) to describe the density structure. The protostar and disk system presented in the previous subsections are used as heating source of the envelope, whose emission is calculated using the two-layer model by Dullemond et al. (2001), and then the output spectrum is used in the 2D radiative transfer calculation for the envelope structure.

Dust opacity
We adopt the dust opacity model used in Ricci et al. (2010). A dust population characterized by a distribution of grains with different sizes was implemented. We used a truncated power law distribution n(a) ∝ a −q , between a minimum and a maximum grain size, a min and a max respectively. We fixed the chemical composition to a silicate, carbonaceous material and water ice in a 1:2:3 volume fractional ratio. Additionally, we set a min =0.01 µm and we use q=3.0. We varied a disk max and a env max according to the range presented in Table 5.

Model fitting
To compare the model with the interferometric observations, we have to create images at the exact wavelengths of our observations. Then, those model images have to be transformed to model visibilities. For that we used the computational library GALARIO (Tazzari et al. 2018). The model image is convolved with the primary beam patterns of the antennas and then Fourier transformed into visibilities. The first step in this modeling is to fit the disk emission. We created a grid of parameters varying M disk , R out and a max to reproduce together F 1.3mm disk and F 2.7mm disk . Once we found the three parameters that match F 1.3mm disk =63.97 mJy and F 2.7mm d =18.8 mJy, we implement these output fluxes (output spectrum) as the heating central source of the envelope. Then, using RADMC-3D (Dullemond et al. 2012), we vary M env and a env max in order to reproduce the interferometric fluxes at 1.3 mm and 2.7 mm. Table 5 gives a complete list of models parameters and indicates whether they are fixed or varied. In Fig. 6 we present the best fit for the observed visibilities at both wavelengths. The set of parameters that provided the best match with the observations are presented in Table 6. The two best fit are discussed in the next section.

Results Two step model
The parameters that provide a good fit respect to the disk emission at both wavelengths are reported in Table 6. The model Article number, page 6 of 17 C. Agurto-Gangas et al.: Revealing the dust grain size in the inner envelope of the Class I protostar Per-emb-50  For the envelope, we explore the effects of changing: R rot , a max and ρ 0 . The envelope inner radius is fixed at the outer radius of the disk model. We tested different R rot between 100-1000 AU to accommodate the total enclosed envelope mass. As mentioned in Crapsi et al. 2008, decreasing the centrifugal radius results in more peaked and spherical envelopes. Using a small centrifugal radius has a significant influence on the amplitude at short baseline length. For example varying the centrifugal radius by a factor of 2 changes the first amplitude point of the model by 20%. We found that a R rot of 600 AU is consistent with the slope at short baselines in both wavelengths. We can constrain the level of grain sizes in the envelope within the framework of the collapsing rotating envelope model. For example, in Fig. 6, if we consider a dust grain size distribution in the envelope with a maximum size of 1 mm, we can reproduce the 1.3 mm observations, but we underestimate the total envelope mass by a factor of 6. In the case of models with 0.1 µm< a max < 100 µm, the flux at 1.3 mm and 2.7 mm matches the observations very well, but the derived envelope masses differ from those derived from observations. The best match with the 2.2M envelope mass derived by Enoch et al. (2009) are those derived from models with dust grain sizes of a max ≤ 50µm (see Table 6). The models with a max = 100 µm recover almost 60% of the envelope mass. Table 7 presents the derived masses for the envelope using different a max in M1 and M2. A distribution of grains with a max ≤50 µm provides a good match with the observations since the flat emission at 2.7 mm matches the observations well and is consistent with the systematic errors due the flux calibration. Based on this model, the maximum grain sizes in the envelope   Table 6). In solid lines we present models with grain sizes of a max ≤ 100 µm. In dashed lines are models with grain sizes of a max =300,1000 µm. The best fits are the models with a distribution of grain sizes with a max ≤ 100 µm. The red shaded region is the uncertainty on the data due to flux calibration. The bottom of each panel shows the residuals between the data and the model with different a max . are unlikely to be larger than a hundred microns. This would imply that the envelope may have gone through a process of grain growth, but there is no evidence that a substantial fraction of grains are large millimeter-sized dust aggregates. As we mention before, the observed flux and spectral index of Per-emb-50 are consistent with a small optically thick disk, in which case, we cannot constrain the spectral index α.
For the envelope we can use our dust model to infer the value of β, which is β env =1.46 and β env =1.63 for a max =10 and a max =50 µm, respectively. In Fig. 7 we compare the different β values for each a max with the value obtained from the parametric model. The β values for 0.1µm<a max <100µm are consistent, within the uncertainties, with the β calculated with the parametric model. In the case of grains larger than 100 µm, the total envelope mass is underestimated.

Full radiative transfer model
In this model we used a system that consists of disk, protostar, envelope and outflow cavity. We used the radiative transfer tool RADMC-3D from Dullemond et al. (2012) to compute the emission from all the contributions. The details of each contribution will be discussed in the following sections.

Disk model
We adopt a disk model heated by its protostellar radiation. The surface density profile Σ(R) was modeled as a truncated power law as in Eq. 4, with a power exponent of the surface density distribution p=1; Σ 0 is scaled to accommodate the total mass of the disk M disk . The 2D volume density with an exponential vertical where H p is the pressure scale height and is defined as H p /r=0.1(r/r h p ) φ , r h p is the reference radius set at 25 AU, and φ is the flaring index of the disk, which in this case is set to 1.14, as an average value according to previous studies on young sources (Pineda et al. 2011;Tobin et al. 2013). We used the disk inclination angle, disk radius, and disk mass presented in Table 1.

Envelope model
For the envelope model we adopted a density profile by Tafalla et al. (2002), which combines a power-law behavior for large radius and a central flattening profile at small radius, i.e., n(r) = n 0 1 + (r/r 0 ) α , where n 0 is the central density, r 0 is the radius of the flat region or truncation radius, and α is the asymptotic power index. The outer radius of the envelope is fixed at 8 800 AU to match the beam of Enoch et al. (2009) observations, in which the rescaled envelope mass is 2.2 M . Additionally, since we have evidence of an outflow in this source (Stephens et al. 2017), we included an outflow cavity with an opening angle of 30 • (M. Dunham, priv. comm.) and a lower density of 1.0 × 10 −30 gr cm −3 for the region inside the cavity and the background.

Backwarming effect
The effects of the envelope thermal emission on disk (i.e., backwarming) have been studied in different environment, as in the case of the heavily embedded source L1551 IRS 5 (Butner et al. 1994).
In the case of an envelope around a disk, the millimeter emission of the disk increases. This is because the envelope acts as a thermal cavity, not letting the temperature within the cavity to fall below the temperature of the envelope wall. Therefore, a substantial backwarming effect on the disk can be present depending on the optical depth and geometry of the cavity. In the previous envelope modeling following Miotello et al. (2014), this effect has been ignored due the geometry of the envelope. Different profiles might heat the disk to a different degree. To explore the effects of backwarming we have computed new models which attempt to take it into account. The net effect of the envelope on the disk temperature is discussed in the Appendix D.

Dust opacity
We used two kind of dust opacities in order to test the model. Firstly, we used the opacity computed in Ossenkopf & Henning (1994) based on a coagulated grain size distribution. In this model, a truncated power law is adopted for the initial dust distribution, n(a) ∝ a −q , where the minimum size of the grain is a min =5 nm, the maximum size is a max =250 nm and the power index q is set to 3.5. The dust distribution is calculated after 10 5 years of coagulation with a gas density of n H =10 5 cm −3 expected in a prestellar core. Secondly, we used the previous dust opacities presented in Section 4.2.3. Since the second dust opacity approach covers maximum grain sizes from small grains of 0.1µm, to big grains of 1 cm, we decided to present here the results with those opacities to compare consistently with the previous modeling.

Model fitting
The free parameters for the disk are the outer radius, r out and the disk surface density Σ 0 . The free parameters for the envelope are its mass M env , its power law density profile α, its flattening envelope radius r 0 and its dust opacity, characterized by a env max . The truncation radius of the envelope is set at the outer radius of the disk parameter. Since the disk parameters estimated by Segura-Cox et al. (2016) are not solid constraints, we test our model using their mass and outer radius values as an upper and lower limit on Σ 0 . The grid of parameters that we test and set are presented in Table 8. Once the dust temperature of the system is calculated from the input parameters of Table 8, we compute the synthetic images, for 1.3mm and 2.7mm, following the same procedure reported in Section 4.2.4. We simultaneously fit the 1.3mm and 2.7mm visibilities by calculating χ 2 values for each model using the equation for the entire set of visibility points between 20 and 110 kλ. The uncertainty in the data, σ i , includes the statistical uncertainty and the absolute flux uncertainty of 10% for 2.7 mm data and 20% for 1.3 mm data, both added in quadrature. Since our observational constraints are dominated by the errors of the data sets, it is possible that the disk and/or envelope structure would be wrong at some level, therefore, our χ 2 value is simply an indicator of an acceptable model, not a best fit. After performing a visual inspection of the models, we report the best match with the observations in the next paragraph and in Table 9. A sample of models with different a max and derived envelope masses are presented in the Appendix C.1.

Results full radiative transfer model
From our interferometric observations, we are limited to study the inner regions of the envelope, from 4 000 AU to 600 AU, therefore, we examine a power-law density profile follow Tafalla et al. (2002). An unresolved component is included to represent a compact disk structure. In the envelope we used a constant dust grain population, in which we vary the maximum grain size from 0.1 µm to 1 000 µm. To study the impact of the maximum grain size in the envelope, a env max , the central density, n 0 , and the density power law index, α, we used the range of parameters reported in Table 8. Table 9 shows the model parameters that provides the best fits to the observations. For the disk properties, we compare our results with the values reported in Table 1. Our disk mass and radius are consistent with the rescaled values reported by Segura-Cox et al. (2016). Both disk models in Per-emb-50 are consistent with a small optically thick disk, but do not allow us to probe if there is grain growth throughout the disk since we are missing very long baselines to resolve the disk.
Similar to the results of the two-step model, the full radiative transfer models suggest a distribution of dust grains in the envelope with maximum size a env max = 50, 100 µm and a resulting envelope mass within 8 800 AU radius of M env ∼2.24, 1.54 M , respectively. In Fig. 8, we present a variety of models with different a max in the envelope that match the visibility data. While all the models match the 1.3 mm data within the flux uncertainty (red region), the 2.7 mm data allow us to determine a good model because the shape of the short baseline emission. Models with a env max < 100 µm follow the flat emission of the 2.7 mm data, while models with a env max > 300 µm overestimate the short baseline emission at 2.7mm and underestimate the envelope mass of Table 1. These results are consistent with our previous modeling and in agreement with the spectral index β that we calculated in Fig. 7. We also reported the 1.1 mm single dish flux (see Table 10) for each model.   Notes. Each mass model is calculated within a 8 800 AU envelope radius. The 1.1 mm fluxes are calculated using an aperture of 30", simulating the diameter aperture of Bolocam.
As discussed by many authors (Draine 2006;Banzatti et al. 2011;Testi et al. 2014), it is quite difficult to explain values of β less than 1 without invoking the presence of millimeter size grains, regardless of the chemical composition, porosity, or grain geometry. In the case of Per-emb-50, the high value of β is compatible with grains no larger than 100 µm, and with values found in Class 0 sources by I-Hsiu Li et al. (2017). The fact that we find grains that have not reached mm sizes in the envelope of Per-emb-50 will be discussed in the following section.

Grain sizes in Class I protostellar envelopes
The presence of millimeter-size grains in envelopes of young protostars, Class 0/I, have been studied and modeled by many authors (e.g Kwon et al. 2009;Chiang et al. 2012b;Tobin et al. 2013;Miotello et al. 2014), but current models cannot easily explain growth at that level (Ormel et al. 2009). This is because the models require high number densities, n H > 10 6 cm −3 , to form such large grains in timescales of 1 Myr. Miotello et al. (2014) found that dust grains start to aggregate up to mm sizes already in the envelope of two Class I protostars, producing a change in the spectral index with values of β env =0.6±0.3 for Elias 29 and β env =0.8±0.7 for WL 12. Those values are smaller than the spectral index for the envelope of Per-emb-50, β env =1.4±0.3, by a factor of 2. The differences between these studies may be associated to the properties of the star forming region. In our case Per-emb-50 is in NGC1333 region in Perseus, which is a very crowded region with young stellar objects, while the sources of Miotello et al. (2014) are isolated and embedded in L1688 in Ophiuchus. With this study on Per-emb-50, we suggest the possibility that: (a) millimeter grains in envelopes of young protostars may not be a common result, or (b) the dust grain growth is not a homogeneous process. Finally, the environment within which a protostar forms could also play a role in the amount of dust coagulation, which significantly affects the future formation and structure of the protoplanetary disk, as shown by Zhao et al. (2016Zhao et al. ( , 2018. The possibility that grains can grow up to mm-size in the envelope of Class 0/I protostars was studied by Wong et al. (2016). They proposed another mechanism to explain the existence of mm-size grains in the envelopes of young protostellar sources that consists of transport of mm-sized grains from dense regions close to the protostar to the envelope via the outflow. This scenario is quite plausible before the central mass of the protostar reaches a mass of 0.1 M with a mass-loss rate of 10 −6 M yr −1 . This could be the case of Per-emb-50, but high resolution data is needed to model the inner regions of this source. The results of our analysis show that dust grains may have grown as large as ∼100 µm in size in the envelope of Per-emb-50 at scales of 4000-2000 AU. This implies that there is a degree of grain growth with respect to the ISM sizes, but not significant enough to lower the value of α. This is also in agreement with the work of Chacón-Tanarro et al. (2017), where they predict grain sizes of a few hundred µm in the central 300 AU of the prestellar core L1544.
Taking this into account, it is crucial to perform surveys for Class I protostars embedded in different environments and at different physical scales to determine the variation of α spectral index and the corresponding amount of grain growth.

The effects of Backwarming
We find that backwarming is important for modeling Per-emb-50. From the previous analysis, using the two-step modeling, it was straightforward to fit the nearly constant emission at long baselines with an unresolved disk. For the full radiative transfer modeling, this was not the case. Considering the full radiative transfer model, the use of a Tafalla et al. (2002) density profile combined with a power-law behavior for large radius and a central flattening profile at small radius shows that backwarming is important since the disk emission is completely affected by the addition of the envelope. This change in emission is discussed in Butner et al. (1994) and in the Appendix D.1. To study the effects of different envelopes geometries on disk emission is beyond the scope of this work. However, backwarming can have other consequences. The change of temperature between a backwarmed disk (∼100 K) and a nonbackwarmed disk with ∼20 K, would affect significantly the gas phase chemistry and the dust mantle chemistry in young disks and envelopes (Butner et al. 1994). Finally, the backwarming in Class I protostars could have an important effect in the thermal history of the outer disks of planetary systems. Detailed studies using proper physical structures and radiative transfer models are necessary to address the backwarming effect present in most young embedded sources.

Conclusions
We present new 1.3 mm data from SMA and 2.7 mm data from NOEMA of the brightest Class I protostar Per-emb-50 in the NGC 1333 cluster in the Perseus star forming region. In the uv plane it is possible to distinguish the presence of a large scale envelope at short baselines and an unresolved and optically thick disk at longer u-v distances. From the data analysis and the different modeling approaches on this source we can conclude: -For the envelope uv analysis we find a spectral index similar to the typical ISM values, α mm =3.3±0.3. -The current observations on Per-emb-50 and the radiative transfer modeling reveal a Class I envelope consistent with maximum sized grains of < 100 µm. This suggests that grain growth has proceeded within the envelope, but not to a level to produce changes in α as the presence of millimeter size grains does. -The presence of grains with a size range of <100 µm in envelopes of Class I protostars may have an impact in our understanding of protostellar evolution. Following the prediction from Chacón-Tanarro et al. (2017), who find that dust grains are expected to grow to sizes of a few hundred µm in the central 300 AU of a pre-stellar core, we could suggest that the larger grains found in the envelope of Per-emb-50 may be inherited from the prestellar phase. -These results show for the first time no evidence of grain growth to millimeter sizes in the inner regions of the envelope of a Class I protostar, providing an interesting case for future studies of the efficiency of grain growth process in these stages. -We also explore the effects of backwarming. The analysis shows that the envelope geometry highly affects the disk temperature. In the collapsing envelope model, the effect is weak, but if a power law envelope is used, the effect is more obvious.
Future high sensitivity data will be needed to allow us to conclusively prove whether there are spectral index variations between the disk and the envelope. Moreover, study of a larger sample of Class I sources in different star forming regions is important to understand how general this process is for grain growth.
The spectral index of the observed flux densities α mm can be approximated using the flux density at two wavelengths. In this appendix, we discuss the error propagation from the observational uncertainty to the deduced α mm value. Let F 1 and F 2 be the flux density at frequencies ν 1 and ν 2 , α mm can be expressed as in Equation (1): We assume that the fluxes F 1 and F 2 are independent and that σ F1 and σ F2 are their standard deviations; using the error propagation we obtain Taking the partial derivative of Equation (   The name of the model and derived envelope mass are in the right panel. The color gradient represent the χ 2 from lower (blue) to high values (yellow), that were used only as reference. After visual inspection, we choose the best models from the green area.