THE JCMT GOULD BELT SURVEY: EVIDENCE FOR DUST GRAIN EVOLUTION IN PERSEUS STAR-FORMING CLUMPS

Michael Chun-Yuan Chen; J. Di Francesco; D. Johnstone; S. Sadavoy; J. Hatchell; J. C. Mottram; H. Kirk; J. Buckle; D. S. Berry; H. Broekhoven-Fiene; M. J. Currie; M. Fich; T. Jenness; D. Nutter; K. Pattle; J. E. Pineda; C. Quinn; C. Salji; S. Tisi; M. R. Hogerheijde; D. Ward-Thompson; P. Bastien; D. Bresnahan; H. Butner; A. Chrysostomou; S. Coude; C. J. Davis; E. Drabek-Maunder; A. Duarte-Cabral; J. Fiege; P. Friberg; R. Friesen; G. A. Fuller; S. Graves; J. Greaves; J. Gregson; W. Holland; G. Joncas; J. M. Kirk; L. B. G. Knee; S. Mairs; K. Marsh; B. C. Matthews; G. Moriarty-Schieven; C. Mowat; S. Pezzuto; J. Rawlings; J. Richer; D. Robertson; E. Rosolowsky; D. Rumble; N. Schneider-Bontemps; H. Thomas; N. Tothill; S. Viti; G. J. White; J. Wouterloot; J. Yates; M. Zhu

doi:10.3847/0004-637X/826/1/95

1. INTRODUCTION

Thermal dust emission is an excellent tracer of cold, star-forming structures. When observed in wide bands, the spectral energy distributions (SEDs) of these structures are dominated by optically thin, thermal dust emission at far-infrared and sub-millimeter wavelengths. The column densities, and consequently the masses, of star-forming structures can thus be estimated from their dust emission by assuming a dust opacity, κ_ν, and a dust temperature, T_d. The dust opacity at frequencies ν ≲ 6 THz, however, has a frequency dependency typically modeled as a power law, characterized by the emissivity spectral index, β (e.g., Hildebrand 1983). Depending on the dust model, κ_ν can be uncertain by up to a factor of 3–7 at a given frequency between 0.3 and 3 THz (see Figure 5 in Ossenkopf & Henning 1994). The ability to precisely determine mass therefore depends heavily on how well β and the normalization reference opacity, ${\kappa }_{{\nu }_{0}}$ , can be constrained.

Measuring β is difficult because it requires the SED shape of the dust emission to be well determined, especially when T_d is not already known from independent measurements. The ability to determine the true SED shape can be particularly susceptible to noise in the data and temperature variations along the line of sight, resulting in erroneous measurements of β and T_d (e.g., Shetty et al. 2009a, 2009b; Kelly et al. 2012). For this reason, a fixed β value of ∼2 has commonly been adopted in the literature, motivated by both observations of diffuse interstellar medium (ISM; e.g., Hildebrand 1983) and grain emissivity models (e.g., Draine & Lee 1984).

Given that β is an optical property of dust grains and can depend on various grain properties such as porosity, morphology, and surface composition, the β value of a dust population may be expected to change as the population evolves. Indeed, observations have shown that dust within protostellar disks have β ≃ 1 (e.g., Beckwith & Sargent 1991), which is substantially lower than the diffuse ISM value (∼2). Furthermore, a wide range of β values (1 ≲ β ≲ 3) has been reported in many observations of star-forming regions on smaller scales (≲2'; e.g., Friesen et al. 2005; Shirley et al. 2005, 2011; Kwon et al. 2009; Schnee et al. 2010; Arab et al. 2012; Chiang et al. 2012), as well as a few lower resolution observations on larger scales of a cloud (e.g., Dupac et al. 2003; Planck Collaboration et al. 2011). These results suggest that dust grains in star-forming regions can evolve significantly away from its state in the diffuse ISM, prior to being accreted onto protostellar disks.

Recent large far-infrared/sub-millimeter surveys of nearby star-forming clouds such as the Herschel Gould Belt Survey (GBS; André et al. 2010) and the James Clerk Maxwell Telescope (JCMT) GBS (Ward-Thompson et al. 2007) have provided unprecedented views of star-forming regions at resolutions where dense cores and filaments are resolved (≲0 farcm 5). Despite this advancement, only a few studies (e.g., Sadavoy et al. 2013, 2016; Schnee et al. 2014) have attempted to map out β at resolutions similar to these surveys over large areas. Such β measurements can be very valuable in providing more accurate mass and temperature estimates, and, consequently, gravitational stability estimates of the structures revealed in these surveys. While the Herschel GBS multi-band data may seem ideal for making β measurements in cold star-forming regions at first, Sadavoy et al. (2013) demonstrated that Herschel data alone are insufficient and longer wavelength data (e.g., JCMT 850 μm observations) are needed to provide sufficient constraints on β.

Here we present the first results of the JCMT GBS 850 μm observation with the Sub-millimetre Common User Bolometer Array 2 (SCUBA-2) toward the Perseus Molecular Cloud as a whole. Employing the technique developed by Sadavoy et al. (2013) for combining the Herschel and JCMT GBS data, we simultaneously derive maps of T_d, β, and optical depth at 300 μm, τ₃₀₀, in the Perseus star-forming clumps by fitting modified blackbody SEDs to the combined data. We investigate the robustness and the uncertainties of our SED fits, and perform a detailed analysis of the derived parameter maps in an attempt to understand the local environment. In particular, we characterize the observed variation in β and discuss the potential physical driver behind such evolution.

In this paper, we describe the details of our observation and data reduction in Section 2, and our SED-fitting method in Section 3. The results, including the parameter maps drawn from SED fits, are presented in Section 4. We discuss the implications of our results in Section 5, with regards to radiative heating (Section 5.1), outflow feedback (Section 5.2), and dust grain evolution (Section 5.3). We summarize our conclusions in Section 6.

2. OBSERVATIONS

2.1. JCMT: SCUBA-2 Data

Wide-band 850 μm observations of Perseus were taken with the SCUBA-2 instrument (Holland et al. 2013) on the JCMT as part of the JCMT GBS program (Ward-Thompson et al. 2007). We included observations that were taken in the SCUBA-2 science verification (S2SV) and the main SCUBA-2 campaign of the GBS program, i.e., in 2011 October, and between 2012 July and 2014 February, respectively. As with the rest of the survey, Perseus regions were individually mapped using a standard PONG1800 pattern (Kackley et al. 2010; Holland et al. 2013; Bintley et al. 2014) that covers a circular region ∼30' in diameter.

Our observations covered the brightest star-forming clumps found in Perseus, namely B5, IC 348, B1, NGC 1333, L1455, L1448, and L1451, listed here in order from east to west. Table 1 shows the names and center coordinates of the observed PONG1800 maps along with the weather grades in which they were observed. Based on priority, each planned PONG target was observed either four times under very dry conditions (Grade 1; τ₂₂₅ < 0.05) or six times under slightly less dry conditions (Grade 2; τ₂₂₅ = 0.05 − 0.07) to reach the targeted survey depth of 5.4 mJy beam⁻¹ for 850 μm. The "northern" PONG region of NGC 1333 is the only exception, containing one extra S2SV PONG900 map observed under Grade 2 weather. As with Sadavoy et al. (2013), we adopted 14 farcs 2 as our effective Gaussian FWHM beam width for the 850 μm data based on the two components of the SCUBA-2 beam obtained from measurements.

Table 1. Details of the Observed PONG Regions

Scan Name	R.A.	Decl.	Clump	Weather Grade^a	Number of Scans
B5	03:47:36.92	+32:52:16.5	B5	2	6
IC348-E	03:44:23.05	+32:01:56.1	IC 348	1	4
IC348-C	03:42:09.99	+31:51:32.5	IC 348	2	6
B1	03:33:10.75	+31:06:37.0	B1	1	4
NGC1333-N	03:29:06.47	+31:22:27.7	NGC 1333	1	4
NGC1333	03:28:59.18	+31:17:22.0	NGC 1333	2	1
NGC1333-S	03:28:39.67	+30:53:32.6	NGC 1333	2	6
L1455-S	03:27:59.43	+30:09:02.1	L1455	1	4
L1448-N	03:25:24.56	+30:41:41.5	L1448	1	4
L1448-S	03:25:21.48	+30:15:22.9	L1451	2	6

Note. The names, center coordinates, targeted clump names, and the weather grades of the individual observations made with the PONG1800 scan pattern. Clumps in the table are ordered from east to west.

^aThe weather grades 1 and 2 correspond to the sky opacity measured at 225 GHz of τ₂₂₅ < 0.05 and 0.05 ≤ τ₂₂₅ < 0.07, respectively.

Download table as: ASCII Typeset image

All SCUBA-2 data observed for the JCMT GBS program were reduced with the makemap task from the Starlink SMURF package (Version 1.5.0; Jenness et al. 2011; Chapin et al. 2013) following the JCMT GBS Internal Release 1 (IR1) recipe, consistent with several earlier JCMT first-look papers (e.g., Salji et al. 2015). The ¹²CO J = 3 − 2 line contribution to the 850 μm data is further removed with the JCMT Heterodyne Array Receiver Program (HARP) observations (see Buckle et al. 2009) as part of the IR1 reduction. Weather-dependent conversion factors calculated by Drabek et al. (2012) were applied to the HARP data prior to the CO subtraction. For more details on the IR1 reduction and CO subtraction, see Sadavoy et al. (2013) and Salji et al. (2015).

2.2. Herschel: PACS and SPIRE Data

The Perseus region was observed with the PACS (Photodetector Array Camera and Spectrometer; Poglitsch et al. 2010) instrument and the SPIRE instrument (Spectral and Photometric Imaging Receiver; Griffin et al. 2010; Swinyard et al. 2010) as part of the Herschel GBS program (André & Saraceno 2005; André et al. 2010; Sadavoy et al. 2012, 2014), simultaneously covering the 70, 160, 250, 350, and 500 μm wavelengths using the fast (60'' s⁻¹) parallel observing mode. The observation of the western and the eastern portions of Perseus took place in 2010 February and 2011 February, respectively, covering a total area of ∼10 deg². We reduced our data with Version 10.0 of the Herschel Interactive Processing Environment (HIPE; Ott 2010), using modified scripts written by M. Sauvage (PACS) and P. Panuzzo (SPIRE) and PACS Calibration Set v56 and the SPIRE Calibration Tree 10.1. Version 20 of the Scanamorphos routine (Roussel 2013) was used in addition to produce the final PACS maps. The final Herschel maps have resolutions³⁷ of 8 farcs 4, 13 farcs 5, 18 farcs 2, 24 farcs 9, and 36 farcs 3, in order of the shortest to longest wavelength. For more details on the Herschel observations of Perseus, see Pezzuto et al. (2012), Sadavoy (2013), and S. Pezzuto et al. (2016, in preparation).

The SCUBA-2 data are spatially filtered to remove slow, time-varying noise common to all bolometers, such as atmospheric emission to which ground-based sub-millimeter observations are susceptible. For effective combinations, we filtered the Herschel data using the SCUBA-2 mapmaker makemap (Jenness et al. 2011; Chapin et al. 2013) following the method described by Sadavoy et al. (2013) to match the spatial sensitivity of the Herschel data to that of the SCUBA-2 observations. Given that any zero-point flux offsets in the Herschel data are removed by the spatial filtering along with other large-scale emission, no offset correction was applied to our Herschel data. Further details on the parameters used for filtering the Herschel data are described in Chen (2015). The reduced SCUBA-2 850 μm and the filtered 160, 250, 350, and 500 μm maps of Perseus are available at https://doi.org/10.11570/16.0004.

3. SED FITTING

We modeled our dust SEDs as a modified blackbody function in the optically thin regime in the form of

$\begin{eqnarray}&&{I}_{\nu }={\tau }_{{\nu }_{0}}{(\nu /{\nu }_{0})}^{\beta }{B}_{\nu }({T}_{d}),\end{eqnarray} \tag{ 1 }$

where ${\tau }_{{\nu }_{0}}$ is the optical depth at frequency ν₀, β is the dust emissivity power-law index, and ${B}_{\nu }({T}_{d})$ is the blackbody function at the dust temperature T_d. By adopting a dust opacity value, ${\kappa }_{{\nu }_{0}}$ , one can further derive the gas column densities as the following using the ${\tau }_{{\nu }_{0}}$ values:

$\begin{eqnarray}&&N({H}_{2})=\displaystyle \frac{{\tau }_{\nu 0}}{\mu {m}_{{H}_{2}}{\kappa }_{\nu 0}},\end{eqnarray} \tag{ 2 }$

where μ is the mean molecular weight of the observed gas and ${m}_{{H}_{2}}$ is the mass of a molecular hydrogen in grams. For this study, we adopted a reference frequency ν₀ = 1 THz (300 μm), μ = 2.8, and ${\kappa }_{{\nu }_{0}}=0.1$ cm² g⁻¹, consistent with assumptions made by the Herschel GBS papers (e.g., André et al. 2010).

We convolved the CO-removed SCUBA-2 850 μm data and the spatially filtered Herschel data with Gaussian kernels to a common resolution of 36 farcs 3 to match that of the 500 μm Herschel map, the lowest resolution of our data. We re-aligned and re-gridded the convolved maps to the original 500 μm Herschel map, which has 14'' × 14'' pixels. The 70 μm data were excluded from our SED fitting because that emission may trace a population of very small dust grains that are not in thermal equilibrium with the dust traced by the longer wavelength emission (Martin et al. 2012).

Following Sadavoy et al. (2013), we applied color correction factors of 1.01, 1.02, 1.01, and 1.03 to the 160, 250, 350, and 500 μm Herschel data, respectively, to account for the spectral variation within each band. The correction factors are taken from the averaged values calculated from SED models with specific T_d and β values (i.e., 10 K ≤ T_d ≤ 25 K and 1.5 ≤ β ≤ 2.5). The color calibration uncertainties associated with these color correction factors are 0.05, 0.008, 0.01, and 0.02 times that of the color corrected values, respectively.

For each pixel of the map with signal-to-noise ratios (S/N) ≥ 10 in all five bands, we fitted a modified blackbody function (Equation (1)) to get best estimates of β, T_d, and ${\tau }_{{\nu }_{0}}$ . We employed the minimization of χ² method using the optimize.curve_fit routine from the Python SciPy software package, which uses the Levenberg–Marquardt algorithm for minimization. The flux uncertainties were calculated as the quadrature sum of the color calibration uncertainties and the map sensitivities (see Table 2), and were adopted as standard errors for the χ² calculation. For each clump, we adopted the 1σ rms noise measured from a relatively smooth, emission-free region within the convolved, filtered maps as the map sensitivity.

Table 2. Noise Levels in the Perseus Data

Band	160 μm	250 μm	350 μm	500 μm	850 μm
B5	60	50	30	20	20
IC 348	150	110	50	20	40
B1	80	90	60	30	30
NGC 1333	50	60	30	20	30
L1455	60	70	50	20	30
L1448	60	50	30	20	30
L1451	60	50	30	20	30

Note. The approximate 1σ rms noise levels (mJy beam⁻¹) of the convolved, spatially filtered maps at a resolution of 36 farcs 3 for different Perseus clumps. The rms noise values were measured in a relatively smooth, emission-free region of the filtered clump maps.

Download table as: ASCII Typeset image

In addition to color calibration uncertainties and map noise, we followed Könyves et al. (2015) in adopting 20% as the absolute flux calibration uncertainties for the PACS 160 μm band and 10% for the SPIRE bands. We also followed Sadavoy et al. (2013) in adopting 10% as the absolute flux calibration uncertainty in the SCUBA-2 850 μm band. We note that 10% is a conservative estimate of the absolute flux calibration uncertainties for extended sources and the calibration accuracy of point sources are much higher (<7% for the PACS bands, Balog et al. 2014; ∼5% for the SPIRE bands, Bendo et al. 2013; <8% for the SCUBA-2 850 μm band).

The flux calibration uncertainties are assumed to be correlated among the bands within an instrument, and we employed a simple Monte Carlo method to determine the most probable fit given the assumed calibration uncertainties based on 1000 instances. The resulting probability distribution for a given pixel was collapsed onto the T_d and β axes separately and fitted with a Gaussian distribution to determine the associated 1σ uncertainties of the derived T_d and β. Due to the asymmetric distribution of the derived τ₃₀₀, we used the best-fit τ₃₀₀ from another SED fit with T_d and β fixed at their mean best-fit values from their respective Gaussian fits. We derived the uncertainty in τ₃₀₀ from the covariance of this final fit. For more details on the SED fitting and how the flux calibration uncertainties are handled, see Sadavoy et al. (2013) and Chen (2015).

The β and T_d derived using the minimization of ${\chi }^{2}$ method can produce artificial anti-correlations due to the shape of the assumed SED becoming degenerate in the presence of moderate noise (Shetty et al. 2009b, 2009a) and flux calibration uncertainties. Detailed analyses of these correlated uncertainties are presented in the Appendix. To ensure the robustness of our results, we rejected pixels with β uncertainties >30%, where the fits to the SED are often poor and the β and T_d uncertainties occupy a distinct part of the parameter space relative to the main population. We also removed isolated regions that contained less than 4 pixels and where Gaussian fits to the derived β and T_d distributions from the Monte Carlo simulation were poor. The final T_d, β, and τ₃₀₀ maps of Perseus, along with their associated error maps, are available at https://doi.org/10.11570/16.0004.

4. RESULTS

4.1. Reduced Data

We present the reduced SCUBA-2 850 μm maps of the B5, IC 348, B1, NGC 1333, L1455, L1448, and L1451 clumps in Figure 1. For a comparison, we also present the reduced, unfiltered SPIRE 250 μm map in Figure 2, cropped to match the same regions shown in Figure 1. For the complete set of reduced, unfiltered PACS and SPIRE maps of Perseus, see S. Pezzuto et al. (2016, in preparation).

**Figure 2.** Unfiltered SPIRE 250 μm maps over the same regions shown in Figure 1. The fluxes are shown on a log scale from 0.05 Jy beam⁻¹ to 110 Jy beam⁻¹.
Download figure:
Standard image High-resolution image

4.2. Derived Dust Temperatures

Figure 3 shows the T_d maps for all seven Perseus clumps overlaid with the positions of embedded young stellar objects (YSOs), B stars, and A stars. The YSOs shown here are Class 0/I (circles) and Flat (squares) protostars identified from the Spitzer Gould Belt catalog of mid-infrared point sources (Dunham et al. 2015). The B and A stars identified in various catalogs (referenced in the SIMBAD database; Wenger et al. 2000) are labeled with stars and triangles, respectively. The contours of the unfiltered Herschel 500 μm emission at the 1.5 Jy beam⁻¹ are overlaid on the maps in gray. As demonstrated in the Appendix, the uncertainty in the T_d measurement due to flux calibration uncertainties is dependent on T_d. The typical uncertainties at T_d values of <12 K, 12–15 K, 15–20 K, and 20–30 K are 0.9 K, 1.5 K, 2.5 K, and 5.0 K, respectively. The derived T_d values at T_d > 30 K are poorly constrained and typically have uncertainties of ≳10 K. Our derived T_d values are systematically lower than those obtained from unfiltered Herschel data due to the preferential removal of the warm, diffuse dust emission by the spatial filtering. See Chen (2015) for further details on the bias introduced by the spatial filtering.

The T_d structures seen in Figure 3 often appear as circular, localized peaks overlaid on relatively cold backgrounds (∼10–11 K). Nearly all localized T_d peaks are coincident with at least an embedded YSO or a B star along their respective lines of sight, indicating that these objects are likely the source of local heating. The few T_d peaks that do not contain an embedded YSO, B star, or A star are all located in IC 348 and partially on the map edges. These T_d peaks may be the result of local heating by sources outside of the map, such as the nearby A star shown in Figure 3 or the nearby star cluster.

The region just southeast of IRAS 2 in NGC 1333 also appears relatively warm (∼14–15 K) with no clear source of localized heating. While the star BD +30 547 (also known as ASR 130, SSV 19) adjacent to IRAS 2B has often been classified as a foreground, late-type star (e.g., G2 IV, Cernis 1990), Aspin (2003) suggested that it may be a late-B V star with a cool faint companion, making BD +30 547 a potential source of heating in this region.

Interestingly, not all embedded YSOs are coincident with localized temperature peaks. This result is in agreement with that found by Hatchell et al. (2013) in NGC 1333. Out of the 61 embedded YSOs located within the temperature maps of Perseus clumps, only 7 are not coincident with local temperature peaks. Some of these YSOs may be too faint, embedded, or young to have warmed their surrounding dust significantly. Alternatively, they could be more evolved, less embedded YSOs misidentified as Class 0/I or Class Flat objects, potentially due to line of sight confusion.

Figure 4 shows histograms of derived dust temperatures, T_d, in each of the Perseus clumps. The light and dark gray histograms represent all the pixels within a clump and the pixels found within a 72 farcs 6 diameter (i.e., two Herschel 500 μm beam widths) area centered on a Class 0/I and Flat YSOs, respectively. Globally, the T_d distributions of Perseus clumps all share two characteristic features: a primary low T_d peak and a high T_d tail. L1451 is the only clump without a high T_d tail. All clumps have their primary T_d peaks located at ∼10.5 K, with the exception of IC 348, which has a temperature peak at ∼11.5 K. The former is consistent with the typical T_d seen in NGC 1333 (Hatchell et al. 2013), the kinetic temperatures observed toward dense cores in Perseus with ammonia lines (Rosolowsky et al. 2008; Schnee et al. 2009), and the isothermal T_d of prestellar cores in other clouds (Evans et al. 2001).

farcs — **Figure 4.** Histograms of derived T_d values in the seven Perseus clumps. All the pixels within a clump are shown in light gray while the pixels within a 726 diameter (i.e., two *Herschel* 500 μm beam widths) centered on a Class 0/I and Flat YSOs are shown in dark gray.
Download figure:
Standard image High-resolution image

As shown in Figure 4, nearly all the pixels found in the high T_d tails of the overall distribution (light gray) of B5, B1, L1455, and L1448 are located within a 72 farcs 6 diameter area centered on an embedded YSO (dark gray), indicating that these YSOs are indeed the main source of heating in these clumps. While protostellar heating also appears to be significant in IC 348 and NGC 1333, as seen in the T_d maps (Figure 3), only slightly less than half of the high T_d pixels in these two clumps are found near embedded YSOs. The remainder of the pixels in the high T_d tails of the distributions are likely from dust externally heated by B stars and nearby star clusters instead. No embedded YSO is found near pixels in L1451 where the SEDs were fitted. This lack of local heating sources likely explains why L1451 does not have a high T_d tail.

Looking from another prespective, nearly all the pixels surrounding embedded YSOs (dark gray) in B5, L1455, and L1448 are found in the high T_d tail of the overall distribution (light gray). This situation, however, is not seen in IC 348, B1, and NGC 1333, where a large fraction of the pixels near the embedded YSOs are found below ∼13 K, near the 10.5 K T_d peak. These low T_d pixels are associated with the embedded YSOs not found toward the locally heated regions, or, in some cases (e.g., B1), toward localized T_d peaks that appear relatively cool.

Our derived T_d values do not consistently drop off near map edges. This result indicates that our T_d derivation is less susceptible to the edge effects caused by the large-scale filtering seen in the T_d map of NGC 1333 derived by Hatchell et al. (2013) using the 450 and 850 μm data from the early SCUBA-2 shared risk observations. The northwestern edge of our NGC 1333 map was the only exception, where the edge of the fitted map lies near the edge of the external mask used for the data reduction and spatial filtering. Nevertheless, the regions in Perseus where SEDs were fit are typically located well within the adopted external masks and are therefore unlikely to have experienced such a problem. We improved T_d derivations in NGC 1333 with respect to Hatchell et al.'s analysis by using four additional bands from Herschel and allowing β to vary. Furthermore, the JCMT GBS data used here were taken with longer integration time with the full SCUBA-2 array instead of just one sub-array and had less spatial filtering. At T_d ∼ 10 K, Hatchell et al.'s T_d uncertainties due to flux calibration alone is ∼1.2 K, whereas our typical uncertainty is ∼0.9 K.

The appearance of our derived T_d map of B1 is morphologically consistent with that derived by Sadavoy et al. (2013) using the same observations reduced with earlier recipes. In other words, the structures in the two T_d maps are qualitatively the same in most places. Due to improvements in data reduction, our T_d values are more reliable than those derived by Sadavoy et al. and tend to be systematically colder by ≲1 K at T_d ≲ 10 K and warmer by ≲1 K at T_d ≳ 11 K. The overall T_d distribution of our map is broader near the 10.5 K T_d peak than that of Sadavoy et al., with a more pronounced higher T_d tail.

4.3. Derived β and τ₃₀₀

Figure 5 shows maps of derived β in the seven Perseus clumps. The 1σ uncertainties in the β measurement due to flux calibration uncertainties are relatively independent of the derived β values and have a median value of ∼0.35 in Perseus. Pixels with similar β values tend to form well defined structures, indicating that β variations seen in these clumps are correlated with local environment and not noise artifacts. Interestingly, low-β structures tend to correlate with local T_d peaks.

To a lesser extent, low-β structures also appear to correlate with outflows traced by CO emission (e.g., see Figure 10 for details). The CO contribution to these maps has been subtracted although the percentage contamination seen in most pixels is less than the flux uncertainties of 10%, with the exception of a few special, local cases (Chen 2015). The additional uncertainty on the 850 μm fluxes associated with the CO removal process is negligible. The resemblance between some of the β minima and outflow structures is therefore unlikely to be an artifact due to errors in the CO decontamination. While no study of free–free emission has been conducted over these outflows to assess whether or not such emission can be a significant contaminant in our data, free–free emission at centimeter wavelengths observed in radio jets is generally <1 mJy (Anglada 1996). Free–free emission is also expected to be relatively weak at 850 μm compared to the rms noise of our 850 μm data, given that these jets have widths much smaller than the JCMT beam. High angular resolution observations of the outflow sources SVS 13 (Rodríguez et al. 1997; Bachiller et al. 1998) and IRAS 4A (Choi et al. 2011) in NGC 1333 with the Very Large Array (VLA) have also shown that free–free emission is negligible at λ ≲ 3 mm at these locations. The VLA survey of the Perseus protostars conducted by Tobin et al. (2016) also found free–free emission near protostars to be faint relative to sub-millimeter dust emission and very compact with respect to our convolved 36 farcs 3 beam. Free–free emission is therefore unlikely to contaminate our data.

Pixels with β ≥ 3 are found only in NGC 1333, mostly on the northwestern edge where the filtering mask may be affecting the emission, as discussed in Section 4.2. Therefore, we do not trust these steep β values. Low-β values (≲1.5), however, are generally well within the filtering mask, and are likely robust against the filtering systematics.

Figure 6 shows the histograms of derived β values from various clumps in Perseus. Unlike T_d, the β distributions of different clumps do not share a similar shape. While four clumps in Perseus appear to have a single peak (i.e., NGC 1333, IC 348, L1455, and B5), the other three clumps (i.e., B1, L1448, and L1451) appear to have two or even three peaks. The separations between most of these peaks are larger than the median β uncertainty (∼0.35), suggesting that these multiple peaks are not artifacts. The β values of the primary peak in each clump range between 1.5 (e.g., L1448) and 2.2 (e.g., B1), and the β values of most pixels range between 1.0 and 2.7. The latter range of β values is similar to that found in nearby star-forming clouds by Dupac et al. (2003; 1.0 ≤ β ≤ 2.5) and for luminous infrared galaxies by (Yang & Phillips 2007; 0.9 ≤β ≤ 2.4).

The appearance of our derived β map of B1 is morphologically consistent with that derived by Sadavoy et al. (2013) using the same observations reduced with earlier recipes. Due to improvements in the data reduction, our derived β values are more reliable. The majority of our derived β values are lower than those derived by Sadavoy et al. by ≲0.7, with the remaining pixels having higher values by ≲0.3. The overall distribution of our derived β is shifted downwards by ∼0.2 relative to that of Sadavoy et al., with the primary peak being skewed toward the higher end of the distribution instead of the lower end.

Figure 7 shows maps of derived τ₃₀₀ in the seven Perseus clumps, overlaid with the same symbols as Figures 3 and 5. Column densities can be further derived from τ₃₀₀ using Equation (2) by assuming a κ_ν0 value. The median uncertainty in the τ₃₀₀ measurement derived from the covariance of the fit with T_d and β fixed at the best determined values is ∼1.5%. The typical uncertainties due to flux calibration uncertainties, however, as seen in the Monte Carlo simulation, is about a factor of 2 (see the Appendix). Higher τ₃₀₀ structures seen in the B1, IC 348, L1448, and NGC 1333 clumps are found along filamentary structures, much like their dust emission counterparts. This similarity is less clear in the B5, L1455, and L1451 clumps where the areas where SEDs were fit (S/N > 10) are relatively small.

Embedded YSOs are preferentially found toward local τ₃₀₀ peaks. Many of these embedded YSOs, however, are spatially offset from the center of these peaks, often by slightly more than a beam width of our maps. The differences between β values found between these offsets are smaller than those expected from the anti-correlated uncertainties discussed in the Appendix. Interestingly, we did not find any embedded YSOs toward the centers of the highest τ₃₀₀ peaks in all Perseus clumps, with the exception of IC 348, which does not have τ₃₀₀ peaks with distinct high values. The T_d values of these starless τ₃₀₀ peaks are also much lower than their surroundings, suggesting that these structures are cores well shielded from the interstellar radiation field (ISRF) due to their high column densities.

The T_d values near embedded YSOs and toward higher ${\tau }_{300}$ regions also tend to be lower than in the lower τ₃₀₀ regions, suggesting that protostellar heated regions can appear cooler due to being more deeply embedded and having more cool dust along the line of sight. Table 3 shows the mean column densities found in each of the Perseus clumps. Indeed, the clumps with significant amounts of cold (≲13 K) pixels near embedded YSOs, i.e., B1, NGC 1333, and IC 348, have the first, third, and fourth highest mean column densities in Perseus, respectively. The B1 clump in particular has a majority of its pixels near embedded YSOs at ≲13 K.

Table 3. Derived Column Densities in Perseus Clumps

Clump	$\bar{N}({H}_{2})$	${\sigma }_{N({H}_{2})}$	Number of Pixels
B5	1.0	0.54	120
IC 348	1.4	1.0	537
B1	2.4	1.8	616
NGC 1333	1.9	1.7	1362
L1455	1.2	0.59	119
L1448	2.0	1.3	449
L1451	0.78	0.31	158

Note. The units of mean and standard deviation of column densities, $N({H}_{2})$ , are both 10²² cm⁻².

Download table as: ASCII Typeset image

4.4. The β, T_d, and τ₃₀₀ Relations

Figure 8 shows scatter plots of β versus T_d for all seven Perseus clumps. At T_d ≲ 16 K, anti-correlations between β and T_d are found in all cases and cannot solely be accounted for by the anti-correlated β —T_d uncertainties discussed in the Appendix. At these T_d values, there appears to be a prominent population of pixels in all seven clumps that exhibits a fairly linear relationship with a slope of ∼–0.3. Three clumps in particular, B1, B5, and L1451, seem to consist only of this population. While this slope is similar to those found in the anti-correlated uncertainties presented in the Appendix (see Figure 14), the range of β found in this population is greater than the calculated β uncertainties.

The other four Perseus clumps contain additional populations that extend into warmer temperatures (T_d ≳ 16 K) and have slopes which are much shallower than −0.3. In particular, the second population found in L1448 at T_d ≳ 13.5 K appears to have a constant β value of ∼1.4 with a spread of ±∼0.15. This spread in β is smaller than the estimated uncertainties based on the assumed absolute flux calibration error, which is systematic across a map. This result suggests that our relative, pixel to pixel uncertainties within a map are indeed relatively small, and further indicates that there indeed is an underlying anti-correlation between β and T_d. The β values found in these warmer populations are generally lower than those found at colder temperatures (T_d ∼ 10 K), but are not necessarily the lowest values found within a clump.

Figure 9 shows scatter plots of derived β versus τ₃₀₀, where these appear to be weakly positively correlated. Given the positive correlation between β and τ₃₀₀ uncertainties discussed in the Appendix, however, we do not consider this correlation to be significant.

5. DISCUSSION

5.1. Radiative Heating

Throughout Perseus, we see evidence of embedded YSOs and nearby B stars heating their local environments. In contrast with the most common T_d found in Perseus, ∼10.5 K, the T_d near embedded YSOs are typically ∼14–20 K and can be well above 30 K near a B star where T_d is not well constrained with our data. While not all embedded YSOs are associated with local T_d peaks, all the local T_d peaks coincide with locations of embedded YSOs and B stars. The only exceptions are the two extended warm regions in IC 348 that appear to be heated by sources outside of our parameter maps. Indeed, the main star cluster in IC 348 is located just northeast of the eastern warm region, while an A star is located just northeast of the western warm region. A higher local ISRF from the nearby young star cluster (Tibbs et al. 2011) likely explains why the most common T_d found in IC 348 (∼11.5 K) is slightly warmer than the other Perseus clumps.

The majority of the heated regions centered on embedded YSOs are larger than the FWHM beam of our map (36 farcs 6), often having ${T}_{d}\geqslant 15$ K out to one full beam width from the center. In NGC 1333, all the locally heated regions with radii nearly one beam width in extent coincide with at least a Class 0 or I YSO with bolometric luminosity, L_bol, of 4 L_⊙–33 L_⊙ (Enoch et al. 2009). In other Perseus clumps, similar regions coincide with at least a Class 0 or I YSO with bolometric luminosity of 1 L_⊙–5 L_⊙. If we assume our distances to the western and eastern sides of Perseus to be 220 pc (Cernis 1990; Hirota et al. 2008) and 320 pc (Herbig 1998), respectively, then these heated regions would have radii of ∼8000 AU and ∼12000 AU, respectively, which is somewhat less than the typical size of a dense core (∼15,000 AU in radius, André et al. 2000, p. 59).

The core T_d profile due to protostellar heating can be modeled with radiative transfer codes. Using Dusty (Ivezic et al. 1999), Jørgensen et al. (2006) modeled such T_d profiles for various core density profiles, protostellar bolometric luminosities, and ISRF, assuming the OH5 dust opacities (Ossenkopf & Henning 1994). For a core modeled after the MMS9 core in the Orion molecular cloud with a density profile of $n={n}_{0}{(r/{r}_{0})}^{-1.5}$ , where r₀ = 50 AU and n₀ = 2.7 × 10⁸ cm⁻³, its T_d drops down to 15 K at a radius of 2000 AU when an ISRF based on the solar neighborhood value (i.e., the standard ISRF) and a L_bol = 10 L_⊙ for the central source are assumed. Despite having assumed an L_bol that is relatively high, the size of this local heating is a factor of 4–6 smaller than what we find in Perseus. Given that the angular size of this modeled T_d profile is smaller than our beam, further analysis that takes beam convolution and its associated errors into account will be needed to make a definitive comparison between the observed T_d profile and the model.

Regions locally heated by B stars are more extended than regions heated by protostars. In NGC 1333, we find regions heated by B stars to be ∼1 farcm 6 in radius, corresponding to ∼ 21000 AU. Since these B stars are not embedded objects, the temperature profile of their surrounding dust can be reasonably approximated in the optically thin limit. An analytic expression of such a temperature profile first derived by Scoville & Kwan (1976) can be generalized as

$\begin{eqnarray}&&{T}_{d}(r)=50{Q}_{\mathrm{abs}}^{-q/2}{\left(\displaystyle \frac{r}{2\times {10}^{15}{\rm{m}}}\right)}^{-q}\left(\displaystyle \frac{{L}_{\mathrm{bol}}}{{10}^{5}{{L}}_{{\odot }}}\right){\rm{K}},\end{eqnarray} \tag{ 3 }$

where Q_abs is the absorption efficient of a dust grain at λ = 50 μm, and $q=2/(4+\beta )$ (Equation (1), Chandler et al. 1998). For the expected range of β values between 1 and 2, the q parameter does not vary significantly. The power-law slopes of the temperature profile at these β values resemble the Dusty profile of Jørgensen et al. (2006) at larger radii where the protostellar envelope is expected to be optically thin, before the ISRF heating starts to become dominant.

We modeled the radial temperature profiles near the two B stars in NGC 1333 by assuming β = 1.6, the typical β found near these two B stars, and normalizing Equation (3) to the Wolfire & Churchwell (1994) dust emission models with ${Q}_{\mathrm{abs}}^{-q/2}=1.2$ , as Chandler et al. (1998), Chandler & Richer (2000), and Jørgensen et al. (2006) did. Given that the B stars BD +30 549 and SVS 3 have luminosities of 42 L_⊙ (Jennings et al. 1987) and 138 L_⊙ (Connelley et al. 2008), respectively, the temperatures of these two stars are expected to drop down to 15 K at radii of ∼13000 AU and 24000 AU, comparable to the projected distances of 21000 AU (i.e., ∼1 farcm 6) that we observed.

The localized radiative heating due to YSOs and B stars may provide star-forming structures with additional pressure support. When treating a star-forming structure as a body of isothermal ideal gas with uniform density, the smallest length at which such an object can become gravitationally unstable due to perturbation, i.e., the Jeans length (λ_J), is related to the object's temperature by λ_J ∝ T^1/2. The total mass enclosed within a Jeans length, i.e., the Jeans mass (M_J), scales as ${M}_{{\rm{J}}}\propto {T}^{3/2}$ .

Even though protostellar heating is the most common form of local heating observed in Perseus, protostellar cores are likely no longer Jeans stable given that they are currently collapsing into protostars. Protostellar heating could also provide additional Jeans stability to inter-core gas, but none of the warm regions in Perseus heated by solar-mass YSOs extends beyond the typical radius of a core (≲15,000 AU; André et al. 2000, p. 59).

While B stars are unlikely to have circumstellar envelopes themselves, they may be capable of heating nearby protostellar cores. Indeed, three SCUBA cores identified by Hatchell et al. (2005) are located in projection within these B-star heated regions: HRF57, HRF56, and HRF54. These cores do not contain known YSOs and have beam averaged T_d between 19 K and 40 K. If we assume the densities of these cores remain unchanged in the process of heating and the gas in these cores are thermally coupled to the dust, then the associated Jeans masses of these cores at their respective temperatures are a factor of ∼3–8 higher than those of their counterparts at 10 K.

5.2. Outflow Feedback

Several studies of dust continuum emission in NGC 1333 (Lefloch et al. 1998; Sandell & Knee 2001; Quillen et al. 2005), as well as other nearby star-forming clumps (e.g., Moriarty-Schieven et al. 2006), have suggested that outflows play a significant role in shaping the local structure of star-forming clumps. Indeed, kinematic studies of outflows in Perseus (e.g., Arce et al. 2010; Plunkett et al. 2013) showed that outflows can inject a significant amount of momentum and kinetic energy into their immediate surroundings. As with these prior dust continuum studies, we find several τ₃₀₀ cavities that coincide with bipolar outflows, particularly toward the ends of outflow lobes, in NGC 1333 and L1448 (see Figures 10). Due to warmer T_d values, which can cause low column density dust to emit brightly, some of these cavities were not revealed directly by single-wavelength continuum emission alone.

**Figure 10.** Maps of derived β in (a) NGC 1333 and (b) L1448 and their τ₃₀₀ counterparts (c) and (d), respectively. These maps are overlaid with the SCUBA cores identified by Hatchell et al. (2005; circles), Class 0/I YSOs identified by Dunham et al. (2015; stars), and contours of spatially filtered, integrated ¹²CO 3–2 emission (gray). Contours of CO emission are drawn at 15 mJy beam⁻¹, 30 mJy beam⁻¹, 80 mJy beam⁻¹, and 110 mJy beam⁻¹ for the NGC 1333 maps and at 10 mJy beam⁻¹, 20 mJy beam⁻¹, and 25 mJy beam⁻¹ for the L1448 maps. The colors in panels (a) and (b) are scaled linearly between 1.0 and 3.0, and logarithmically between 0.002–0.2 and 0.002–0.1 for panels (c) and (d), respectively. Only the cores mentioned in the text are identified with numbers designated by Hatchell et al.
Download figure:
Standard image High-resolution image

In addition to low-τ₃₀₀ regions, we also find two high-τ₃₀₀ regions in Perseus that may have been formed from compression by nearby outflows. For example, the high-τ₃₀₀ region near the starless core HRF51 in NGC 1333 has the highest τ₃₀₀ values in the clump, and is also located toward the southeastern end of the outflow associated with the HH7-11 outflows driven by SVS 13 (Snell & Edwards 1981). The high-τ₃₀₀ region surrounding the HRF27 and HRF28 cores in L1448 also has the second highest τ₃₀₀ values found in L1448, and is located on the northeastern end of the outflow originating from L1448C (HRF29, also known as L1448-mm) near HH197 (e.g., Bachiller et al. 1990; Bally et al. 1997). The presence of Herbig-Haro objects in these high-τ₃₀₀ regions indicates that outflows are indeed interacting with ambient gas at these locations and perhaps compressing the gas (and dust) in the process, resulting in the observed high-τ₃₀₀ values. Given that the τ₃₀₀ values we find around the starless core HRF51 are unusually high in comparison with the rest of Perseus, HRF51 would be a good candidate for follow up observations to test the limits of starless core stability.

We do not find evidence for shock-heating from outflows in our T_d maps, which is unsurprising given that shock-heated gas is not well coupled to dust thermally (e.g., Draine et al. 1983; Hollenbach & McKee 1989). The dominant local influence of outflows appears to be their ability to reshape local structures such as cores and filaments. Indeed, gas compression driven by outflows can stimulate turbulence to form starless cores beyond simple Jeans instability. For example, the extra material being pushed toward a starless core by outflows can provide the needed perturbation to trigger core collapse, as well as providing the core with extra mass to accrete. On the other hand, these outflows can, in principle, also inject turbulence into the structures they compress, providing additional pressure support for these dense structures to remain temporarily stable against Jeans instability. The mechanical feedback from outflows therefore is likely important in regulating star formation in clustered environments, where it can enhance or impede the processes through compression or clearing, respectively.

5.3. Beta and Clump Evolution

The global properties of a star-forming clump are expected to evolve significantly as star formation progresses within it. Given that dust grains can undergo growth inside cold, dense environments and potentially be fragmented or destroyed in energetic processes, dust grains themselves are expected to evolve over time. If grain evolution in a clump is dominated by processes that result in a net change in β, such as growth in maximum grain sizes (Testi et al. 2014, p. 339), then we should be able to observe the global values of β varying from clump to clump due to the different evolutionary stages of these clumps.

5.3.1. Beta Variation and its Relation to Temperature

In Perseus, we find significant β variations within individual clumps. As shown in Figure 6, the majority of β values in these clumps range from 1.0 to 2.7. The spatial distribution of β is not erratic, and appears fairly structured on a local scale, suggesting that β is likely dependent on the properties of the local environment. Furthermore, the spatial transition between pixels with β > 2 and <1.6, though relatively sharp, is not abrupt with respect to our beam width (36 farcs 3 FWHM). The overall appearance of the β structures, as well as their T_d and τ₃₀₀ counterparts, is fairly smooth and coherent on scales larger than the beam, suggesting that our derived β values are robust against random noise.

In Figure 8, our derived β and T_d values appear to be anti-correlated, similar to behavior found in prior studies of nearby star-forming regions at lower angular resolutions (e.g., Dupac et al. 2003; Planck Collaboration et al. 2011). As detailed in the Appendix, this anti-correlation cannot be solely accounted for by the anti-correlated β—T_d uncertainties associated with our SED fits. The fact that lower β regions tend to be spatially well structured and coincide with locations that are locally heated by known stellar sources further demonstrates that there indeed exists an anti-correlation between β and T_d that is not an artifact of the noise.

Shetty et al. (2009b) have shown that line of sight T_d variations can lead to an underestimation of β and an overestimation of density-weighted T_d (column temperature) by fitting a single-component modified blackbody curve to the SEDs of various models of cold cores with warm envelopes that have a single β value. When these synthetic SEDs are noiseless and are equally well sampled in both the Wien's and Rayleigh–Jeans' regimes, Shetty et al. found that systematic exclusion of shorter wavelength data, i.e., the bands toward the Wien's regime, will allow the β and T_d values drawn from fits to approach the true β values and column temperatures. A single-component modified blackbody curve simply does not well fit SEDs with multiple temperature components in the Wien's regime.

In Perseus, we find the median T_d of starless cores to be 9.8 K when cores in IC 348 are excluded. The median T_d of starless cores in IC 348 is a bit higher (11.7 K), likely due the clump having a higher ISRF. The former T_d is consistent with both the typical core T_d found by Evans et al. (2001; ∼10 K) and the modeled column temperature of 9.6 K from Shetty et al. (2009b), which was itself based on the core density and T_d profiles of Evans et al. The agreement between our derived T_d and Shetty et al.'s modeled column temperature suggests that our fits to the SEDs of these cold cores are robust against line of sight T_d variation. If the true β values along the lines of sight of our observations are fairly uniform, then our derived β values should approach those of the true β, as demonstrated by Shetty et al. Figure 11(a) shows a typical example of an SED found in cold core that is well fit by a single-component modified blackbody curve.

**Figure 11.** Best SED fits to four selected pixels in NGC 1333 with T_d, β, and τ₃₀₀ as free parameters (solid), and with β fixed at a value of 2 (dash). The best-fit T_d values in panels (a)–(d) with β as a free parameter are 10.6 K, 11.4 K, 16.3 K, and 16.8 K, respectively. Their corresponding best-fit β values are 2.0, 2.4, 1.6, and 1.4, respectively. The error bars represent the flux calibration uncertainties in each bands (i.e., 20% for the SPIRE 160 μm band and 10% for the other bands).
Download figure:
Standard image High-resolution image

In warmer regions, some of our SED fits have trouble reconciling the shortest wavelength band (i.e., 160 μm), suggesting that line of sight T_d variation may be a concern in these specific regions. Figures 11(b)–(d) show some examples. Given that preferential sampling toward the Rayleigh–Jeans tail will allow SED fits to place a tighter constraint on β, as Shetty et al. demonstrated, the β values we measured across 160–850 μm toward the warmer regions should remain fairly robust against T_d variation along the line of sight as the peaks of these SEDs shift further toward shorter wavelengths. Indeed, the fact that our sampled SED fits only have trouble fitting the 160 μm bands in warmer regions suggests that these SED fits are preferentially anchored toward the longer wavelength points where β can be more robustly determined. Measurements of the flux ratio between the longer wavelength bands in NGC 1333 by Hatchell et al. (2013; SCUBA-2 450 and 850 μm bands) and in Perseus by Chen (2015; SPIRE 500 μm and SCUBA-2 850 μm bands) have further confirmed that these regions are indeed warm and the T_d drawn from the SED fits are not erroneous due to poor fits to the shorter wavelength bands.

5.3.2. The Cause Behind β Variations

The β value of a grain population depends on various physical properties of its grains. The β value we derived in a given pixel should therefore be viewed as an effective β of the various underlying populations seen along the line of sight. If these effective β values are dominated by the β of a single grain population, then a simple comparison with laboratory measurements may reveal some physical properties of this particular population.

Laboratory measurements have shown that β values can be intrinsically temperature-dependent (e.g., Agladze et al. 1996; Mennella et al. 1998; Boudet et al. 2005). Many of these measured dependencies, however, are not significant enough to explain by themselves the anti-correlation seen here. Agladze et al. (1996) found two species of amorphous grains with β values that decrease from ∼2.5 to ∼1.7–2 as the T_d increases from 10 to 25 K. While these measured temperature dependencies are not strong enough to be the sole cause of the anti-correlation seen here, they are perhaps significant enough to provide a partial contribution.

Dust coagulation models have shown that grain growth can cause β values to decrease significantly from the typical diffuse ISM value (β ∼ 2) to as low as zero (e.g., Miyake & Nakagawa 1993; Henning et al. 1995). Since dust coagulation is favored in a denser environment, where the dust collision rates are higher, lower β values may be expected toward the densest regions.

As shown in Figure 9, we do not find evidence that β is anti-correlated with τ₃₀₀ and thus column density. In contrast, we find that the highest τ₃₀₀ structures in Perseus tend to have β ≳ 2, especially toward the colder regions not associated with localized heating. The lack of decrease in β toward these dense regions, however, is not sufficient to rule out grain growth through coagulation. The presence of ice mantles on the surface of dust grains, for example, can increase the β value from that of bare grains, i.e., ∼2, to ∼3 for dust typically found in diffuse ISM (Aannestad 1975) and from ∼1 to ∼2 for modeled dust in the densest star-forming environments (Ossenkopf & Henning 1994). Since the growth of ice mantles is also favored in the coldest and densest environments, such growth can, in principle, compete with coagulated grain growth by modifying β in opposite directions, causing the dust in some of the highest τ₃₀₀ regions to have β ∼ 2 instead of higher or lower values.

Given that the presence of ice mantles on grain surfaces may be capable of suppressing the decrease in β caused by coagulated grain growth, we speculate that the β − T_d anti-correlation may be explained by the sublimation of ice mantles in warmer environments. This hypothesis is useful in explaining why internally heated protostellar cores have lower β values than their starless, non-heated counterparts. Measuring the depletion levels of molecules of which ice mantles typically composed of toward heated and non-heated cores could test this hypothesis, provided these molecules can themselves survive in these heated environments.

While lower β regions tend to coincide with locally heated regions, these β regions also tend to be more spatially extended than the heated regions. The low-β region surrounding IRAS 4³⁸ in NGC 1333 is a prominent example of this behavior. Interestingly, many low-β regions also coincide with the locations of bipolar outflows, prompting us to suspect that outflows may be capable of carrying lower β grains from deep within a protostellar core out to a scale large enough to be observed with our spatial resolution. Recent high-resolution interferometric studies on the inner regions of a small number of protostellar envelopes have found that β ≤ 1 (Chiang et al. 2012; Miotello et al. 2014), suggesting that dense, inner regions of cores may indeed be a site for significant grain growth. Further high angular resolution observations of β within protostellar and prestellar cores may verify if low-β grains are preferentially found deep within a dense core, and whether or not their production is associated with the presence of protostars or disks.

Interestingly, the dependence of β on the maximum grain size, a_max, assuming a power-law grain size distribution, i.e., $n(a)\propto {a}^{-q}$ , behaves similarly regardless of other β dependencies (Testi et al. 2014, p. 339; see their Figure 4). When a_max is less than ∼10 μm, the β values computed by Testi et al. are all between ∼1.6 and 1.7, independent of the dust grain's chemical composition and porosity. As a_max increases, β values for all three models first increase upwards to a maximum value of ∼2.5 for compact grains, or ∼2 for porous icy grains, and then quickly drop down to values less than 1 starting at a_max ≳ 10³ μm. The β values computed by Testi et al. are for wavelengths measured between 880 and 9000 μm, which do not overlap with our observed wavelengths between 160 and 850 μm. Nevertheless, given that a dust grain has less accessible modes to emit photons with wavelengths longer than or comparable to its physical size, we expect the same β dependency on a_max to hold in our observed wavelength range with the corresponding a_max shifted toward smaller values. Complementary β measurements at wavelengths greater than 850 μm will test this idea.

If we assume the β dependency on a_max to be the dominant, global driver of β evolution in a clump, as shown by Testi et al. (2014, p. 339), while the other β dependencies are only locally significant, then β values observed on the global scale of a clump should be a good proxy to measure the relative evolutionary stage of a clump in terms grain growth.

5.3.3. Beta Variations Between Clumps

As demonstrated in Figure 6, the β distributions in Perseus differ significantly from clump to clump, suggesting that the dust grains are co-evolving with their host clumps. To make the measurement of β variation more sensitive, we binned our β values into three categories, with β values of 1.5–2, 2–3, and 0–1.5 to represent pristine, transitional, and evolved values of a_max, and thus β, respectively. This choice is based on the Testi et al. (2014, p. 339) modeled β − a_max relations, and we extended the upper limit of the pristine, i.e., medium, β bin to a value of 2 to include the observed β values of ∼2 in the diffuse ISM. While β values of dust do go through the "pristine" category again as they migrate from the transition category to the evolved category, this typically only occurs over a very narrow range of a_max values relative to the overall a_max values considered. To quantify the advancement of grain growth in each clump, we devise the growth index, G, as follows:

$\begin{eqnarray}&&G=100-{P}_{\mathrm{pri}}+{P}_{\mathrm{evo},}\end{eqnarray} \tag{ 4 }$

where P_pri and P_evo are the percentage of pixels in each clump with the pristine and evolved values, respectively. As the typical a_max in a clump increases, the $100-{P}_{\mathrm{pri}}$ term increases as β migrates away from the pristine values, resulting in an increase in G. When the grain growth enters an advanced stage, the increase in a_max will also increase the value of P_evo, resulting in further increase in G.

Figure 12 shows the percentage of pixels in each of the β categories for the seven Perseus clumps and their corresponding G values. The β bins are arranged in the order of pristine, transitional, and evolved categories from left to right. The order of the clumps here has been rearranged in the order of increasing G, corresponding to the perceived increase in a_max based on calculations by Testi et al. (2014, p. 339). As seen in Figure 12, the dust population in each of the Perseus clumps appears to be in a different evolutionary stage. For a reference, the median β uncertainty is ∼0.35. In terms of their relative dust evolution, the Perseus clumps are ordered as L1455, L1451, IC 348, B5, B1, NGC 1333, and L1448. We note that our sample sizes vary significantly from clump to clump, ranging from 1362 pixels in NGC 1333 to only 119 pixels in L1455. Our ability to make comparison between clumps is therefore limited by these sample size differences and the associated uncertainties in measuring β.

It may seem surprising at first that the measured β values would suggest IC 348 to be a less evolved clump in Perseus, given that IC 348 has often been considered to be an older star-forming clump near the end of its star-forming phase (Bally et al. 2008). The region covered by our SED fits, however, actually traces a filamentary ridge of higher density gas just southwest of the main optical cluster. This sub-region hosts a higher fraction of Class 0 and I YSOs and a lower fraction of Class II and III YSOs relative to the rest of IC 348 (Herbst 2008). The relatively higher abundance of embedded YSOs in comparison to their more evolved counterparts in this area of IC 348 suggests that this sub-region is indeed quite young. Due to low number statistics (≤10 YSOs), we did not perform a similar consistency check on the youngest clumps, L1455 and L1451.

L1448, NGC 1333, and B1 are likely the more evolved clumps in Perseus based on their G values. Indeed, the higher fractions of Class II and III YSOs relative to Class 0, I, and Flat YSOs seen in NGC 1333 (Sadavoy 2013) suggest that NGC 1333 is one of the most evolved clumps in Perseus. Interestingly, B1 has the highest fraction of transitional β values and the lowest fraction of evolved β values in Perseus. Given that B1 also has the highest mean column density in Perseus (see Table 3), B1 may potentially be a younger clump with an enhanced grain growth due to its higher density.

In general, we do not expect the fraction of embedded YSOs and the relative dust evolutionary state to be tightly correlated given that the fraction of embedded YSOs is more indicative of the current star-forming rate, whereas the dust evolutionary state measures the accumulated grain growth throughout a clump's history. Considering that grain growth is subjected to factors such as density and temperature over time, which may also differ between clumps, we do not expect the relative evolution we inferred from a_max to be strictly proportional to time either.

As mentioned in Section 4.4, we find a significant anti-correlation between β and T_d in all seven Perseus clumps that cannot be solely accounted for by the anti-correlated uncertainties driven by noise or T_d variation along the line of sight. Interestingly, there appears to be a population of pixels in all seven clumps that are located at the bottom left corner of the β − T_d scatter plots with T_d less than 20 K (see Figure 8). The anti-correlation displayed by this particular population appears fairly linear, and has a slope of ∼−0.3. This population is primarily responsible for the wide range of β values observed in each clump and its presence in all seven clumps may suggest a universal dust evolution that is common to all star-forming clumps. The β − τ₃₀₀ distributions of Perseus clumps (see Figure 9) all behave fairly similarly and do not show obvious trends associated with clump evolution.

6. CONCLUSION

In this study, we fit modified blackbody SEDs to combined Herschel and JCMT continuum observations of star-forming clumps in the Perseus molecular cloud and simultaneously derived dust temperature, T_d, spectral emissivity index, β, and optical depth at 300 μm, τ₃₀₀, using the technique developed by Sadavoy et al. (2013). We performed a detailed analysis of the derived T_d, β, and τ₃₀₀ values in seven Perseus clumps, i.e., B5, IC 348, B1, NGC 1333, L1455, L1448, and L1451, and investigated anti-correlated uncertainties between β and T_d. We found that the β − T_d anti-correlations seen in Perseus are significant and are not artifacts of noise or calibration errors. In addition, we demonstrated that our final SED fits are robust against systematic uncertainties associated with temperature variation along lines of sight.

We summarize our main findings as follows.

1.
The most common (i.e., mode) T_d seen in Perseus clumps is ∼10.5 K, except in IC 348 where it is ∼11.5 K. The former T_d is consistent with the kinetic temperatures of dense cores observed with ammonia emission lines in Perseus (∼11 K; Rosolowsky et al. 2008; Schnee et al. 2009), and the isothermal T_d of prestellar cores seen in various clouds (∼10 K; Evans et al. 2001). The IC 348 clump being slightly warmer than the rest of Perseus may be due to stellar heating from the nearby young star cluster.
2.
Nearly all the local T_d peaks found in Perseus coincide with locations of embedded YSOs (and nearby B stars in the case of NGC 1333), which indicates local heating by these YSOs. The immediate region surrounding an embedded YSO can have T_d of up to ∼20 K, whereas that near a B star can have T_d well above ∼30 K.
3.
We found significant β variations over individual star-forming clumps. Most β values found in Perseus are in the range of 1.0 ≲ β ≲ 2.7, similar to that found in nearby star-forming clouds by Dupac et al. (2003; 1.0 ≤ β ≤ 2.5) and for luminous infrared galaxies by Yang & Phillips (2007; 0.9 ≤ β ≤ 2.4). Maps of β, in general, appear well structured on a local scale, indicating that β is likely a function of its local environment.
4.
We found β and T_d to be anti-correlated in all Perseus clumps, and that lower β regions tend to coincide with local T_d peaks. These anti-correlations cannot be solely accounted for by anti-correlated β and T_d uncertainties associated with our SED fitting. The dust grains' intrinsic β dependency on temperature (e.g., Agladze et al. 1996) may be partially responsible for the anti-correlation, but is not strong enough to be the sole cause. The sublimation of surface ice mantles, which can increase β when present on a dust grain (e.g., Aannestad 1975; Ossenkopf & Henning 1994), may provide an explanation for the observed anti-correlation.
5.
Similar to prior studies (e.g., Sandell & Knee 2001), we found cavities and peaks in maps of derived τ₃₀₀, and thus column density, that coincide with the ends of some outflow lobes, suggesting that these structures were cleared out or compressed by local outflows, respectively. The fact that Herbig-Haro objects were found in locations where outflows meet the highest column density structures in individual clumps further supports this view. Creating column density maps from fitting SEDs with T_d, β, and τ is important for this analysis.
6.
The β distribution found in Perseus differs from one clump to another (see Figure 6). Binning β values into pristine (medium), transitional (high), and evolved (low) values based on their corresponding maximum grain size, a_max, calculated by Testi et al. (2014, p. 339) revealed that the typical a_max value may also differ from clump to clump (see Figure 12). This result suggests that dust grains can grow significantly as a clump itself evolves.

Following Sadavoy et al. (2013), we mapped β values over several star-forming clump at angular resolutions of ∼0.5 arcminutes. Our study is the first of its kind to extend such an analysis to the major star-forming clumps across a molecular cloud. In Perseus, we found β values which are significantly different from the diffuse ISM values (i.e., ∼2), indicative of dust grain evolution in the denser, star-forming ISM. The discovery of coherent small-scale, low-β regions, and their coincidence with local temperature peaks have provided some hints on the origins of these low-β value grains. Further high-resolution observations of β toward protostellar cores and starless cores would be most welcome. Observations of potential sublimation of ice mantles may also help to constrain the role of ice mantles.

This work was possible with funding from the Natural Sciences and Engineering Research Council of Canada (NSERC) Postgraduate Scholarships. We acknowledge the support by NSERC via Discovery grants, and the National Research Council of Canada (NRC). We wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

The James Clerk Maxwell Telescope has historically been operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the United Kingdom, the NRC, and the Netherlands Organization for Scientific Research. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation. We thank the JCMT staff for their support of the GBS team in data collection and reduction efforts.

This research has made use of NASA's Astrophysics Data System (ADS) and the facilities of the Canadian Astronomy Data Centre (CADC) operated by the NRC with the support of the Canadian Space Agency (CSA). We have made use of the SIMBAD database (Wenger et al. 2000), operated at CDS, Strasbourg, France. The Starlink software (Currie et al. 2014) is supported by the East Asian Observatory. This research made use of APLpy, an open-source plotting package for Python hosted at http://aplpy.github.com.

Facilities: Herschel (SPIRE and PACS) - , JCMT (SCUBA-2). -

Software: APLpy (Robitaille & Bressert 2012), Astropy (Astropy Collaboration et al. 2013), NumPy (van der Walt et al. 2011), SciPy (Jones et al. 2001), Starlink (Currie et al. 2014), Matplotlib (Hunter 2007), Python.

APPENDIX: ANTI-CORRELATED β – T_d UNCERTAINTIES

Anti-correlation between β and T_d can be an artifact of SED fitting using the minimization of χ² method due to β −T_d degeneracy. Indeed, Shetty et al. (2009a, 2009b) demonstrated that SED fitting using the minimization of χ² method can lead to such an anti-correlation when modest amounts of noise are present in data. Such anti-correlations due to erroneous fittings have also been noted in previous works (e.g., Keene et al. 1980; Blain et al. 2003; Sajina et al. 2006; Kelly et al. 2012; Sadavoy et al. 2013).

The anti-correlated β − T_d uncertainties arise naturally from both β's and T_d's ability to shift the peak of a modified blackbody SED (see Equation (1)). While β is responsible for producing the power-law slope of the Rayleigh–Jeans tail of an SED, it can also shift the SED peak independent of T_d. An underestimation of β, for example, leads to a shallower power-law tail that pulls the SED peak lower in frequency. T_d, on the other hand, shifts the SED peak through Wien's displacement law, which acts on the blackbody component of the SED. The position of the SED peak can thus be maintained at a constant value by anti-correlating β and T_d, resulting in a degeneracy in reduced χ² SED fitting when the height of the SED can be arbitrary scaled by a free τ parameter.

Figure 13 shows the 50% and 75% probability contours of T_d and β drawn from SED fits to noisy modified blackbody models using the Monte Carlo approach with flux uncertainties of σ = 5% and 10%. The true T_d and β values behind the models are marked with "X" symbols. The modeled SEDs are sampled at the same wavelengths as our Herschel and JCMT bands for the fits shown in panels (a) and (c), and at the wavelengths used by Shetty et al. (2009a, i.e., 100, 200, 260, 360, and 580 μm) based on observations made by Dupac et al. (2003) for the fits shown in panels (b) and (d). Consistent with the results of Shetty et al. (2009a), anti-correlated uncertainties between T_d and β are significant for SED fits to these wavelengths, even when the flux noise is relatively modest (σ = 5%). These uncertainties tend to increase with temperature as the SED peak shifts upwards away from the sampled wavelengths. The fits to SEDs are poorly constrained at T_d ≳ 30 K for a noise level of σ = 5% and at T_d ≳ 20 K for σ = 10%.

By having a broader wavelength coverage and a better sampling of the longer wavelengths than those used by Dupac et al. (2003), we are able to achieve smaller T_d and β uncertainties at lower temperatures (T_d ∼ 10 K; see Figure 13) at a given noise level. With 160 μm as our shortest wavelength band instead of the 100 μm used by Dupac et al., however, our uncertainties due to flux noise at higher temperatures are larger than those of Dupac et al. due to having fewer samples on the SED peak. The preferential sampling of the Rayleigh–Jeans tail of the SED instead of the peak, nevertheless, is advantageous in making SED fits more robust against temperature variation along the line of sight (Shetty et al. 2009b).

The median flux noises in our data over where SED fits are made are small (<1% in Herschel bands and ∼5% in SCUBA-2 850 μm band) relative to our assumed flux calibration uncertainties (20% for PACS 160 μm band and 10% for the other bands). Indeed, we found the inclusion of measured rms noise in addition to the flux calibration uncertainties (correlated within each instrument) to have a negligible effect on the T_d and β uncertainties derived using the Monte Carlo method with the NGC 1333 data. We therefore assumed flux noise to be negligible with respect to the flux calibration uncertainties, which are systematic across a map, in our final uncertainties analysis. Due to the color calibration uncertainties being of the same order as the flux noise, we also assumed the color calibration error to be negligible.

Figure 14 shows the 50% and 75% probability contours of T_d and β drawn from SED fits to a few selected pixels in our NGC 1333 data using the Monte Carlo approach. Instead of adding random noise to each individual band, random flux calibration corrections are generated for each instrument and the bands within each instrument received the same flux correction at each iteration of the Monte Carlo simulation. We adopted a flux calibration error of 20% for PACS, 10% for SPIRE, and 10% for SCUBA-2. The best-fit T_d and β values are marked with "X" symbols. The T_d and β uncertainties due to flux calibration we quoted in this paper are effectively derived from fitting a Gaussian to these probability contours collapsed onto one of these axes (see Sadavoy et al. 2013 for details).

As shown in Figure 14, even though the probability contours of these selected pixels collectively resemble the shape of the T_d − β scatter plot of NGC 1333 shown in Figure 8(d), the probability contours of each pixel are only a fraction in size of the overall distribution. Furthermore, due to the flux calibration uncertainties being systematic across a map, erroneous offsets in our derived parameters should be spatially correlated on a scale of ∼30', which is the size of our smallest maps (i.e., SCUBA-2 maps). The relative pixel to pixel uncertainties between the derived parameters in each map should therefore be smaller than the stated absolute uncertainties demonstrated in Figure 14.

The derived T_d values are well constrained at low temperatures (≲20 K), which include most of our pixels (see Figure 4). At these low temperatures, the uncertainties in the derived β are relatively large with respect to the range of derived β values, and decrease only slightly at higher temperatures. Despite these large uncertainties, low-, intermediate-, and high-β values can still be distinguished from each other at a 1σ confidence interval.

Our adopted flux calibration uncertainties for extended sources are relatively conservative compared to the values estimated for point sources (<7% for PACS bands, Balog et al. 2014; ∼5% for SPIRE bands, Bendo et al. 2013; <8% for SCUBA-2 850 μm band). The most common T_d values found in all Perseus clumps, except for IC 348, agree with each other and with the ammonia gas temperatures of Perseus dense cores (Rosolowsky et al. 2008) to within 0.5 K. In comparison, our estimated absolute T_d uncertainties at these temperatures is ∼0.8 K. This result suggests that the flux calibration uncertainties we assumed are indeed relatively conservative.

Figure 15 shows the 50% and 75% probability contours of τ₃₀₀ and β drawn from SED fits to the same pixels as in Figure 14. As expected, τ₃₀₀ and β are positively correlated, given that the uncertainties between T_d and β and between T_d and τ are anti-correlated. The anti-correlation found in the latter is due to the flux of a modified blackbody being positively correlated with T_d (see in Equation (1)). We find the slopes of the correlated uncertainties to be relatively well defined in the log(τ) − β space for well constrained SED fits (i.e., at T_d ≲ 30 K), and relatively constant regardless of the fitted parameters. Given these correlated uncertainties, we did not find a significant underlying correlation between β and τ₃₀₀ shown in Figure 9.

THE JCMT GOULD BELT SURVEY: EVIDENCE FOR DUST GRAIN EVOLUTION IN PERSEUS STAR-FORMING CLUMPS

Article metrics

Permissions

Author affiliations

ORCID iDs

Dates

ABSTRACT

1. INTRODUCTION