Systematic analysis of spectral energy distributions and the dust opacity indices for Class 0 young stellar objects

TL;DR

This study analyzes the dust properties of nine Class 0 young stellar objects using archival data, revealing evidence of early grain growth and discussing its implications for star formation models.

Contribution

It provides a systematic analysis of dust opacity indices in Class 0 YSOs, highlighting potential early grain growth before the Class I stage.

Findings

Four sources show beta < 1.7, indicating possible grain growth.

Remaining sources have beta near 1.7, suggesting grain growth occurs later.

Results support early dust evolution in star formation process.

Abstract

We are motivated by the recent measurements of dust opacity indices beta around young stellar objects (YSOs), which suggest that efficient grain growth may have occurred earlier than the Class I stage. The present work makes use of abundant archival interferometric observations at submillimeter,millimeter, and centimeter wavelength bands to examine grain growth signatures in the dense inner regions (<1000 AU) of nine Class 0 YSOs. A systematic data analysis is performed to derive dust temperatures, optical depths, and dust opacity indices based on single-component modified black body fittings to the spectral energy distributions (SEDs). The fitted dust opacity indices (beta) are in a wide range of 0.3 to 2.0 when single-component SED fitting is adopted. Four out of the nine observed sources show beta lower than 1.7, the typical value of the interstellar dust. Low dust opacity index (or spectral index) values may be explained by the effect of dust grain growth, which makes beta<1.7. Alternatively, the very small observed values of beta may be interpreted by the presence of deeply embedded hot inner disks, which only significantly contribute to the observed fluxes at long wavelength bands. This possibility can be tested by the higher angular resolution imaging observations of ALMA, or more detailed sampling of SEDs in the millimeter and centimeter bands. The beta values of the remaining five sources are close to or consistent with 1.7, indicating that grain growth would start to significantly reduce the values of beta no earlier than the late-Class 0 stage for these YSOs.

Tables

Table 1. Table 1 List of our targets

Source Name	R.A. and Dec.	Region	Distance (pc)	Class
L1448CN	03:25:38.80, +30:45:05.0	Perseus	240	0
NGC1333 IRAS2A	03:28:55.70, +31:14:37.0	Perseus	240	0
NGC1333 IRAS4 (A1, A2)	03:29:10.50, +31:13:31.0	Perseus	240	0
Barnard 1b (N, S)	03:33:21.14, +31:07:35.3	Perseus	240	0
L1527	04:39:53.90, +26:03:10.0	Taurus	140	0/I
Serpens FIRS 1	18:29:49.80, +01:15:20.6	Serpens	415	0/I
L1157mm	20:39:06.19, +68:02:15.9	Isolated	325	0

Table 2. Table 2 Summary of Observations

Source	Frequency	Telescope	UV Coverage	Observing Date
	(GHz)		(k$\lambda$)
L1448CN	37	VLA	10–1290	2013 Oct. 28
	43	VLA	10–1300	1996 Feb. 23, 24, 25
	230	SMA	3–390	2004 Nov. 7
				2011 Sep. 12
				2012 Jan. 7
	340	SMA	21–579	2010 Jul. 19, Sep. 17
NGC1333 IRAS 2A	37	VLA	20–1390	2013 Nov. 4
	43	VLA	5–483	2004 Mar. 28
	100	CARMA	5–125	2008 Oct. 14, 20
	194	SMA	6–50	2008 Aug. 19
	340	SMA	10–256	2010 Oct. 14
				2013 Oct. 21
NGC1333 IRAS4 (A1,A2)	37	VLA	21–1420	2013 Oct. 21
	43	VLA	4–500	2014 Oct. 13
	230	SMA	8–394	2004 Nov. 22
				2006 Jan. 17
				2011 Sep. 12
	340	SMA	10–530	2004 Dec. 05, 06
				2009 Jan. 27, Feb. 24
				2011 Jan. 24, 25
Barbard 1B (N, S)	37	VLA	10–1290	2013 Nov. 7
	100	NMA	0–100	1998 Jan. – May
	230	SMA	5–170	2013 Nov. 19, 23, Dec. 6
				2014 Jan. 11
	280	SMA	3–75	2008 Sep. 03
	340	SMA	12–88	2012 Nov. 10
L1527	43	VLA	10–2000	1996 Feb. 23, 24, 25
				2004 Feb. 11
	90	CARMA	5–560	2009 Oct. 11
				2010 Dec. 02
	230	SMA	4–400	2011 Sep. 02, 09
				2012 Jan. 08
	345	ALMA	12–432	2012 Aug. 29
Serpens FIRS 1	43	VLA	3–426	2004 Mar. 11
	230	SMA	37–386	2009 Jul. 28
	340	SMA	6–45	2012 Jul. 01
L1157	43	VLA	3–492	2001 Nov. 29
				2004 Mar. 14
				2004 Nov. 26
	230	SMA	4–52	2005 Jul. 06
				2009 Sep. 15
	340	SMA	7–79	2005 Jul. 03, Aug. 22

Table 3. Table 3 Gaussian decomposition in the u ​ v 𝑢 𝑣 uv -amp domain

Source	Frequency	Extended Component		Compact Component
		Flux (mJy)	Uncertainty (mJy)	Flux (mJy)	Uncertainty (mJy)
L1448CN	37	0.22	0.06	2.01	0.05
	43	0.73	0.39	1.95	0.02
	230	344.46	15.16	136.43	3.93
	340	275.98	49.28	305.01	6.1
NGC1333 IRAS2A	37	–	–	1.55	0.08
	43	1.18	0.47	1.17	0.08
	100	29.16	2.92	33.04	3.04
	194	12808.92	28263.68	323.17	32.06
	340	826.11	45.34	554.21	34.68
NGC1333 IRAS4A1	37	5.75	0.65	7.83	0.29
	41	7.86	0.51	8.83	0.3
	43	6.4	0.95	11.26	0.54
	45	9.85	0.89	11.81	0.81
	47	10.14	1.29	13.61	1.25
	230	2535.59	114.3	588.62	32.62
	340	2177.26	216.03	897.91	46.18
NGC1333 IRAS4A2	37	5.99	3.01	0.74	0.05
	230	2821.82	165.04	225.48	31.57
	340	2735.19	557.3	539.36	47.32
Barnard 1b-N	37	1.21	0.06	0.3	0.03
	230	67.17	24.88	80.61	7.63
	280	196.26	42.56	232.16	20.17
	340	667.1	445.19	415.05	27.3
Barnard 1b-S	37	0.67	0.05	0.24	0.03
	100	19.04	7.57	18.17	2.54
	230	68.03	38.22	249.89	7.03
	280	1174.25	397.04	378.92	8.15
	340	337.44	101.03	811.71	30.71
L1527	43	2.33	0.8	1.14	0.11
	90	14.37	3.84	13.45	4.05
	230	170.65	6.47	27.05	4.75
	340	283.00	24.58	267.05	26.3
Serpens FIRS1	43	7.33	1.48	1.5	0.23
	230	732.26	123.28	319.47	74.44
	340	3057.06	91.43	1869.91	33.67
L1157	43	6.11	0.99	0.96	0.17
	230	90.68	42	134.82	13.02
	340	412.99	37.42	183.23	9.81

Table 4. Table 4 The results of our SED fitting

Source	$\beta$				Temperature				${\tau }_{1mm}$				${\beta }_{RJ}$
	best fit	median	$16\%$	$84\%$	best fit	median	$16\%$	$84\%$	best fit	median	$16\%$	$84\%$	best fit	median	$16\%$	$84\%$
L1448CN	0.58	0.75	0.56	0.94	102	99	97	103	0.03	0.03	0.03	0.04	0.54	0.71	0.52	0.90
IRAS2A	1.20	1.41	1.16	1.75	36	36	34	37	0.24	0.24	0.17	0.32	1.03	1.22	0.98	1.47
IRAS4A1	0.51	0.72	0.50	0.98	41	41	40	42	0.37	0.33	0.24	0.44	0.34	0.57	0.38	0.81
IRAS4A2	1.23	1.34	0.99	1.62	40	40	38	40	0.18	0.16	0.09	0.21	1.11	1.24	1.00	1.49
B1bN	1.53	1.63	1.36	1.92	18	18	17	19	0.32	0.30	0.21	0.39	1.29	1.40	1.15	1.67
B1bS	2.02	2.19	1.95	2.51	22	22	19	23	0.61	0.60	0.45	0.84	1.72	1.86	1.65	2.18
L1527	0.30	0.63	0.17	1.15	69	58	51	76	0.01	0.02	0.01	0.03	0.26	0.49	0.04	0.91
SerFIRS1	1.74	1.88	1.60	2.27	51	51	49	52	0.39	0.36	0.25	0.48	1.57	1.73	1.43	2.16
L1157	0.82	1.00	0.71	1.48	43	42	41	44	0.06	0.06	0.04	0.08	0.72	0.89	0.60	1.23

Table 5. Table 5 Parameters for three-components SED models

	$T_{di}$	${\Omega }_{di}$	${\kappa }_{\text{230 GHz}}$${\Sigma }^{\text{di}}$	$T_{\text{do}}$	${\Omega }_{\text{do}}$	${\kappa }_{\text{230 GHz}}$${\Sigma }^{\text{do}}$	$T_{\text{ei}}$	${\Omega }_{\text{ei}}$	${\kappa }_{\text{230 GHz}}$${\Sigma }^{\text{ei}}$
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
	(K)	(sr)		(K)	(sr)		(K)	(sr)
L1527	650	1.0$\times$${10}^{-13}$	9.0	12	8.0$\times$${10}^{-13}$	2.8	35	8.6$\times$${10}^{-11}$	6.4$\times$${10}^{-3}$
NGC 1333 IRAS4A1	$\cdots$	$\cdots$	$\cdots$	45	1.2$\times$${10}^{-11}$	13	20	2.0$\times$${10}^{-11}$	1.5
L1157	$\cdots$	$\cdots$	$\cdots$	35	5.8$\times$${10}^{-12}$	1.7	15	6.3$\times$${10}^{-12}$	1.5
L1448CN	$\cdots$	$\cdots$	$\cdots$	40	1.4$\times$${10}^{-12}$	1.5	20	6.6$\times$${10}^{-12}$	1.5
NGC 1333 IRAS2A	$\cdots$	$\cdots$	$\cdots$	150	5.2$\times$${10}^{-12}$	0.39	15	2.0$\times$${10}^{-11}$	1.0

Equations

$$\frac{0.1\times \text{mbox}uv-distance{}^2}{max(\text{mbox}uv-distance)}.$$

$$F_{\nu }=\Omega (1-e^{-{\tau }_{\nu }})B_{\nu }(T),$$

$${\tau }_{\nu }={\tau }_0{\left(\frac{\nu }{{\nu }_0}\right)}^{\beta }\left(\text{mbox}or{\kappa }_{\nu }={\kappa }_0{\left(\frac{\nu }{{\nu }_0}\right)}^{\beta }\right),$$

$$F_{\nu }\propto {\nu }^{2+{\beta }_{RJ}}$$

$${\beta }_{RJ}=(0.89\pm 0.002)\times \beta -(0.03\pm 0.004).$$

$$\begin{array}{r}{F}_{\nu }={\Omega }_{di}(1-e^{-{\tau }_{\nu }^{di}})B_{\nu }(T_{di})e^{-{\tau }_{\nu }^{do}-{\tau }_{\nu }^{ei}}+\\ {\Omega }_{do}(1-e^{-{\tau }_{\nu }^{do}})B_{\nu }(T_{do})e^{-{\tau }_{\nu }^{ei}}+\\ {\Omega }_{ei}(1-e^{-{\tau }_{\nu }^{ei}})B_{\nu }(T_{ei})+\\ {\Omega }_{ff}(1-e^{-{\tau }_{\nu }^{ff}})B_{\nu }(T_{ff})e^{-{\tau }_{\nu }^{di}-{\tau }_{\nu }^{do}-{\tau }_{\nu }^{ei}},\end{array}$$

$$\begin{split}\tau_{\nu}^{ff}=\\ 8.235\times 10^{-2}\left(\frac{T_{e}}{\mbox{K}}\right)^{-1.35}\left(\frac{\nu}{\mbox{GHz}}\right)^{-2.1}\left(\frac{\mbox{EM}}{\mbox{pc\,cm${}^{-6}$}}\right),\\ \end{split}$$

$$\left(\frac{m_{gas}}{M_{\odot }}\right)=4.8\times 10{}^5\left(\frac{\Omega \times {\kappa }_{\text{mbox}230GHz}\Sigma }{{\kappa }_{\text{mbox}230GHz}}\right)\times {\left(\frac{\text{mbox}distance}{\text{mbox}pc}\right)}^2.$$

Full text

1. Introduction

When do the Keplerian rotating protoplanetary disks form? When and how does the maximum dust grain size grow from micron-size to millimeter-size in the protoplanetary disks? These are crucial questions when addressing the formation of planetesimals and planetary systems in star-forming regions. Observationally, the maximum size of dust grains may be constrained by the observations of spectral energy distributions (SEDs) and the dust opacity indices $\beta$ (Draine, 2006). The typical $\beta$ value for dust in the optically thin interstellar medium (ISM) is $\sim$1.7 (e.g., Lin et al. 2016); grain growth up to millimeter sizes will result in $\beta<1$ (e.g., Beckwith & Sargent, 1991; Draine, 2006; Ricci et al., 2010b; Kataoka et al., 2014). Previous observations have shown that some Class I and II sources have $\beta$ values smaller than that of ISM (e.g., Beckwith & Sargent, 1991; Ricci et al., 2010a, b; Miotello et al., 2014). On the other hand, Forbrich et al. (2015) studied the starless core FeSt 1-457 in the pipe nebula and found that $\beta$ in the inner core region shows no significant difference from local clouds, which is similar to the values associated with ISM. Therefore, studying Class 0 YSOs may provide some clues to the starting point when grain growth begins to significantly reduce the value of $\beta$. A comparison of sources at various evolutionary stages will further constrain the timescale of this process.

Derivations of $\beta$ from deeply embedded Class 0 YSOs with observations merely at shorter than $\sim$1 mm wavelengths can be degenerated due to the high optical depth. Jørgensen et al. (2007) have performed a study of 20 Class 0 and I YSOs with the 0.88 and 1.3 mm observations of SMA, and derived spectral indices $\alpha$ around 2.4. Their interpretation was either that significant grain growth has occurred in the observed sources with $\beta$ lower than 2.0, or the dust emission is optically thick and the Rayleigh-Jeans approximation breaks down. Based on CARMA observations, Kwon et al. (2009) and Chiang et al. (2012) found \beta$$\sim1 in L1448 IRS2, L1448 IRS3B and L1157 in the envelope ($>10^{3}$ AU scale). They claimed a radial dependence of $\beta$ in these sources, which might suggest that grain growth is quicker in denser regions. However, most past studies assumed that Class 0 envelopes are optically thin and fall in the Rayleigh-Jeans regime at observed (sub)millimeter wavelengths, which may potentially result in underestimating $\beta$.

In this work, we utilized the archival data of National Radio Astronomy Observatory (NRAO) Very Large Array (VLA) [now the Karl G. Jansky Very Large Array (JVLA)], the Submillimeter Array (SMA), the Combined Array for Research in Millimeter-wave Astronomy (CARMA), the Nobeyama Radio Observatory (NRO) Nobeyama Millimeter Array (NMA), and the Atacama Large Millimeter Array (ALMA), for nine Class 0 YSOs. The good sampling of SEDs from the centimeter to the submillimeter bands, incorporated with archival infrared photometry, allows simultaneously fitting the dust temperature, dust opacity and $\beta$. Observations and our data reduction are described in Section 2. Section 3 summarizes our procedure of data analysis. The results are presented in Section 4. We discuss our results, their implications on grain growth and alternative explanations in Section 5. Our conclusions are given in Section 6.

2. Observations and Data Reduction

In the following sections, we describe how our target sources were selected based on the spatial resolution and spectral band completeness of the archival observations. Also, we briefly summarize the data calibration procedures. Our targets are listed in Table 1. The observations are tabulated in Table 2.

2.1. Target sources

To simultaneously determine dust temperature, optical depth and $\beta$ using a single component modified black body fitting (Section 3.3), we select Class 0 sources that have archival data at no less than three frequency bands in the millimeter and submillimeter wavelength range. In addition, we require the selected sources to present at least one significant detection at $>$3 mm wavelengths. This is because the dust emission tends to be optically thinner at these wavelengths, so that we can obtain a better constrain on the value of $\beta$. Limited by the (low) angular resolution of the archival SMA observations (Section 2.2.4), we only consider nearby targets that are located within 420 pc from the Sun for better spatial resolution.

Based on these criteria, nine well-studied class 0 YSOs are selected and investigated. We quote their distances from Chen et al. (2013), except for Serpens FIRS1 (also known as Serpens SMM1), for which we adopt the most recent measured distance (see Dzib et al., 2010, for further information). Seven out of nine of the sources are classified as Class 0 YSOs, and the remaining two are assumed to be in the transition phase between the Class 0 and Class I stages (Chen et al., 2013; Froebrich, 2005). Basic properties of our samples are summarized in Table 1. We refer the reader to Chen et al. (2013) for other detailed information about these YSOs.

2.2. Submillimeter, millimeter, and radio interferometric data

We retrieve centimeter band observations from VLA, millimeter band observations from VLA, CARMA, NMA and SMA, and submillimeter band observations from SMA and ALMA. Each set of data is calibrated using their official data reduction software packages. We use the AIPS (Greisen, 2002) software package to reduce the VLA data taken before 2010, the CASA (McMullin et al., 2007) software package to reduce (J)VLA data taken after 2010 and ALMA data, the MIR IDL software package (Qi, 2005) for the SMA data reduction, the UVPROC2 software package for the NMA data reduction, and the MIRIAD (Sault et al., 1995) software package for the CARMA data reduction. Finally, MIRIAD is used to evaluate visibility amplitudes as a function of uv distance. When multiple observations are accessible at a certain wavelength for a target source, we combine all the available visibility data in order to maximize the signal-to-noise ratio (SNR).

Some of our targets are binaries. In such cases, we first generate the clean component maps using the MIRIAD task CLEAN. Then we use the MIRIAD task UVMODEL for modeling and subtracting the secondary component from the calibrated visibility data.

A brief summary of each observation and additional treatments in the data reduction is provided below.

2.2.1 VLA Observations

We have performed deep $\sim$43 GHz observations towards NGC1333 IRAS4A using JVLA in 2014 (see Liu et al. (2016) for the details). For all the Perseus sources in our sample, we obtain the public data from the VLA Nascent Disk and Multiplicity (VANDAM) survey (Cox et al., 2015; Tobin et al., 2014), which observed all known protostars in the Perseus molecular cloud in 2013–2014. We only include the upper side band (USB, centered at 37 GHz) data for our analysis, owning to the poor data quality in the lower side band (LSB, centered at 29 GHz). These data taken with the upgraded JVLA have a very high sensitivity and broad bandwidth coverages, and hence allow us to derive good constraints on $\beta$ in SED fitting. We obtain 43GHz data from the VLA data archive for all the sources except B1b (N, S), which has a low SNR. These observations were conducted between 1996 and 2004.

We combined all available array configurations (including A, B, C, BnC and D) with good SNRs. The combined visibility data in general provide visibility data up to few thousands $k\lambda$ uv distances. To compare the (J)VLA data with the SMA data, we trim all the data at $\gtrsim$800 $k\lambda$ uv distances.

2.2.2 CARMA observations

The 90 GHz data for L1527 and 100 GHz data for NGC1333 IRAS2A are obtained from the CARMA data archive. We combine data taken with the A and C array configurations to match the uv distance of other observations for L1527.

2.2.3 NMA observations

The 100 GHz data for Barnard 1b-N/S were observed with the NMA. The details of the observations are described in Hirano & Liu (2014).

2.2.4 SMA observations

We obtain 183 GHz, 230 GHz, 280 GHz and 340 GHz data from the SMA data archive. The antenna configurations include subcompact, compact, extended and very extended. The uv coverage of the observations is typically from a few to few hundred $k\lambda$, which corresponds to ${10}^{3}$-${10}^{5}$AU in physical scale for our targets. We use the MIRIAD task UVLIN to generate the spectral line-free continuum data. For the case that UVLIN cannot provide a good fit to the continuum level, we use the MIRIAD UVAVER routine to average the continuum data over wide channel ranges.

2.2.5 ALMA observations

For L1527, both ALMA and SMA 340 GHz data exist. We choose to use the ALMA data only, because of their better sensitivity and angular resolution. The observation was conducted in ALMA Cycle 0 in 2012 (see, e.g., Sakai et al., 2014, for the details). We separate the continuum data using the CASA task UVCONTSUB.

2.3. FIR data

In order to resolve the degeneracy between $\beta$ and the dust temperature, we include the Spitzer Multiband Imaging Photometer (MIPS) 24 and 70 $\micron$ data to our SED fitting. Considering that the typical dust temperature for our targets is about 10-100 K, these short wavelength observations can cover the high frequency wing in SEDs. Note that the Spitzer observations are lower in spatial resolution (6*′′* at 24 $\micron$ and 18*′′* at 70$\micron$), compared with radio observations. The low resolution of these observations may lead to underestimate temperatures and dust opacity indices. A quantitative discussion is provided in Section 4.2.

3. Flows of Data Analysis Procedures

This section introduces the procedure we used for a systematic data analysis. It is motivated by a widely adopted geometric picture of YSOs (Shu et al., 1987); YSOs can be modeled as spatially compact disk(-like) structures surrounded by large-scale extended envelopes. We first perform an automatized binning for the azimuthally averaged visibility data in the uv distance domain (Section 3.1). We subsequently perform Gaussian decomposition for the binned visibility amplitudes to separate the compact and extended components (Section 3.2). Finally, we conduct the modified black body SED fittings for these components to derive dust temperature, dust opacity, and opacity index (Section 3.3). A number of error sources are considered in the SED fitting, which are described in detail in Section 3.4.

3.1. Visibility Amplitude Analysis: Optimizing Bin Sizes

We dynamically bin the visibility amplitudes to properly weight the observations at different uv distance ranges. We first compute the azimuthally averaged visibility amplitudes at each uv distance using the MIRIAD task UVAMP. Once the visibility amplitudes are obtained, optimal bin sizes are searched such that the following two criteria will be fulfilled: (1) low SNR data points ($<$0.8 of the median value of unbinned data points) are binned with adjacent data points, and (2) long baseline data points with low SNR are binned with wider bins under the condition that the maximum bin width is restricted to the empirically assigned upper limit:

[TABLE]

The former is to ensure good SNRs, while the latter is to avoid over smearing the spatial information.

3.2. Gaussian Decomposition

The typical sizes of circumstellar envelopes and (proto-)disks of Class 0 YSO are of a few thousands AU and $<$100 AU (Williams & Cieza, 2011), respectively. Most of our retrieved radio data cover uv distances of from a few to few hundreds k$\lambda$, which corresponds to ${10}^{5}$-${10}^{2}$AU in physical scale at the distance of $\sim$ 300 pc. Therefore, these observations are, in most cases, insufficient to resolve, or at most resolve marginally the hypothesized compact disk-like component. On the other hand, they are good enough to resolve the extended circumstellar envelops. Nevertheless, the assumed compact and extended components should have very different visibility amplitude distributions, and hence would be possible to systematically decompose the visibility data into these two components. We have inspected all the data of our targets and confirmed that both the compact and the extended components are present clearly in their visibility amplitude distributions.

In the actual decomposition process, we adopt two kinds of fitting profiles: the binned visibility amplitudes are fitted by two Gaussian profiles if the compact component is marginally resolved, while they are fitted by the combination of one Gaussian profile and a point source (i.e. a straight line with a constant amplitude) if the compact component is entirely unresolved. The best fit was searched by evaluating the reduced $\chi^{2}$. As an example, Figure 1 shows the results for L1157. The filled squares with error bars denote the binned visibility amplitudes, and the solid curve is for the best fit. For this target, the compact component is fitted well by the straight line (i.e., a point source), while the extended component is by a Gaussian profile.

We found that the FWHM values derived from Gaussian fitting are highly uncertain and cannot be used for quantitative analysis. The large uncertainties in the size of the compact component originate from the lack of long baseline data that are necessary to resolve the source. Nonetheless, the obtained FWHM is useful for confirming that the emission of the compact component originates from the roughly same physical scale at all the frequency band observations.

3.3. SED fitting

Our assessment for the results of Gaussian decomposition is that the extended components do not have a well defined size and geometry, and are seriously subject to different degrees of missing fluxes at different frequency bands. Therefore, the determination of the total flux for the extended components is highly sensitive to the shortest uv spacing in the interferometric observations. Since we do not have short spacing data from single dish observations for most of our targets, we only make use of the compact components for all the targets in the following SED fitting.

The compact components in general have much higher brightness than the extended components, and their observed fluxes are not sensitive to the missing flux problem. Furthermore, the compact components correspond to much higher density regions of YSOs and thereby a higher degree of grain growth can be expected there, compared with the extended components. Thus, analyzing the compact components is more relevant to our purpose of studying the initial dust grain growth in star-forming regions (e.g., Wong et al. 2016). Therefore, we perform modified black body fittings to the fluxes of the compact components to derive the dust temperature, optical depth and $\beta$ simultaneously.

To proceed, we assume that the flux ($F_{\nu}$) observed at a frequency ($\nu$) can be expressed by the modified black body emission (Hildebrand, 1983):

[TABLE]

where $\Omega$ is the solid angle related to the physical size of a target source ($r=\sqrt{ln2\times\Omega/\pi}$), $\tau_{\nu}$ is the optical depth measured at $\nu$ from the observers towards the source, and $B_{\nu}(T)$ is the Planck function for a certain temperature ($T$) of the source. We assume that compact components at all frequency bands share an identical sold angle of $\Omega$=0$\farcs$5, which can be factored out during the derivation of beta. A simple power-law profile is adopted for $\tau_{\nu}$:

[TABLE]

where $\tau_{0}$ is the optical depth at $\nu=\nu_{0}$. We set $\nu_{0}=$ 300 GHz in our fittings. In the case of spatially uniform $\kappa_{\nu}$, $\tau_{\nu}$=\kappa_{\nu}$$\int\rho d\ell=\kappa_{\nu}$$\Sigma, where $\Sigma$ is dust mass surface density, and $\ell$ is the thickness of structure along the line-of-sight.

In the case where dust emission is optically thin (i.e. $\tau_{\nu}\ll 1$), equation (2) can be simplified according to ${e}^{-\tau_{\nu}}\simeq 1-\tau_{\nu}$. In addition, when the observed frequency bands are in the Rayleigh-Jeans limit ($h\nu/KT\ll 1$), then $B_{\nu}(T)$ becomes approximately proportional to $\nu^{2}$. Consequently, $F_{\nu}$ can be simplified as

[TABLE]

We label $\beta$ as $\beta_{RJ}$ to emphasize that the value of $\beta_{RJ}$ is derived under the assumption of the optically thin emission and the Rayleigh-Jeans limit. Based on these assumptions, equation (4) permits deriving $\beta_{RJ}$ from observations at merely two frequency bands. The approach of $\beta_{RJ}$ has been widely adopted in the literature for examining grain growth. However, when using only two frequency bands (e.g. $\sim$230GHz and $\sim$340GHz), the small separation in the frequency range and uncertainties from flux measurements would yield a large uncertainty in the derivation of beta. For instance, considering same uncertainty estimation as described in Section 3.4, $\beta$ derived from $\sim$230GHz and $\sim$340GHz bands give uncertainties of $>$1 for all sources in our sample, and can not be used as a tool to diagnose grain growth as in literatures.

Also, it should be noted that many Class 0 objects tend to be highly obscured even at (sub)millimeter wavelengths and can have a relatively low average dust temperature (e.g., $<$ 20 K). In these cases, both the optically thin and the Rayleigh-Jeans assumptions break down. This may be particularly significant for the observations at $\sim$340 GHz or higher frequencies. In this work, we hence quantify the difference between the values of $\beta$ and $\beta_{RJ}$ to compare them with each other.

3.4. Uncertainty estimation

To realistically determine the uncertainty of $\beta$ and $\beta_{RJ}$, we perform a Monte Carlo (MC) error estimation. In this procedure, we first estimate the uncertainty of the observed fluxes that is caused by errors arising from a number of sources. In this paper, we consider two plausible sources: the intrinsic error coming from observations and data analyses, and the confusion between the dust and the non-dust emissions. Once the flux uncertainty is determined for these two sources (see below), we then alter the measured flux randomly based on the flux uncertainty estimation. We use the perturbed fluxes to perform SED fitting and obtain the resultant value of $\beta$ and $\beta_{RJ}$. The above procedure is repeated for 1000 times to compute the median values of $\beta$ and $\beta_{RJ}$. We define the median value as the most likely value of $\beta$ and $\beta_{RJ}$ to explore a possibility of dust grain growth in Class 0 YSOs (see Section 4). The median is preferred over the mean because the distributions from MC error estimation are not always symmetric and the median is less affected by the MC realization with extremely steep or flat SEDs. Generally, the difference between the mean and median values are within the uncertainty given for $\beta$ and $\beta_{RJ}$. For the uncertainty, we adopt the 16 and 84 percentile values from the 1000 repetitions as the lower and upper limits. Furthermore, we conduct the reduced $\chi^{2}$ analysis to find out the best fit for $\beta$ under that assumption that free-free emission contributes to 25% of the observed flux at longer wavelengths (see Table 4).

3.4.1 Intrinsic errors

We adopt the flux measurement uncertainty from the Gaussian decomposition procedure ($\sim$10%) as the statistical flux measurement error (not including absolute flux uncertainty). In addition, we consider the nominal, 20% absolute flux calibration uncertainty for all of our observations by adding it onto the statistical flux error. For each MC repetition, we incorporate the intrinsic errors by drawing the perturbed fluxes from Gaussian distributions with the measured flux as the mean and the uncertainties as the standard deviation.

3.4.2 Contribution from non-dust emission

For Class 0/I objects, free-free emission may be produced by the $\sim$100 AU scale thermal radio jets (e.g., Rodriguez, 1997), and occasionally by the spatially unresolved non-thermal flares (Forbrich et al., 2006) for which spectral indices are not yet certain (Liu et al., 2014). The confusion between the dust and the non-dust emissions due to these phenomena may be relatively important at 43 GHz or lower observing frequencies especially for the compact component. Without considering such confusion, SED fitting may lead to underestimating $\beta$. The SED measured at centimeter wavelengths can constrain such confusion. However, this is not always possible when utilizing the archival data which were taken at different times, since the fluxes of radio jets and the non-thermal flares both vary with time (for the characteristic timescale, see Liu et al., 2014).

For a Class 0/I object at $d\sim 140$ pc, the typical radio jet flux may be on the order of 0.1–1 mJy, and the typical non-thermal emission (if present) may be $\sim$0.1 mJy or fainter (e.g., Dzib et al., 2013; Liu et al., 2014; Pech et al., 2016). Our selected target sources appear to be embedded in relatively abundant circumstellar gas and dust, and present brighter 37 GHz and 43 GHz emission than the expected fluxes of optically thin non-dust emission (see Table 3). For the case of L1448 CN, Hirano et al. (2010) found that free-free emission can contribute to up to $\sim$50% of the flux at 30–40 GHz. Based on these observations, we empirically assume in our SED fitting that the non-dust emission may contribute to 0–50$\%$ of the total 37 GHz or 43 GHz fluxes at the epoch of those observations. For each MC repetition, we randomly assign a non-dust emission factor of 0.0–0.5 to the 30–43 GHz fluxes and subtract the computed non-dust emission to evaluate a potential bias in our determination of $\beta$. Since the properties of non-thermal emission in many of our sources are uncertain, the non-dust emission factor is drawn from a uniform distribution. Non-dust emission is negligible at higher frequencies. This correction only lowers the fluxes by a factor of two at most. As a result, (potential) contributions from the non-dust emission yield negligible effects in the derived $\beta$ thanks to the wide frequency and flux range of the SEDs.

4. Results

4.1. Dust opacity indices $\beta$

Table 3 summarizes the fluxes of the spatially compact components that are derived from our Gaussian decomposition (Section 3.2). The results of SED fitting for the compact components are plotted in Figure 2. The obtained values of $\beta$, $\beta_{RJ}$, the dust temperature, and the optical depth are listed in Table 4. A comparison between the values of $\beta$ and $\beta_{RJ}$ is shown in Figure 3.

Our results show that both $\beta$ and $\beta_{RJ}$ distribute widely, with the average values of $\beta=1.3$ and $\beta_{RJ}=1.1$ (Figure 3). The derived $\beta$ of SerFIRS 1, B1-bN, IRAS 2A and IRAS 4A2 are comparable to the interstellar value of $\sim$1.7, while the $\beta$ of B1-bS is slightly larger than 2.0. For the rest of our targets, the resultant values of $\beta$ are $<$1.0. In addition, Figure 3 shows that the derived $\beta_{RJ}$ and $\beta$ tend to be well correlated with each other. A linear regression gives

[TABLE]

We also find that the fitted optical depths at \lambda$$=1 mm range from 0.02 to 0.61 and the resultant dust temperatures are in the range of from $\sim$10 K to $\sim$100 K for our targets (see Table 4).

The fitted values of $\beta$ for our samples are broadly in agreement with the previous studies. For instance, our results indicate that B1-bN and B1-bS have larger values of $\beta$ and the lowest dust temperatures in our targets, which are both consistent with the previous estimate done by Hirano & Liu (2014). As other examples, the previous derivations based on CARMA observations have suggested that L1157 has \beta$$\sim0.5 (Kwon et al., 2009; Chiang et al., 2012), and that L1527 has \beta$$\sim0 (Tobin et al., 2013). Jørgensen et al. (2007) also investigated the dust opacity indices for L1448C, IRAS 2A, IRAS 4A, IRAS 2B, L1157 and L1527 using the SMA data at 0.88 and 1.3 mm, and found that the calculated $\beta$ are around 0–1 at long baselines ($>$ 40k$\lambda$) for all their sources. Such low values of $\beta$ are consistent with our single-component fittings.

We caution that the derived, lower than 1.7 values of $\beta$ here do not immediately imply the evidence of grain growth. The values of $\beta$ can be underestimated when the assumption of a single uniform emission component breaks down. In fact, a signature of this possibility can be seen in Figure 2; the observed $\sim$340 GHz fluxes of L1448CN, IRAS 2A, IRAS 4A1, and L1157 are smaller than the regression lines derived from the optically thin and Rayleigh-Jeans assumptions. If such differences are not entirely due to measurement errors, then these differences suggest that the dust emission at $\gtrsim$230 GHz should already become optically thick for these targets. This is however inconsistent with the low $\tau_{\mbox{\scriptsize{1\,mm}}}$ values given by the full single-component modified black body fittings, especially for the cases of L1448CN and L1157 (Table 4). In addition, we have to admit that the observed SED of L1527 is not fitted well by either the optically thin, Rayleigh-Jeans approximation or the full modified black body fitting with the single emission component (see Figure 2). Furthermore, while some of the previous studies concluded that the small $\beta$ value indicates (a significant degree of) grain growth, Jørgensen et al. (2007) also pointed out that resolved observations of the compact region and modeling of the envelope and disk structure are necessary to put tighter constrains on $\beta$ (see Section 4.2 for more quantitative discussion). Finally, the recent, higher angular resolution VLA observations show that the $\beta$ of IRAS 4A1 and 4A2 are much closer to $\sim$1.7 and $\sim$1.2 between 8 and 10 millimeter wavelengths, respectively (Cox et al., 2015). It is obvious that there is a big difference between their estimate and the one done by us and Jørgensen et al. (2007) for IRAS 4A1. We will discuss below a potential origin of why such inconsistencies emerge (see Section 4.3).

4.2. Degeneracy between dust temperature and $\beta$

FIR data points are included in our SED fitting to provide a better constraint on the estimate of the dust temperature ($T$). However, due to the low spatial resolution of available observations at these wavelengths, we cannot distinguish the compact and extended components for the observations, as done for our radio/submillimeter data. The mixture of the flux from both the components can result in underestimating $T$ and $\beta$ for the compact component. We quantify this effect on $\beta$ by fitting SEDs at centimeter to submilimeter wavelengths with $T$ fixed at 100K. Note that the value of $T=100$ K is a conservative upper limit since the averaged dust temperature around these sources is generally $<$60K. We find that for this case $\beta$ is underestimated only by $\sim$10$\%$ compared to $\beta_{RJ}$, which is on the same scale of the our error estimation. Thus, we can conclude that the low spatial resolution observations at FIR do not affect our estimate of $\beta$ very much.

4.3. Effects of multiple emission components to the $\beta$ measurements

Figure 3 shows that grain growth has not yet significantly lowered the values of $\beta$ for five of the nine selected target sources. Here, we construct a toy model to discuss whether or not the small values of the derived $\beta$ for the remaining four sources (L1448CN, IRAS4A1, L1157, L1527) indicate dust grain growth. We will also attempt to interpret some features in the observed SEDs of IRAS2A, which are difficult to be explained only by the effects of grain growth or optical depth when assuming single emission component.

Many of the observed sources show unexpectedly low spectral slopes ($\alpha$) in certain frequency ranges, which correspond to \alpha$$<2.0. Examples include the observations of L1527 in between 90 GHz and 230 GHz, which is consistent with the data presented in Tobin et al. (2013). Others are L1157, IRAS2A, IRAS4A1, in between $\sim$200 GHz and 345 GHz (see Table 3). Some low values of $\alpha$ ($<$2.0) were also reported by Jørgensen et al. (2007) and Miotello et al. (2014).

One possibility to explain the low $\alpha$ values may be a very low dust temperature (e.g., $\lesssim$10 K), which shifts SED peaks to submillimeter bands. For example, Hirano & Liu (2014) have demonstrated that superimposing the fluxes of an 11 K extended dust component with a 23 K compact dust component may explain the observed SED of Barnard 1-bN on $\sim$2*′′* angular scales. For our present observations of the compact emission of L1527, IRAS 4A1, L1157, and L1448CN, it is unlikely that the dominant emission at any band is significantly contributed from a $\lesssim$10 K dust component, given the detected high brightness on small spatial scales. It is strictly impossible for the case of IRAS 4A1, since the observed brightness temperature at 44 GHz with $\sim$0$\farcs$5 resolution is already 41 K (Liu et al., 2016). A high dust temperature may also present in the inner region of IRAS2A, given the $\gtrsim$100 K gas excitation temperature reported in the previous spectral line observations (e.g., Coutens et al. 2015). Therefore, it is less likely that the dominant contribution of the (sub)millimeter band continuum emission of IRAS2A is from $<$10 K dust, in spite of the observed steeply decreasing spectral index in between 230 GHz and 345 GHz frequencies. From the astrochemistry point of view (e.g., Sakai et al. 2014), L1527 appears warm (30–90 K) from the 0$\farcs$5-1*′′* resolution observations, and may require at least three density/temperature components to explain the observed abundance distributions. The previous Herschel observations of high-$J$ (J_{u}$$>13) CO lines towards the Class 0 object BHR 71 also indicated a very hot ($\sim$400-1000 K) gas component near the embedded protostar (Yang et al., 2017).

On the other hand, an obscured, inner hot dust component may only significantly contribute to the observed fluxes at long wavelength bands when it is embedded by an optically thicker, cooler outer component. This can be a mechanism to produce an excess of strong emission at long wavelength bands, which is alternative to the common interpretation with significant dust grain growth (i.e., small $\beta$ values), very cold dust, or very optically thick dust. In realistic Class 0 YSOs, such a hypothesized geometric picture may be realized as an inner hot (pseudo)disk of few AU scales, obscured by an optically thick outer cool (pseudo)disk to few tens of AU. For a deeply embedded, high density accretion flow, important heating mechanisms can include compression work, viscous heating, and radiative heating. However, most of the previous radiative transfer modeling for fitting observed SEDs only consider radiative heating. An example of temperature and column density profile for embedded YSOs based on numerical simulations can be found in Vorobyov et al. (2013) and reference therein. Furthermore, Class 0 YSOs may be further embedded within the more extended, optically thinner inner envelope. Such hypothesized geometric and temperature profiles may be particularly relevant for the case of L1527, since its geometry has been resolved approximately as an edge-on accretion flow (Tobin et al., 2013). The rotating disk was also identified from the high angular resolution VLA observations of NH3 lines towards IRAS 2A (Choi et al., 2010). Here we cross reference to a parallel work reporting a spatially resolved case of self-obscured disk at 345 GHz, which is approximately in an edge-on projection (Lee et al., 2017).

As motivated by the above consideration, we undertake SED fitting with a toy model of the multiple emission components. In the model, we suppose that an inner hot pseudo-disk with a scale of a few AU is obscured by an outer, optically thick, cool pseudo-disk with a scale of a few tens of AU. After incorporating the potential contribution of free-free emission from an ionized disk or a jet ($F^{ff}_{\nu}$), the overall integrated flux of such a system may be approximated as

[TABLE]

where $\Omega_{ff,di,do,ei}$ are the solid angles of the free-free emission source, the inner disk, outer disk, and inner envelope components (assuming \Omega_{ei}$$>$$\Omega_{do}$$>$$\Omega_{di}$$>$$\Omega_{ff}); $\tau^{ff,di,do,ei}_{\nu}$ and $T_{ff,di,do,ei}$ are the optical depths and dust temperatures of these components. Following Keto (2003) and Mezger & Henderson (1967), we approximate $\tau_{\nu}^{ff}$ by:

[TABLE]

where $T_{e}$ is the electron temperature, and EM is the emission measure defined as EM=$$\int n_{e}^{2}d\ell, with $n_{e}$ being the electron number density. Here we assume that the free-free emission source is obscured by all dust components. The following discussion will not be sensitive to whether or not the free-free emission source is obscured.

Our goal is to test if (a) the SEDs of L1448CN, IRAS41, L1157, L1527 can be explained without having to assume a significant affectation by dust grain growth, and (b) if the $\alpha<$2.0 spectral index of IRAS2A at 230–345 GHz can be explained in terms of obscuration. We construct toy SED models for all these sources based on Equation 6 and assuming $\beta$=2.0. The parameters we used are listed in Table 5 and the results are presented in Figure 4. Notice that the number of free parameters in our toy model is larger than the number of measurements we have. Therefore, it is impossible to avoid a degeneracy in the product of the dust opacity at 230 GHz and column density, \kappa_{\mbox{\tiny{230 GHz}}}$$\Sigma and the dust temperature, T. Nevertheless, we argue that these toy models only use the minimum numbers of dust or free-free emission components and thus the minimum numbers of free parameters that are required to explain the observed SED. For example, two dust emission components are required to reproduce the \alpha$$<2.0 results of L1527 between 90 GHz and 230 GHz, and yet another optically thinner dust component is needed to reproduce its high flux at 345 GHz. For this particular case, the SED peak of the hot inner disk component (i.e., Figure 4, dash-dot cyan curve) is shifted to a lower frequency than what is expected from a single emission component, modified black body profile (Equation 2, 3) due to the obscuration by other components. For L1448CN, L1157, IRAS2A, and IRAS4A1, a cool inner envelope component is required to obscure an embedded warm disk component to reproduce the \alpha$$<2.0 results in between 230 GHz and 345 GHz. Similarly, the SED peaks of the outer disk components of these four sources (i.e., Figure 4, red dotted curve) are shifted to the lower frequencies, due to the obscuration by the envelope components. The SEDs of these four sources can be modeled without considering the innermost hot disk component. The contribution of the innermost hot disk component is uncertain for L1448CN due to the confusion by free-free emission. We thus assume that the hot inner disk component is negligible in our present toy models. We caution that in order to minimize the number of free parameters, our toy models may be yet over simplified. In reality, dust temperature may have rather continuous spatial distributions.

We find that our simplified toy models can reasonably reproduce the observed SEDs of all five sources. This implies that the actual $\beta$ may be consistent with 2.0, which is higher than the fitted values presented in Figure 2, 3, and Table 4. Assuming a gas-to-dust mass ratio of 100.0, the implied gas mass of individual dust components can be approximated by

[TABLE]

The estimated integrated gas masses for L1527, IRAS 4A1, L1157, L1448CN, and IRAS2A are 0.084/\kappa_{\mbox{\tiny{230GHz}}}$$M_{\odot}, 5.2/\kappa_{\mbox{\tiny{230GHz}}}$$M_{\odot}, 0.99/\kappa_{\mbox{\tiny{230GHz}}}$$M_{\odot}, 0.34/\kappa_{\mbox{\tiny{230GHz}}}$$M_{\odot}, and 1.0/\kappa_{\mbox{\tiny{230GHz}}}$$M_{\odot}, respectively. Assuming \kappa_{\mbox{\tiny{230GHz}}}$$\sim1, the estimated gas masses are reasonable for young, low-mass star-forming sources.

5. Discussion

5.1. Alternative explanation for low $\beta$

If the simplified geometric picture and temperature profile of our toy SED models are reasonable, then the previously reported, radially increasing $\alpha$ (Chiang et al., 2012) may be alternatively interpreted by that the embedded hot dust component progressively becomes more prominent towards the projected central region (i.e., becoming warmer or denser), instead of by radially increasing $\beta$. On the other hand, perhaps not all the inner dust components are fully obscured by the outer ones. When dust emission from the inner region is not obscured (e.g., when the disk components are approximately in face-on projection), the sparsely sampled SED will look similar to a single modified black body profile unless the temperature changes drastically.

Comparing the toy models presented in Figure 4 with the systematic fittings presented in Figure 2 leads us to think that when the SEDs cannot be well represented by a single modified emission profile (Equation 2), the full single-component modified black body fittings does not provide a more physical derivation of $\beta$ than the derivation based on the optically thin and Rayleigh-Jeans assumption. Instead, they may converge to equally poor (or good) fits. One may face the situation that the assumption of single dominant emission component breaks down when a too wide frequency range is considered, and that the high frequency emission and the low frequency emission do not necessarily trace the same emission component (e.g., Figure 4). The $\beta$ value of the dominant contributor for low frequency (e.g., $\lesssim$50 GHz) emission may be more accurately constrained when there are two independent measurements from the optically thinner part of the spectrum, which may be (close to) the case of IRAS 4A2 reported by Cox et al. (2015).

5.2. Implications for dust grain growth in Class 0 YSOs

We are now in a position to discuss implications of our results that are derived from the SED analysis for our targets.

Our SED fitting with the single emission component indicates that the resultant value of $\beta$ becomes comparable to the value of the interstellar dust for five of the nine targets (IRAS 2A, IRAS4A2, B1bN, B1bS, and SerFIRS1, see Figure 3). This suggests that dust grain growth is unlikely to proceed significantly in these targets. For the remaining four sources (L1448CN, NGC 1333 IRAS4A1, L1157, L1527), the single component SED fitting shows a lower ($<1$) value of $\beta$, which can be attributed to a higher degree of dust growth. However, by invoking a multiple emission component scenario, we can account for the observed SEDs reasonably well under the condition that $\beta=2$. Thus, we argue that the values of $\beta$ and $\beta_{RJ}$ for these four targets should be regarded as lower limits. As a consequence, our presented SED analysis does not provide firm evidence of significant dust grain growth to reduce the value of $\beta$ down to an order of unity or less in our nine Class 0 YSO targets. Note that the additional information of the target sources, especially time variability, would be valuable to examine how our interpretations can be affected by such effects. All kinds of pure radiative transfer simulations/modeling (including our model) generally assume equilibrium states for gas and dust in the targets.

How is our interpretation compatible with the previous studies? A lower degree of dust growth in Class 0 YSOs is indeed in favor for the results of theoretical studies. Some studies show that the collisional growth of dust grains generally does not proceed up to the grain size of $\sim 1$mm due to the low density, short-lived environment there (e.g., Ormel et al., 2009; Hirashita & Li, 2013). Even if some degrees of grain growth can occur, the number density of larger grains may be low (e.g., Wong et al., 2016) and/or the resultant dust particles would be porous aggregates (e.g., Ormel et al., 2007; Okuzumi et al., 2012), so that a signature of grain growth may not be captured well by our SED fitting. Thus, our interpretation can hold in the current literature.

The above interpretation enables in turn, to imply that dust grain growth may not become significant enough to lower the value of $\beta$ in a stage earlier than late-Class 0. On the other hand, multiple studies have reported that the values of $\beta$ already go down to an order of unity or less at the Class II stage. Comprehensive surveys (Ricci et al., 2010a, b) on 38 flux selected class II sources in the $\rho$-Ophiuchi and Taurus-Auriga star forming region (SFR) revealed $\langle\beta\rangle$=0.5–0.6. They fitted the SEDs of resolved sources at the wavelengths of 0.5–3mm with a radiative transfer model of two-layers disks, proposed by Chiang & Goldreich (1997), and concluded that maximum grain size $a_{max}$ $>$ 1cm is necessary to explain the derived $\beta$ in both SFR considering a reasonable range of dust opacity and power-law size distribution (see Figure 3 of Ricci et al., 2010a). The distribution of $\beta_{RJ}$ of their Class II sources and our sources are compared in Figure 5. Taking into account these results, it may be crucial to perform a survey for Class I YSOs to detect the evolutionary history of $\beta$ due to grain growth.

Finally, Figure 5 shows a relative uniform distribution of the $\beta$ derived from our systematic analysis of the SED fitting for our Class 0 smaple. If this flat distribution would arise from the effect of the multiple emission components, then this trend may imply that our Class 0 YSOs harbor randomly orientated disk-like structures, such that the hot inner parts of only some of them are obscured (at shorter than millimeter wavelength bands), or that the heating mechanism(s) for the innermost disk is episodic (e.g. Vorobyov et al. 2013). These hypotheses may be tested by observing a statistical sample of Class 0 and Class I objects, and by exploring their gas velocity profiles.

6. Conclusions

We have analyzed the archival SMA, (J)VLA, NMA, CARMA, and ALMA data for nine well studied Class 0 YSOs, which cover $\sim$30 GHz to $\sim$345 GHz frequency ranges. We systematically perform single-component modified black body fittings to their SEDs at radio and (sub)millimeter wavelengths, and find that five of them show a dust opacity spectral index \beta$$\sim1.7, which implies no significant signatures of dust grain growth.

We cannot rule out that dust grain growth has significantly lowered the values of $\beta$ in the remaining four sources. However, our radiative transfer toy models assuming \beta$$=2.0 show that their SEDs can alternatively be interpreted by the presence of hot inner dust components that are obscured by the cooler and optically thicker outer components at short wavelengths (e.g., $<$1 mm). In fact, this mechanism is favored to explain the spectral slope (alpha) values lower than 2.0 observed in certain frequency ranges of some sources. Therefore, the presented data in this work do not show robust evidence that a higher degree of dust grain growth has taken place in any of our observed sources. By comparing with the previous observations of Class II sources, we suggest that dust grain growth may significantly reduce the values of $\beta$ in circumstellar (pseudo-)disks or envelopes no earlier than the late-Class 0 stage but before the Class II stage of YSOs. A future survey of radio and (sub)millimeter SEDs for Class I YSOs may improve the constraints on the timescales of grain growth and the starting point where a lower ($\lesssim 1$) value of $\beta$ is achieved, and thereby improving our understanding of planet-forming environments and the formation of planetesimals.

The authors thank the anonymous referee for useful comments on our manuscript. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics, and is funded by the Smithsonian Institution and the Academia Sinica (Ho et al., 2004). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.00210.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. Part of this research was carried out at Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. Y.H. is supported by JPL/Caltech. N.H. is supported by the grant from Ministry of Science and Technology of Taiwan (MoST 105-2112-M-001-026).

Reference

• Beckwith & Sargent (1991)

Beckwith, S. V. W., & Sargent, A. I. 1991, ApJ, 381, 250

• Chen et al. (2013)

Chen, X., Arce, H. G., Zhang, Q., et al. 2013, The Astrophysical Journal, 768, 110

• Chiang & Goldreich (1997)

Chiang, E., & Goldreich, P. 1997, The Astrophysical Journal, 490, 368

• Chiang et al. (2012)

Chiang, H.-F., Looney, L. W., & Tobin, J. J. 2012, The Astrophysical Journal, 756, 168

• Choi et al. (2010)

Choi, M., Tatematsu, K., & Kang, M. 2010, ApJ, 723, L34

• Coutens et al. (2015)

Coutens, A., Persson, M. V., Jørgensen, J. K., Wampfler, S. F., & Lykke, J. M. 2015, A&A, 576, A5

• Cox et al. (2015)

Cox, E. G., Harris, R. J., Looney, L. W., et al. 2015, The Astrophysical Journal Letters, 814, L28

• Draine (2006)

Draine, B. 2006, The Astrophysical Journal, 636, 1114

• Dzib et al. (2010)

Dzib, S., Loinard, L., Mioduszewski, A. J., et al. 2010, The Astrophysical Journal, 718, 610

• Dzib et al. (2013)

Dzib, S. A., Loinard, L., Mioduszewski, A. J., et al. 2013, ApJ, 775, 63

• Forbrich et al. (2015)

Forbrich, J., Lada, C. J., Lombardi, M., Román-Zúñiga, C., & Alves, J. 2015, A&A, 580, A114

• Forbrich et al. (2006)

Forbrich, J., Preibisch, T., & Menten, K. M. 2006, A&A, 446, 155

• Froebrich (2005)

Froebrich, D. 2005, The Astrophysical Journal Supplement Series, 156, 169

• Greisen (2002)

Greisen, E. W. 2002, in Information Handling in Astronomy-Historical Vistas (Springer), 109–125

• Hildebrand (1983)

Hildebrand, R. H. 1983, QJRAS, 24, 267

• Hirano et al. (2010)

Hirano, N., Ho, P. P. T., Liu, S.-Y., et al. 2010, ApJ, 717, 58

• Hirano & Liu (2014)

Hirano, N., & Liu, F.-c. 2014, ApJ, 789, 50

• Hirashita & Li (2013)

Hirashita, H., & Li, Z.-Y. 2013, MNRAS, 434, L70

• Ho et al. (2004)

Ho, P. T. P., Moran, J. M., & Lo, K. Y. 2004, ApJ, 616, L1

• Jørgensen et al. (2007)

Jørgensen, J. K., Bourke, T. L., Myers, P. C., et al. 2007, ApJ, 659, 479

• Kataoka et al. (2014)

Kataoka, A., Okuzumi, S., Tanaka, H., & Nomura, H. 2014, A&A, 568, A42

• Keto (2003)

Keto, E. 2003, ApJ, 599, 1196

• Kwon et al. (2009)

Kwon, W., Looney, L. W., Mundy, L. G., Chiang, H.-F., & Kemball, A. J. 2009, The Astrophysical Journal, 696, 841

• Lee et al. (2017)

Lee, C.-F., Li, Z.-Y., Ho, P. T. P., et al. 2017, Science Advances, 3

• Lin et al. (2016)

Lin, Y., Liu, H. B., Li, D., et al. 2016, ApJ, 828, 32

• Liu et al. (2014)

Liu, H. B., Galván-Madrid, R., Forbrich, J., et al. 2014, ApJ, 780, 155

• Liu et al. (2016)

Liu, H. B., Lai, S.-P., Hasegawa, Y., et al. 2016, ApJ, 821, 41

• McMullin et al. (2007)

McMullin, J., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Data Analysis Software and Systems XVI, Vol. 376, 127

• Mezger & Henderson (1967)

Mezger, P. G., & Henderson, A. P. 1967, ApJ, 147, 471

• Miotello et al. (2014)

Miotello, A., Testi, L., Lodato, G., et al. 2014, Astronomy & Astrophysics, 567, A32

• Okuzumi et al. (2012)

Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106

• Ormel et al. (2009)

Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. 2009, A&A, 502, 845

• Ormel et al. (2007)

Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215

• Pech et al. (2016)

Pech, G., Loinard, L., Dzib, S. A., et al. 2016, ApJ, 818, 116

• Qi (2005)

Qi, C. 2005, Harvard-Smithonian Center for Astrophysics

• Ricci et al. (2010b)

Ricci, L., Testi, L., Natta, A., & Brooks, K. J. 2010b, Astronomy & Astrophysics, 521, A66

• Ricci et al. (2010a)

Ricci, L., Testi, L., Natta, A., et al. 2010a, A&A, 512, A15

• Rodriguez (1997)

Rodriguez, L. F. 1997, in IAU Symposium, Vol. 182, Herbig-Haro Flows and the Birth of Stars, ed. B. Reipurth & C. Bertout, 83–92

• Sakai et al. (2014)

Sakai, N., Oya, Y., Sakai, T., et al. 2014, ApJ, 791, L38

• Sault et al. (1995)

Sault, R. J., Teuben, P. J., & Wright, M. C. 1995, in Astronomical Data Analysis Software and Systems IV, Vol. 77, 433

• Shu et al. (1987)

Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23

• Tobin et al. (2013)

Tobin, J. J., Hartmann, L., Chiang, H.-F., et al. 2013, The Astrophysical Journal, 771, 48

• Tobin et al. (2014)

Tobin, J. J., Dunham, M. M., Looney, L. W., et al. 2014, The Astrophysical Journal, 798, 61

• Vorobyov et al. (2013)

Vorobyov, E. I., Baraffe, I., Harries, T., & Chabrier, G. 2013, A&A, 557, A35

• Williams & Cieza (2011)

Williams, J. P., & Cieza, L. A. 2011, arXiv preprint arXiv:1103.0556

• Wong et al. (2016)

Wong, Y. H. V., Hirashita, H., & Li, Z.-Y. 2016, PASJ, 68, 67

• Yang et al. (2017)

Yang, Y.-L., Evans, II, N. J., Green, J. D., Dunham, M. M., & Jørgensen, J. K. 2017, ArXiv e-prints, arXiv:1701.00803