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

TL;DR

This study uses multi-wavelength millimeter observations and modeling to determine dust grain sizes in the inner envelope of a Class I protostar, revealing grains no larger than 100 micrometers, which informs early star formation and planet formation theories.

Contribution

First to constrain dust grain sizes in the inner envelope of a Class I protostar using combined observational and radiative transfer modeling.

Findings

Envelope dust grains are no larger than 100 micrometers.

The envelope's structure influences millimeter emission from the disk.

Dust spectral index in the envelope is similar to ISM values.

Abstract

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.

Figures

Tables

Table 1. Table 1: Parameters from literature.

Source	Per-emb-50	New value	Ref
RAJ2000	03:29:07.76	–	1
DecJ2000	+31:21:57.2	–	1
$L_{\mathrm{bol}}$ $\left(L_{\odot }\right)$	$10\pm$3.0	13.7$\pm$3.0	1
$M_{\mathrm{env}}$ $\left(M_{\odot }\right)$	$1.62\pm$0.16	2.2$\pm$0.16	1
$\mathrm{PA}$ (deg)	$170\pm$0.3	–	2
$i$ (deg)*	$67\pm$10	–	2
$M_{\mathrm{disk}}$ $\left(M_{\odot }\right)$	0.18 – 0.36	0.28 – 0.58	2
$R_{\mathrm{disk}}$ (au)	21.9 – 25.7	27.3 – 32.1	2
F1.1mm (mJy)	$612\pm 18$	–	3

Table 2. Table 2: Summary of Observations

Observatory	Representative	Date	Flux	Bandpass	Flux	Array	Synthesized	P.A.
	Frequency		Calibrator	Calibrator		Configuration	Beam
	(GHz)				(Jy)		(′′)	(∘)
SMA	220	Nov 27 2014	Uranus	3c84	11.64	Subcompact	1.2$\times$0.96	86.9
				3c454.3	16.71	Subcompact
		Sep 15 2015	Uranus	3c84	11.64	Extended
		Oct 29 2016	Neptune	3c84	13.8	Extended
				3c454.3
NOEMA	109$$$$footnotemark:	Nov 6 2016	0333+321	3c454.3	14.03	8C	2.2$\times$1.7	35
		Nov 12 2016	0333+321	3c84	24.80	6C

Table 3. Table 3: Best Parametric Model

Fit parameters
$F_e$ (mJy)	$17.0\pm$1.1
${\alpha }_e$	$3.3\pm$0.3
${\alpha }_d$	$1.7\pm$0.3
$\sigma$ (k$\lambda$)	$43.1\pm$0.5

Table 4. Table 4: Derived parameters from Parametric model

$F_d^{2.7\mathrm{mm}}$ [mJy]	$18.4\pm 0.7$
$F_{\mathrm{zero}}^{2.7\mathrm{mm}}$ [mJy]	$19.8\pm 2.2$
$F_{\mathrm{zero}}^{1.3\mathrm{mm}}$ [mJy]	$81.3\pm 2$
$F_{\mathrm{zero}}^{1.1\mathrm{mm}}$ [mJy]	$127.4\pm 2$

Table 5. Table 5: Two–step model grid parameters

Parameter	Description	Values	Parameter Use
Stellar model parameters
d (pc)	Distance	293	fixed
$L_{\star }$ ($L_{\odot }$)	Photosphere luminosity	13.7	fixed
$T_{\star }$ (K)	Effective temperature	5 011	fixed
$R_{\star }$ (AU)	Stellar radius	0.025	fixed
$M_{\star }$ ($M_{\odot }$)	Stellar mass	2.9	fixed
Disk model parameters
$R_{\mathrm{in}}$ (AU)	Disk inner radius	0.1	fixed
$R_{\mathrm{out}}$ (AU)	Disk outer radius	25 27 30 32 34 36	varied
$M_{\mathrm{disk}}$ ($M_{\odot }$)	Disk mass	0.05 0.1 0.15 0.2 0.3 0.4	varied
${\Sigma }_{\mathrm{disk}}$ (gr ${\mathrm{cm}}^{-2}$)	Disk surface density	54.7–908	varied
$a_{\max }^{disk}$ ($\mu$m)	Disk maximum grain size	500 1 000 5 000 10 000 20 000	varied
RADMC-3D / Envelope parameters
$r_{\mathrm{in}}$ (AU)	Envelope inner radius	25 27 30 32 34 36	varied
$r_{\mathrm{out}}$ (AU)	Envelope outer radius	8 800	fixed
$R_{\mathrm{rot}}$ (AU)	Centrifugal radius	100–1000	varied
${\rho }_0$ (gr ${\mathrm{cm}}^{-3}$)	Density in the equatorial plane at $R_{\mathrm{rot}}$	0.5$\times {10}^{-20}$–20.0$\times {10}^{-20}$	varied
$a_{\max }^{env}$ ($\mu$m)	Envelope maximum grain size	0.1–1 000	varied

Table 6. Table 6: Two–step Model best-fit parameters

Parameter	Description	Best–fit M1	Best–fit M2
Disk model parameters
$R_{\mathrm{out}}$ (AU)	Disk outer radius	32	34
$M_{\mathrm{disk}}$ ($M_{\odot }$)	Disk mass	0.4	0.2
${\Sigma }_{\mathrm{disk}}$ (gr ${\mathrm{cm}}^{-2}$)	Disk surface density	554.13	245.38
$a_{\max }^{\mathrm{disk}}$ ($\mu$m)	Disk maximum grain size	10 000	10 000
RADMC-3D / Envelope parameters
$R_{\mathrm{rot}}$ (AU)	Centrifugal radius	600	600
${\rho }_0$ (gr ${\mathrm{cm}}^{-3}$)	Density in the equatorial plane at $R_{\mathrm{rot}}$	¡8.5$\times {10}^{-20}$	¡6.0$\times {10}^{-20}$
$a_{\max }^{\mathrm{env}}$ ($\mu$m)	Envelope maximum grain size	¡100	¡100

Table 7. Table 7: Derived envelope masses

Model	$a_{\max }$ ($\mu m$)	$M_{\mathrm{env}}$ ($M_{\odot }$)
M1	0.1,10	1.73
	50	1.53
	100	1.15
	300	0.23
	1000	0.34
M2	0.1	1.22
	10,50	1.38
	100	0.85
	300	0.16
	1000	0.29

Table 8. Table 8: Full radiative transfer model grid parameters

Parameter	Description	Values	Parameter Use
Stellar model parameters
$M_{\star }$ ($M_{\odot }$)	Stellar mass	2.9	fixed
$R_{\star }$ ($R_{\odot }$)	Stellar radius	5.0	fixed
$T_{\star }$ (K)	Effective temperature	5011	fixed
Disk parameters
${\Sigma }_{bkg}$ (gr ${\mathrm{cm}}^{-3}$)	Background density	$1.0\times {10}^{-30}$	fixed
$M_{\mathrm{dust}}/M_{\mathrm{gas}}$	Dust-to-gas mass ratio	0.01	fixed
$R_{H_p}$ (AU)	Reference radius at which $H_p/R$ is taken	25	fixed
$m_{\mathrm{disk}}$ ($M_{\odot }$)	Mass of the disk	0.18–0.36	varied
$\varphi$	Flaring index	1.14	fixed
$p$	Power exponent of the surface density distribution	1.0	fixed
$r_{\mathrm{out}}$ (AU)	Disk outer radius	25,27,30,32	varied
$r_{\mathrm{in}}$ (AU)	Disk inner radius	1.0	fixed
$a_{\max }^{\mathrm{disk}}$($\mu$m)	Disk maximum grain size	10000	fixed
Envelope parameters
$R_{\mathrm{out}}$ (AU)	Envelope outer radius	8,800	fixed
$\alpha$	Power exponent of the radial density distribution	-1.1,-1.5,-1.8	varied
$n_0$ (gr ${\mathrm{cm}}^{-3}$)	Central density	$1.0,1.5,2.0\times {10}^{-16}$	varied
$r_0$ (AU)	Within this radius the density profile is flat	25,27,30,32	varied
$\theta$ (∘)	Opening angle of the outflow	30	fixed
$a_{\max }^{\mathrm{env}}$($\mu$m)	Envelope maximum grain size	0.1–1000	varied

Table 9. Table 9: Full radiative transfer best fit models

Parameter	Description	Best-Model 1	Best-Model 2
Disk parameters
mdisk ($M_{\odot }$)	Mass of the disk	0.20	0.24
${\Sigma }_{\mathrm{disk}}$ (gr ${\mathrm{cm}}^{-2}$)	Disk surface density	362.2	364.2
$r_{\mathrm{out}}$ (AU)	Disk outer radius	25	27
$a_{\max }^{\mathrm{disk}}$($\mu$m)	Disk maximum grain size	10 000	10 000
Envelope parameters
$\alpha$	Power exponent of the radial density distribution	-1.1	-1.1
${\rho }_0$ (gr ${\mathrm{cm}}^{-3}$)	Central density	$2.0\times {10}^{-16}$	$1.5\times {10}^{-16}$
$M_{\mathrm{env}}$ ($M_{\odot }$)	Envelope mass	2.24	1.54
$r_{\mathrm{trun}}$ (AU)	Truncation radius	25	27
$r_0$ (AU)	Within this radius the density profile is flat	25	27
$a_{\max }^{\mathrm{env}}$($\mu$m)	Envelope maximum grain size	50	100

Table 10. Table 10: Derived 1.1  mm mm {\rm mm} fluxes & envelope mass

$a_{\max }$ ($\mu$m)	$M_{\mathrm{env}}$ ($M_{\odot }$)	$F_{1.1\mathrm{mm}}$ (Jy)
0.1	2.24	1.87
50	2.24	1.76
100	1.54	3.18
300	0.04	0.37
1000	0.08	0.35

Equations

$${\alpha }_{\mathrm{mm}}=\frac{\ln {\mathrm{F}}_1-\ln {\mathrm{F}}_2}{\ln {\nu }_1-\ln {\nu }_2}$$

$$f($uv$,\mathrm{\nu})=F_{e}\left(\frac{\mathrm{\nu}}{\nu_{\mathrm{1.3mm}}}\right)^{\alpha_{e}}\mathrm{exp}\left({\frac{-{(\mathrm{uv})}^{2}}{2\sigma^{2}}}\right)+F_{d}\left(\frac{\nu}{\nu_{\mathrm{1.3mm}}}\right)^{\alpha_{d}},$$

$$R_{\mathrm{eff}}={\left(\frac{L_{\mathrm{bol}}}{4\pi \sigma T_{\mathrm{eff}}^4}\right)}^{1/2},$$

$$\Sigma (R)={\Sigma }_0{\left(\frac{r_{out}}{r_{{\Sigma }_0}}\right)}^{-p},$$

$${\rho }_{\mathrm{env}}(r,\theta )={\rho }_0{\left(\frac{R_{\mathrm{rot}}}{r}\right)}^{3/2}{\left(1+\frac{\cos \theta }{\cos {\theta }_0}\right)}^{-1/2}{\left(\frac{\cos \theta }{2\cos {\theta }_0}+\frac{R_{\mathrm{rot}}}{r}\mathrm{co}{\mathrm{s}}^2{\theta }_0\right)}^{-1},$$

$$\frac{r(\cos {\theta }_0-\cos \theta )}{(R_{\mathrm{rot}}\cos {\theta }_0{\sin }^2{\theta }_0)}=1$$

$$\rho (r,z)=\frac{\Sigma (r)}{H_p\sqrt{2\pi }}\exp \left(-\frac{z^2}{2H_p^2}\right),$$

$$n(r)=\frac{n_0}{1+(r/r_0{)}^{\alpha }},$$

$${\chi }^2=\sum _{i=1}^N\frac{(F_{\nu ,\mathrm{observed},i}-F_{\nu ,\mathrm{model},i}{)}^2}{{\sigma }_i^2},$$

$${\alpha }_{\mathrm{mm}}=\frac{\ln {\mathrm{F}}_1-\ln {\mathrm{F}}_2}{\ln {\nu }_1-\ln {\nu }_2}.$$

$$\sigma_{\mathrm{\alpha}}^{2}=\Bigg{|}\frac{\partial\alpha}{\partial\mathrm{F_{1}}}\Bigg{|}^{2}\sigma_{F1}^{2}+\Bigg{|}\frac{\partial\alpha}{\partial\mathrm{F_{2}}}\Bigg{|}^{2}\sigma_{F2}^{2}\,$$

$$\frac{\partial \alpha }{\partial {\mathrm{F}}_1}=\frac{1}{(\ln {\nu }_1-\ln {\nu }_2){\mathrm{F}}_1},$$

$$\frac{\partial \alpha }{\partial {\mathrm{F}}_2}=\frac{1}{(\ln {\nu }_1-\ln {\nu }_2){\mathrm{F}}_2}.$$

$$\sigma_{\mathrm{\alpha}}^{2}=\Bigg{(}\frac{1}{\mathrm{ln}\mathrm{\nu}_{1}-\mathrm{ln}\mathrm{\nu}_{2}}\Bigg{)}^{2}\Bigg{(}\frac{\sigma_{F1}^{2}}{\mathrm{F_{1}^{2}}}+\frac{\sigma_{F2}^{2}}{\mathrm{F_{2}^{2}}}\Bigg{)}.\,$$

$$X_k(t+1)=X_j+Z(X_k(t)-X_j),$$

Full text

*Aims. *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.

*Methods. *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.

*Results. *We find a spectral index in the envelope of Per-emb-50 of $\alpha_{\rm env}$=$3.3\pm 0.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}}$$\leq100 $\mu$m. Additionally, we explore the backwarming effect, where we find that the envelope structure affects the millimeter emission of the disk.

*Conclusions. *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.

1 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 (Testi et al. 2014).

Thanks to detailed studies of protoplanetary disks at several sub-mm 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 $\mu$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 timescales of \sim$$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$$\mum.

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 adversely affect the growth from sub-micron particles, and large dust aggregates could form (Birnstiel et al. 2010).

The advantage of studying protostars at millimeter wavelengths is that the dust emission from the envelope and the disk is mostly optically thin. In this wavelength range, the dust opacity coefficient $\kappa_{\mathrm{\nu}}$, can be approximated by a power law $\kappa_{\mathrm{\nu}}\propto\nu^{\beta_{\mathrm{mm}}}$, where $\beta_{\mathrm{mm}}$ is the millimeter dust opacity spectral index, and is directly related to the maximum size of the grain (Natta et al. 2007). In the presence of very large grains, much larger than the observing wavelength, the opacity becomes gray (only the geometrical cross section of the grains is relevant) and $\beta_{\mathrm{mm}}$=  0. Values of $\beta_{\mathrm{mm}}$ can be estimated by measuring the slope $\alpha_{\mathrm{mm}}$ of the sub-mm spectral energy distribution (SED), $F_{\mathrm{\nu}}\propto\nu^{\alpha_{\mathrm{mm}}}$. When the Rayleigh–Jeans approximation is applicable, the spectral index of the observed flux densities, $\alpha_{\mathrm{mm}}$, would translate to a power law index of the dust opacity $\beta_{\mathrm{mm}}=\alpha_{\mathrm{mm}}$ - 2. While values around $\alpha_{\mathrm{mm}}\sim$3.7 represent size distribution similar to interstellar medium (ISM) particles (Natta et al. 2007; Testi et al. 2014). Classical protoplanetary disks around Class II objects, present clear signs of dust coagulation, with \alpha_{\mathrm{mm}}$$\leq 3 (Testi et al. 2014, and references therein).

Previous observations of Class 0 protostars by Chiang et al. (2012a), Jørgensen et al. (2007) and Kwon et al. (2009), indicate spectral indexes $\alpha_{\mathrm{mm}}\sim 3$, which is shallower than the ISM, but not quite as steep as Class II disks. However, Class 0 objects, are affected by the presence of powerful accretion of material from the envelope and jets, e.g. (Tobin et al. 2013), making them difficult to observe and model. In contrast, Class I protostars have less massive envelopes, which provides a more cleaner analysis of the dust properties since the envelope and disk emission can be separated.

Here we present a dual wavelength analysis and modeling on the Class I protostar Per-emb-50, in the Perseus star forming region. Observations and data reduction are described in Sec. 2. In Sec. 3 we present our observational analysis. The modeling and discussion are presented in Section 4 and 5 respectively. Conclusions and future work are in Section 6.

2 Observations and Data reduction

2.1 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 Tbol=254$\pm$23 K. The rescaled bolometric luminosity is Lbol=13.7$\pm$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: Mdisk=0.28–0.58 M*⊙* and rout=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.

2.2 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 $\tau_{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. (2018).

The quasars 3C454.3 and 3C84 were used as bandpass and phase calibrators. The absolute flux was calibrated on Uranus, with $\sim$20% of flux calibration uncertainty. For the purpose of this work, we use the 1.3 mm data in the Subcompact and Extended array configurations, with projected baselines in the range of 23-119 k$\lambda$. The resulted combined beam was 1.′′72 $\times$ 1.′′40 at P.A. 50.80∘.

2.3 NOEMA observations

The 2.7 mm data presented in this work were obtained with NOEMA, the IRAM222http://www.iram-institute.org/ NOrthern Extended Millimeter Array. The observations were performed on November 6th and 12th, 2016. The array was in the C compact configuration, with 8 antennas (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*333Stations with * correspond to antennas not available for the second day track, N17 and E04. The projected baselines were from 7.8 k$\lambda$ and 102 k$\lambda$.

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 $\times$ 1.′′6 at P.A. 34.84∘. The clean map has a $rms$ noise level of 2.1 mJy beam*-1*.

3 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 ($\mathrm{uv}$-distance). The deprojected uv-distances are given by $R$=$\sqrt{\smash[b]{d_{a}^{2}+d_{b}^{2}}}$, where $d_{a}$=$\sqrt{u^{2}+v^{2}}\sin\phi$ and $d_{b}$=$\sqrt{u^{2}+v^{2}}\cos\phi\cos i$, $\phi$=$\arctan(v/u)-PA$ (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$\lambda$ (see Fig.2) for both wavelengths. We assume that the emission from those baselines belongs to the embedded disk, where the average values at 2.7 mm and 1.3 mm are: $F_{\mathrm{disk}}^{\mathrm{2.7mm}}$=$18.82\pm 0.13$ mJy and $F_{\mathrm{disk}}^{\mathrm{1.3mm}}$=$63.85\pm 4.2$ mJy.

The spectral index $\alpha_{\mathrm{mm}}$ can be calculated through the flux ratio between two wavelengths,

[TABLE]

Using the fluxes between the u–v ranges 47–80 k$\lambda$ at 1.3 and 2.7 mm, we obtain therefore the average value $\alpha_{\mathrm{mm}}$ in the unresolved disk, which is $\alpha_{\mathrm{disk}}$=1.71$\pm$0.3.

As shown in the Fig.2, an excess of emission is present at short baselines ($<$47 k$\lambda$) 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: $F_{\mathrm{ex}}^{\mathrm{2.7mm}}$=$0.52\pm 0.1$ mJy and $F_{\mathrm{ex}}^{\mathrm{1.3mm}}$=$10.1\pm 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 $\alpha_{\mathrm{env}}$ bigger than the typical ISM values, $\alpha_{\mathrm{env}}$=$4.0\pm 0.8$.

The uncertainty in $\alpha_{\mathrm{disk}}$ and $\alpha_{\mathrm{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 $\alpha_{\mathrm{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 \beta_{\mathrm{env}}$$\sim2.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 $\alpha_{\mathrm{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 $\alpha_{\mathrm{mm}}$.

4 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 $\alpha$ 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.

4.1 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 $\mathrm{uv}$ distance, defined as $\mathrm{\sqrt{\mathrm{u}^{2}+\mathrm{v}^{2}}}$, and frequency, $\nu$, is described by:

[TABLE]

where $\mathit{F_{\mathrm{e}}}$ and $\mathit{F_{\mathrm{d}}}$ are the flux density from the Gaussian emission (extended envelope) and point source emission (unresolved disk) respectively, $\alpha_{\mathrm{e}}$ and $\alpha_{\mathrm{d}}$ are the spectral indexes of the two components, and $\sigma$ is the width of the Gaussian given by $\sigma\approx\mathrm{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, $\mathit{F_{\mathrm{d}}}$=63 mJy. Then, four parameters are explored: $F_{\mathrm{e}}$, $\alpha_{\mathrm{e}}$, $\alpha_{\mathrm{d}}$ and $\sigma$. 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.

4.1.1 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 $\alpha_{\mathrm{mm}}$ for the envelope and disk using a simplistic approximation for non correlated errors is $\pm 0.3$.

Additionally, from this simple model we can constrain the size of the region where the envelope emission arises, which the 1-sigma width (from Table 3) 43 k$\lambda$ (1405 AU). From the model we can derive the flux from the disk at 2.7 $\mathrm{mm}$ and the prediction of the total flux at baseline=0 k$\lambda$ or zero spacing flux, $\mathit{F^{2.7\mathrm{mm}}_{\mathrm{zero}}}$ and $\mathit{F^{1.3\mathrm{mm}}_{\mathrm{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, $\mathit{F_{\mathrm{zero}}^{\mathtt{1.1mm}}}$. We find the flux is only 127.4 mJy beam*-1* which is much lower than the single dish flux of 612$\pm$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.

4.2 Two–step Model

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).

4.2.1 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_{\star}$, effective temperature $T_{\mathrm{eff}}$, and mass $M_{\star}$. To obtain $T_{\mathrm{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_{\rm bol}$ reported in Table 1, we estimate $T_{\rm eff}$=5011 K. With $L_{\rm bol}$ and $T_{\rm eff}$, we can estimate $R_{\rm eff}$ using the Stefan-Boltzmann law:

[TABLE]

Then, with $R_{\rm eff}$=5.01 $R_{\odot}$, we use the mass vs. radius relation for a spherical protostar accreting at a rate of $10^{-5}$ $M_{\odot}$ yr*-1* from Palla & Stahler (1991), to deduce an effective mass of $M_{\rm eff}$=2.9 $M_{\odot}$ (Table. 5). Additionally, we add a disk structure defined by an inner and outer radius, $r_{\mathrm{in}}$ and $r_{\mathrm{out}}$, an inclination angle $i$, and a dust surface density profile that follows a simple power law,

[TABLE]

where $\Sigma_{0}$ is the surface dust density fixed at r${}_{\Sigma_{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_{\mathrm{in}}$, we set $r_{\mathrm{in}}$=0.1 AU. $r_{\mathrm{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.

4.2.2 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,

[TABLE]

where $\rho_{\mathrm{0}}$ is the density in the equatorial plane at the centrifugal radius $R_{\mathrm{rot}}$ of the envelope, and $\theta_{\mathrm{0}}$ is the solution of the parabolic motion of an infalling particle given by:

[TABLE]

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 $\rho_{0}$by imposing a total envelope mass Menv, 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.

4.2.3 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$) $\propto$ $a$-q, between a minimum and a maximum grain size, $a_{\mathrm{min}}$ and $a_{\mathrm{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_{\mathrm{min}}$=0.01 $\mu$m and we use q=3.0. We varied $a_{\mathrm{max}}^{\mathrm{disk}}$ and $a_{\mathrm{max}}^{\mathrm{env}}$ according to the range presented in Table 5.

4.2.4 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 Mdisk, $R_{\mathrm{out}}$ and $a_{\mathrm{max}}$ to reproduce together $F_{\mathrm{disk}}^{\mathrm{1.3mm}}$ and $F_{\mathrm{disk}}^{\mathrm{2.7mm}}$. Once we found the three parameters that match F${}_{\mathrm{disk}}^{\mathrm{1.3mm}}$=63.97 mJy and F${}_{\mathrm{d}}^{2.7\mathrm{mm}}$=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_{\mathrm{env}}$ and $a_{\mathrm{max}}^{\mathrm{env}}$ 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.

4.2.5 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 M1 with a 32 AU disk radius and $M_{\mathrm{disk}}=0.4\rm M_{\odot}$ is consistent with the rescaled values reported by Segura-Cox et al. (2016). While all the disk models match the long baselines 1.3 mm data, the disk emission at 2.7 mm is 15% lower than the data. On the other hand, the disk model M2, with a 34 AU disk radius and $M_{\mathrm{disk}}=0.2\rm M_{\odot}$ matches very well the observations at both wavelengths, but compared with values of Table 1, the disk radius is slightly larger.

The differences in the disk models may be due to the assumed values of $\kappa_{\mathrm{\nu}}$= 0.00146 $\mathrm{cm}^{2}$ $\mathrm{g}^{-1}$ and disk temperatures of 20 K and 40 K in Segura-Cox et al. (2016). Therefore, higher resolution millimeter observations that would resolve the disk are needed to put much stronger constraints on Per-emb-50.

For the envelope, we explore the effects of changing: $R_{rot}$, $a_{\mathrm{max}}$ and $\rho_{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_{\mathrm{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 $\mu$m¡ $a_{\mathrm{max}}$ $<100$ $\mu$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.2\rm M_{\odot}$ envelope mass derived by Enoch et al. (2009) are those derived from models with dust grain sizes of $a_{\mathrm{max}}\leq 50\mu\mathrm{m}$ (see Table 6). The models with $a_{\mathrm{max}}=100\leavevmode\nobreak\ \mu\mathrm{m}$ recover almost 60% of the envelope mass. Table 7 presents the derived masses for the envelope using different $a_{\mathrm{max}}$ in M1 and M2.

A distribution of grains with a_{\mathrm{max}}$$\leq$$50 $\mu\mathrm{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 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 $\alpha$.

For the envelope we can use our dust model to infer the value of $\beta$, which is $\beta_{\mathrm{env}}$=1.46 and $\beta_{\mathrm{env}}$=1.63 for $a_{\mathrm{max}}$=10 and $a_{\mathrm{max}}$=50 $\mu$m, respectively. In Fig. 7 we compare the different $\beta$ values for each $a_{max}$ with the value obtained from the parametric model. The $\beta$ values for 0.1$\mu$m¡$a_{\mathrm{max}}$¡100$\mu$m are consistent, within the uncertainties, with the $\beta$ calculated with the parametric model. In the case of grains larger than 100 $\mu$m, the total envelope mass is underestimated.

4.3 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.

4.3.1 Disk model

We adopt a disk model heated by its protostellar radiation. The surface density profile $\Sigma(R)$ was modeled as a truncated power law as in Eq. 4, with a power exponent of the surface density distribution p=1; $\Sigma_{0}$ is scaled to accommodate the total mass of the disk $M_{\mathrm{disk}}$. The 2D volume density with an exponential vertical profile is defined by:

[TABLE]

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 $\phi$ 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.

4.3.2 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.,

[TABLE]

where $n_{0}$ is the central density, $r_{0}$ is the radius of the flat region or truncation radius, and $\alpha$ 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\times 10^{-30}$ $\mathrm{gr\,cm^{-3}}$ for the region inside the cavity and the background.

4.3.3 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.

4.3.4 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$) $\propto$ $a^{-q}$, where the minimum size of the grain is $a_{\mathrm{min}}$=5 nm, the maximum size is $a_{\mathrm{max}}$=250 nm and the power index q is set to 3.5. The dust distribution is calculated after 105 years of coagulation with a gas density of $n_{\rm H}$=105 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$\mu$m, to big grains of 1 cm, we decided to present here the results with those opacities to compare consistently with the previous modeling.

4.3.5 Model fitting

The free parameters for the disk are the outer radius, $r_{\mathrm{out}}$ and the disk surface density $\Sigma_{0}$. The free parameters for the envelope are its mass $M_{env}$, its power law density profile $\alpha$, its flattening envelope radius $r_{\mathrm{0}}$ and its dust opacity, characterized by $a_{\mathrm{max}}^{\mathrm{env}}$. 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 $\Sigma_{\mathrm{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 $\chi^{2}$ values for each model using the equation

[TABLE]

for the entire set of visibility points between 20 and 110 k$\lambda$. The uncertainty in the data, $\sigma_{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 $\chi^{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 10.

4.3.6 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 $\mu$m to 1 000 $\mu$m. To study the impact of the maximum grain size in the envelope, $a_{\mathrm{max}}^{\mathrm{env}}$, the central density, $n_{0}$, and the density power law index, $\alpha$, 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_{\mathrm{max}}^{\mathrm{env}}=50,100$ $\mu$m and a resulting envelope mass within 8 800 AU radius of M_{\mathrm{env}}$$\sim2.24, 1.54 M*⊙*, respectively. In Fig. 8, we present a variety of models with different $a_{\rm 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_{\mathrm{max}}^{\mathrm{env}}$ ¡ 100 $\mu$m follow the flat emission of the 2.7 mm data, while models with $a_{\mathrm{max}}^{\mathrm{env}}$ ¿ 300 $\mu$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 $\beta$ that we calculated in Fig. 7. We also reported the 1.1 mm single dish flux (see Table 10) for each model.

As discussed by many authors (Draine 2006; Banzatti et al. 2011; Testi et al. 2014), it is quite difficult to explain values of $\beta$ 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 $\beta$ is compatible with grains no larger than 100 $\mu$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.

5 Discussion

5.1 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_{\mathrm{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 $\beta_{\mathrm{e}nv}$=0.6$\pm$0.3 for Elias 29 and $\beta_{\mathrm{e}nv}$=0.8$\pm$0.7 for WL 12. Those values are smaller than the spectral index for the envelope of Per-emb-50, $\beta_{\mathrm{e}nv}$=1.4$\pm$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. (2016, 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_{\mathrm{\odot}}$ with a mass-loss rate of 10*-6* $M_{\mathrm{\odot}}$ 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 $\sim$100 $\mu$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 $\alpha$. This is also in agreement with the work of Chacón-Tanarro et al. (2017), where they predict grain sizes of a few hundred $\mu$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 $\alpha$ spectral index and the corresponding amount of grain growth.

5.2 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 11. 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 ($\sim$100 K) and a nonbackwarmed disk with $\sim$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.

6 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 u–v 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 $\mathrm{uv}$ analysis we find a spectral index similar to the typical ISM values, $\alpha_{\mathrm{mm}}$=3.3$\pm$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 $\mu$m. This suggests that grain growth has proceeded within the envelope, but not to a level to produce changes in $\alpha$ as the presence of millimeter size grains does.

• •

The presence of grains with a size range of ¡100 $\mu$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 $\mu$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.

Acknowledgements.

C.A.G. acknowledges support from CONICYT-Becas Chile (grant 72160297). P.C, J.P and L.S acknowledge the financial support of the European Research Council (ERC; project 320620). L.T. acknowledges the financial support of the Italian Ministero dell’Istruzione, Università e Ricerca through the grant Progetti Premiali 2012 – iALMA (CUP C52I13000140001), and by the Deutsche Forschungs-gemeinschaft (DFG, German Research Foundation) – Ref no. FOR 2634/1 TE 1024/1–1. MT has been supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC-2013-ADG. A.M. acknowledges an ESO Fellowship.

Appendix A Error estimate of the spectral index of the observed flux densities

The spectral index of the observed flux densities $\alpha_{\mathrm{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 $\alpha_{\mathrm{mm}}$ value. Let $\mathrm{F_{1}}$ and $\mathrm{F}_{2}$ be the flux density at frequencies $\nu_{1}$ and $\nu_{2}$ , $\alpha_{\mathrm{mm}}$ can be expressed as in Equation (1):

[TABLE]

We assume that the fluxes $\mathrm{F_{1}}$ and $\mathrm{F}_{2}$ are independent and that $\sigma_{F1}$ and $\sigma_{F2}$ are their standard deviations; using the error propagation we obtain

[TABLE]

Taking the partial derivative of Equation (A.1), we obtain

[TABLE]

[TABLE]

Substituting equation A.3 and A.4 in equation A.2, the uncertainty of the derived $\alpha_{\mathrm{mm}}$ is then:

[TABLE]

Appendix B $emcee$ implementation

To compute the posterior distribution for all the free parameters, we use a variant of the Markov Chain Monte Carlo (MCMC) (Mackay 2003; Press et al. 2007) algorithm, which is widely known and efficient in finding a global maximum for a range of posteriors. We follow the affine-invariant ensemble sampler for MCMC by Goodman & Weare (2010), which basically transforms highly anisotropic and difficult-to-be-sampled multivariate posterior probability distribution function (PDFs) into isotropic Gaussians. The immediate advantage is that it is possible to simultaneously run many Markow chains ($walkers$) that will interact in order to converge to the maximum of the posterior.

This algorithm involves an ensemble $S=\left\{X_{k}\right\}$ of simultaneously evolving $K$ walkers, where the transition distribution for each walker is based on the current position of the other $K-1$ walkers belonging to the complementary ensemble $S_{k}=\left\{X_{j},\ \forall j\neq k\right\}$. The position of a walker $X_{k}(t)$ is updated as follows:

[TABLE]

where $X_{j}\in S_{k}$ and $Z$ is a random variable drawn from a distribution that does not depend on the covariances between the parameters.

In this study we adopted an ensemble of 400 walkers, and let MCMC evolve for an initial burn–in phase. The burn–in phase is needed to allow MCMC to perform a consistent sampling of the space of parameters and to find the posterior maximum. To achieve the posterior maximum is needed to introduce the term: autocorrelation-time888Note: The longer the autocorrelation time, the larger the number of the samples we must generate to obtain the desired sampling of the posterior PDF., which is a direct measurement of the number of the posterior PDF evaluations needed to produce independent samples of the target density. For the analysis of Per-emb-50, 750 burn-in steps were performed to achieve convergence. Fig. B.1 presents a staircase plot, using the Python module $corner$ by Foreman-Mackey (2016), showing the marginalized and bi-variate probability distributions resulting from the fit for Per-emb-50.

Appendix C Full Radiative transfer models

The models presented in this appendix were created using a simple python module to set up RADMC-3D for disk plus envelope systems, SimpleDiskEnv999https://gitlab.mpcdf.mpg.de/szucs/SimpleDiskEnv. Fig. C.1. show the best 36 models from a total of 288 for each maximum grain size (0.1,50,100,300,1000 $\mu$m).

Appendix D Backwarming effect

We studied the net effect of the envelope on the disk temperature using a RADMC-3D toy model of a Class I protostar. As mentioned in Butner et al. (1994), the envelope can have an important backwarming effect on the disk, affecting the outer edges of the disk with a flat temperature distribution.

To probe this effect, we first modeled a disk of 25 AU without an envelope and with a distribution of dust grains in the disk with a maximum size $a_{max}^{disk}$=1 cm. Then we add a 1.3 M*⊙* envelope, with a Tafalla et al. (2002) density profile and grain sizes with $a_{max}^{env}$=100 $\mu$m. The inner edge of the envelope and the outer radius of the disk are the same. To compare and quantify the effect, we model a disk with the same characteristics but with a density profile of a collapsing envelope defined by Ulrich (1976). Fig D.1. shows the temperature structure (in cylindrical coordinates) for both of these cases. In the left panels (disk only) we can see that the outer regions of the disk are around 20–30 K. In the right upper panel (disk+envelope) using the Ulrich envelope structure, the temperature increases to 40–60 K. In the case of the model with a Tafalla et al. envelope structure the effect is quite strong, reaching disk outer temperatures of 120–140 K.