© ESO, 2012

1. Introduction

During the formation of a low-mass star, the first hydrostatic core (FHSC) phase is one of the shortest phases, lasting only 10²−10³ yr, depending on how the collapse proceeds (Bate 2011). The FHSC phase is characterized by the presence of a central object with the mass of a giant planet that has a size of just a few astronomical units. This phase lasts until the central object’s temperature reaches  ~2000 K. At that point, the molecular hydrogen dissociates, causing a second collapse that brings the central object to the typical size of a protostar, starting the Class 0 phase (André et al. 1993). The short duration of the FHSC phase makes it difficult to observe an FHSC in a star-forming region.

Before the FHSC phase, the star-forming core (a “pre-stellar core”) has a spectral energy distribution (SED) that can be modelled as a greybody spectrum. After this FHSC phase, the SED in the far-infrared (FIR) still resembles a greybody spectrum, but the emission of the central forming protostar starts to be visible at shorter wavelengths, i.e., λ < 70 μm. Theoretical models (e.g., Omukai 2007; Saigo & Tomisaka 2011) show that the SED deviates from a pure greybody shape during the FHSC phase at wavelengths λ ≲ 200 μm. At even shorter wavelengths, e.g., for λ ≲ 30 μm, the emission is still too faint to be detected. Taking into account that the photometric bands ot the Spitzer MIPS instrument (Rieke et al. 2004) were centred in the FIR at 24 μm, 70 μm, and 160 μm, the expected signature for a 1 M_⊙ FHSC SED includes a detection at 70 μm and the non-detection at 24 μm (Commerçon et al. 2012).

None of the FHSC candidates meets these photometric requirements. In the Perseus star-forming region, the candidate [EYG2006] Bolo 58 (hereafter Per 58) (Enoch et al. 2010; Dunham et al. 2011) is visible in the 24 μm Spitzer map. Similarly, the candidate [RNC96] Cha-MMS 1 (herafter MMS1) in Chamaeleon (Belloche et al. 2006; Cordiner et al. 2012) is observed both at 24 μm and at 70 μm. The proposed FHSC candidates LDN 1451-mm (Pineda et al. 2011), [EYG2006] Bolo 45 (Schnee et al. 2012), and [CB88] 17 (Chen et al. 2012) are indeed not visible at 24 μm, but they remain undetected also at 70 μm in the Spitzer maps, making the SED-based classification less firm. Finally, L1448-IRS2E (Chen et al. 2010) was not even detected at 160 μm. It must be said, however, that the knowledge of the SED alone is not enough to firmly identify an FHSC. Spectral lines and interferometric observations are necessary and, indeed, these FHSC candidates were not proposed based on continuum data only.

Furthermore, the (non)detection in a certain band clearly depends on the distance to the source. To predict the expected flux from an FHSC, the distance of the star-forming region in Taurus (150 pc) is often taken as reference distance (e.g., Tomida et al. 2010; Commerçon et al. 2012). Moving to larger distance, however, strengthens the requirement on non detection at 24 μm, while the detection at 70 μm may become less stringent. Indeed, the non-detection at 70 μm has been explained by Enoch et al. (2010) on the basis that an FHSC could, in principle, be detected with Spitzer, but the “cores to disks” (c2d) survey (Evans et al. 2003) was not sufficiently sensitive to detect these sources. The recent launch of the Herschel Space Observatory satellite (Pilbratt et al. 2010) has now opened up the possibility to observe in the FIR bands between 70 μm and 500 μm, with unprecedented sensitivity and spatial resolution. Among the observing programs, the Herschel Gould Belt survey (GBS; André et al. 2010) aims to obtain a complete census of pre-stellar cores and Class 0 sources in the closest star-forming regions. Since FHSCs can be considered to be extremely young Class 0 sources, we expect to find new FHSC candidates among the numerous objects discovered with Herschel.

One of the target clouds of the GBS is the star-forming region in Perseus molecular cloud, at an average distance of  ~235 pc (Hirota et al. 2008). Perseus hosts low- and intermediate-mass young stellar objects (YSOs), which places it roughly between a low-mass star-forming cloud like Taurus, and a high-mass star-forming cloud like Orion.

As a starting point for our analysis of these observations, we looked for sources detected in the Herschel map at 70 μm without any counterpart in the corresponding Spitzer 24 μm map. We found one such source and in this paper discuss the physical properties that can be derived from the analysis of its SED. This object is spatially associated with the source [HKM99] B1-bS (hereafter B1-bS) discovered by Hirano et al. (1999). They found this source through interferometric observations carried out with the Nobeyana Millimeter Array (NMA) at 3 mm. They also found a second source, [HKM99] B1-bN (hereafter B1-bN),  ~20′′ north of B1-bS, clearly visible in the Herschel PACS maps at λ ≥ 100 μm. Combining the NMA data with other single-dish observations between 350 μm and 2 mm, and spectral observations of H¹³CO⁺, Hirano et al. (1999) concluded that these two sources were younger than already known Class 0 objects. Their observations, however, were all taken in the Rayleigh-Jeans part of their SEDs, causing a large uncertainty in the determination of the temperature and, in turn, of the mass. Herschel data now extend the knowledge of the SEDs exactly where the peak of the emission falls.

On the basis of Herschel fluxes, which are incompatible with a greybody, we reconsider the physical properties of these two sources, and we propose they are new promising FHSC candidates. The rest of the paper is organized as follows: in Sect. 2, the observations are presented and the algorithm of source detection, CuTEx, is briefly described. In Sect. 3, we fit the SEDs with a two-component model that describes the FHSC emission in a simple way. We show that alternative approaches fail to reproduce the SEDs. In Sect. 4, we discuss the results, and, in Sect. 5, our work is summarized. Appendix A gives a detailed description of the SED fitting, including ancillary explanations on the colour corrections and on the Gaussian fitting of sources detected with Herschel. We also present a new formula to derive the mass of an object whose SED is fitted with a greybody. Appendix B describes an alternative approach to the sigma-clipping algorithm for outlier detection, where the threshold for the detection depends on the size of the sample.

2. Observations and data analysis

The observations of Perseus cover a total area of about 10 deg², mapped with the PACS (Photodetector Array Camera & Spectrometer; Poglitsch et al. 2010) and SPIRE (Spectral and Photometric Imaging Receiver; Griffin et al. 2010; Swinyard et al. 2010) instruments using the parallel mode with the telescope scanned at 60′′/s. In this observing mode, the two instruments observe almost the same field with SPIRE, delivering data in its three bands (250 μm, 350 μm, and 500 μm), and PACS in two out of three bands, 70 μm and 160 μm in our case. A smaller area, 6.7 deg², was mapped again with PACS at 100 μm and 160 μm only, with the telescope scanned at 20′′/s, to reach a higher spatial resolution and a better sensitivity.

The observations we consider in this paper cover a region of about 5.6 deg², centred on NGC 1333. Table 1 summarizes the observations. Preliminary results from the parallel mode observations were recently presented by Sadavoy et al. (2012), and by Bressert et al. (2012). A complete description of this region based on Herschel data will be given in a future paper (Pezzuto et al., in prep.).

Data were processed with version 8.0.1559 of HIPE (Ott 2010), except for the deglitching, which was performed with a sigma-clipping algorithm with a variable threshold detailed in Appendix B, and the map reconstruction, which was made using the ROMAGAL code (Traficante et al. 2011). The zero point of the flux calibration for extended emission was established by comparing Herschel data with Planck and IRAS data over the same area (Bernard et al. 2010). The result of the comparison is a set of offsets, one for each band, algebraically added to the observed images. For photometric measurements, however, where a background is estimated and subtracted from the source flux, as in this work, the zero points are not important.

2.1. Source detection and photometry

To identify sources in the field, we adopted the CuTEx algorithm (Molinari et al. 2011). This method computes the second derivative of the input image along four different directions. In the “curvature” image so obtained, we looked for pixels where the second derivative is negative and its absolute value is higher than a user defined threshold with respect to the local fluctuactions. Isolated high-curvature pixels are rejected, and only groups of contiguous pixels, whose number depends on the sampling of the beam, were considered in the following analysis. The curvature in each group was analyzed to verify the statistically significance of the local peak. The detection was made independently at each wavelength. Once a preliminary list of candidate sources was derived, photometry was recovered by fitting a 2D Gaussian profile. The position of the peaks of the candidate sources and their sizes, as derived from the curvature map, were used as initial guesses for the parameters of the 2D Gaussian. If two or more sources were closer than twice the width of the point spread function (PSF) at a given wavelength, their Gaussians were fitted simultaneously.

For λ ≥ 160 μm, the diffuse emission contributes a strong, variable signal across all spatial scales. An important question is then how to estimate the contribution of the diffuse background at the source position. CuTEx models this emission by fitting a linear 2D-polynomial (i.e., a plane) inside a subregion centred on each source, with a size of six times the PSF at the specific wavelength. The plane and the Gaussian(s) are fitted simultaneously. There is, however, no general consensus on how the background can be estimated: different approaches can give different results and one of the main, and to large extent unknown, uncertainties in the derived fluxes, and in turn on the shape of the SED, comes from the background subtraction. For this reason is important to compare the results of flux extraction with more than one method. To this aim, we also used getsources (Men’shchikov et al. 2012) to derive the fluxes of our sources.

Among the sources detected at 70 μm, we found that only one was not detected in the corresponding Spitzer 24 μm map. This source is one of a pair of very young objects first discovered at millimetre wavelengths by Hirano et al. (1999) and dubbed B1-bS and B1-bN. Our 70 μm source corresponds positionally to B1-bS, while B1-bN is clearly detected between 100 μm and 250 μm. The separation between the two sources is about 20′′, larger than the beam of PACS in all bands. It is also larger than the beam of SPIRE at 250 μm (~18′′), but smaller than the beams at 350 μm and 500 μm (~25′′and  ~36′′, respectively). We attempted to fit two Gaussians to the SPIRE data using the positions of B1-bS and B1-bN at 160 μm as initial guesses. During the fitting procedure the initial distance between the two objects was kept constant.

Herschel data were complemented with archival Spitzer data from the c2d survey (Evans et al. 2003), SCUBA Legacy Catalogue data (Di Francesco et al. 2008) at 450 μm and 850 μm, and the Bolocam data (Enoch et al. 2006) at 1.1 mm. We ran CuTEx on all these maps except for the Spitzer maps (see below). Again, the positions of the sources at 160 μm were used as initial guesses. The effective FWHM resolutions of the SCUBA Legacy Catalogue data at 450 μm is 17$\stackrel{\mathrm{″}}{.}$3, so that the two sources can be considered to be fairly separated. At 850 μm, however, the effective FWHM is 22$\stackrel{\mathrm{″}}{.}$9 and encompasses both B1-bS and B1-bN, causing the photometry to be more uncertain. At 1.1 mm, the Bolocam beam of 31′′ makes it impossible to distinguish the two objects.

Spitzer source [EDJ2009] 295 (hereafter S295) is very close to B1-bS. It is visible in Herschel maps at λ ≤ 100 μm. S295 is commonly associated with the millimetre-wavelength source [EYG2006] Bolo 81 (Enoch et al. 2006). Since B1-bS is close to the wings of the Spitzer PSF of S295, it could be that B1-bS is not seen at 24 μm because of its proximity to the bright S295. To test this possibility, we used DAOPHOT on the Spitzer map to perform PSF photometry of S295, a method only possible due to the low background level at 24 μm. B1-bS remains undetected also after S295 is subtracted; from the analysis of the detections, we concluded that any source with a flux  ≥0.2 mJy, corresponding to a  ≥5 σ detection, should have been detected, and so we assumed 0.2 mJy flux as an upper limit at 24 μm for both B1-bS/bN. At shorter wavelengths, the distance between the B1-b sources and S295 is large enough and the background is so low and smooth, that to estimate an upper limit for the fluxes in the IRAC bands we just computed the rms of the background over a region of 11 × 11 pixels, centred at the 160 μm positions. We used this value as the 1σ upper limit for the peak flux which, along with the instrumental FWHM, can then be used to derive an upper limit at the IRAC wavelengths for a 2D Gaussian profile. The upper limits for B1-bS are systematically higher than those of B1-bN. Since the former is close to the brighter S295, it is possible that the B1-bS upper limits suffer some flux leakage from S295. For this work, however, shifting the upper limits by some percentage does not impact the fitting result (see below). The full set of measured fluxes is reported in Table 2. Image cut-outs at all wavelengths centred on B1-bS, 1$\stackrel{\mathrm{\prime }}{.}$5  ×  1$\stackrel{\mathrm{\prime }}{.}$5 in size, are shown in Fig. 1.

The positions of the B1-b sources from 24 μm to 1.1 mm are shown in Fig. 2. The triangle denotes the position of S295 taken from the Spitzer c2d catalogue (Evans et al. 2003), which agrees with our S295 positions at 70 μm and 100 μm. The positions of the B1-b sources given in Hirano et al. (1999) are precessed to J2000. The position of Bolo 81 is taken from Enoch et al. (2006). This source was later detected by Rosolowsky et al. (2008) in a survey of ammonia cores, and dubbed [RPF2008] NH3SRC 123. The cross labelled CuTEx Bolocam is the position we find using CuTEx on the 1.1 mm map. We also report in the figure the positions of the “core” and “outflow” recently found by Öberg et al. (2010) from CH₃OH observations. The position they name “core” is halfway between Bolo 81 and B1-bN, while their “outflow” is slightly SE of B1-bS.

To derive the Herschel coordinates of the B1-b sources we proceeded as follows. First, we averaged the PACS positions of S295 and compared them with the coordinates reported in the c2d catalogue, finding an offset of $\mathrm{\Delta }\mathit{\alpha }\mathrm{=}\mathrm{1}\stackrel{\mathrm{″}}{.}\mathrm{4}$ and $\mathrm{\Delta }\mathit{\delta }\mathrm{=}\mathrm{0}\stackrel{\mathrm{″}}{.}\mathrm{7}$. Assuming that the c2d coordinates are more accurate, these offsets can be considered estimates of the absolute value of the PACS pointing errors in our maps. Next, we averaged the PACS coordinates of B1-bS after adding the offset found for S295, finding α = 3^h33^m21$\stackrel{\mathrm{″}}{.}$3, δ =  + 31^d7′27$\stackrel{\mathrm{″}}{.}$4, with a total uncertainty (1σ) of 1$\stackrel{\mathrm{″}}{.}$1. We did not average the SPIRE positions because at 350 μm and 500 μm the two sources are more blended and the Gaussian fits are less reliable. We used the SPIRE 250 μm position to estimate the offset between PACS and SPIRE: $\mathrm{\Delta }\mathit{\alpha }\mathrm{=}\mathrm{4}\stackrel{\mathrm{″}}{.}\mathrm{2}$ and $\mathrm{\Delta }\mathit{\delta }\mathrm{=}\mathrm{3}\stackrel{\mathrm{″}}{.}\mathrm{0}$. Finally, for B1-bN, we applied to the PACS 100 μm and 160 μm coordinates the same offsets found for S295; and to the SPIRE position at 250 μm the offsets found between PACS and SPIRE for B1-bS. The three pairs of coordinates were averaged finding α = 3^h33^m21$\stackrel{\mathrm{″}}{.}$2, δ =  + 31^d7′44$\stackrel{\mathrm{″}}{.}$2, with a total uncertainty (1σ) of 3$\stackrel{\mathrm{″}}{.}$7. Note that all the offsets are much smaller than the separations between the objects, so that source confusion due to a mis-pointing is excluded.

Although it was already noted in the past that S295 and B1-bS are spatially separated, in the literature there is some tendency to confuse the two sources. From Figs. 1 and 2, it is clear that these two objects are different. The association of Bolo 81 with S295 is very dubious since the latter is not detected longward of 100 μm. The large offset,  ~20′′, between Bolo 81 and the B1-b sources also casts some doubts on the possible association between these objects.

3. Results of SED fitting

The simplest model that can be fitted to the data is an optically thin, isothermal greybody, which is appropriate for starless and pre-stellar cores ${\mathit{F}}_{\mathit{\nu }}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathit{\nu }}\mathit{B}\mathrm{\left(}\mathit{T,\nu }\mathrm{\right)}}{{\mathit{D}}^{\mathrm{2}}}\mathit{M}$where κ_ν is the opacity (see the discussion in Appendix A.3), B(T,ν) is a blackbody at temperature T, D is the distance, M is the mass. Following the approach used for other GBS obervations (e.g., Könyves et al. 2010), the dust opacity index was fixed to 2. The best-fit models, considering only data at λ ≥ 160 μm, have T = 9 K, M = 7.3 M_⊙ for B1-bS, and T = 8 K, M = 9.4 M_⊙ for B1-bN.

We had to exclude the fluxes at short wavelengths from the greybody fitting because in a two-colour diagram of [100−160] vs. [250−350] the two B1-b sources have colours that are bluer than those of a blackbody. Such a colour could only be reproduced with a greybody given an unrealistic situation where the dust opacity increases with (far-IR/submm) wavelength. If the entire SEDs are fitted with a greybody, as in Eq. (A.1) and leaving all the parameters free to vary, the best-fit models have β = 0, and T = 18 K for B1-bS; and β = 0 and T = 15 K, for B1-bN. During the fitting procedure, we imposed the condition β ≥ 0, and the value β = 0 best matches the condition β < 0, corresponding to a dust opacity increasing with wavelength.

Another way of fitting an SED with colours that are bluer than those of a blackbody is by assuming that the observed emission is caused by two components (see, e.g., Fig. 2 in Pezzuto et al. 2002). For this reason, we fitted the B1-b derived SEDs by adding the contributions of a blackbody embedded in a dusty envelope whose emission is modelled with a greybody. We give a physical interpretation of this two-component model below.

The details of the fitting procedure, as well as some additional information on colour corrections applied to the PACS and SPIRE data, are given in Appendix A. The best result of the two-component fitting procedure is shown in Fig. 3 for both B1-b sources, and the corresponding parameters of the best-fit models, T_b and Ω_b for the blackbody component, and T_g, λ₀, β, and Ω_g for the greybody component, are listed in Table 3.

Fig. 1 Maps of 1$\stackrel{\mathrm{\prime }}{.}$5  ×  1$\stackrel{\mathrm{\prime }}{.}$5 of B1-bS and B1-bN at all wavelengths, with sizes and positions derived from the Gaussian fits. The FWHM of the respective instrumental beam is shown as a circle in the bottom-left corner of each respective panel. At λ ≤ 100 μm, the source S295 is also visible. Colour bars are in Jy/beam for the SCUBA maps, and MJy/sr in all the others. In the Spitzer images, positions and sizes of the sources are copied from the 70 μm map. IRAC maps are in the top row, from left to right: a) 3.8 μm; b) 4.5 μm; c) 5.8 μm; d) 8.0 μm. Second row: e) Spitzer 24 μm; f) PACS 70 μm, B1-bN is not really detected and its flux corresponds to a 1σ rms point source; g) PACS 100 μm; h) PACS 160 μm, S295 is no longer visible. Third row: i) SPIRE 250 μm; j) SPIRE 350 μm; k) SCUBA 450 μm, where the beam shown is 11′′, see Appendix A.2; l) SPIRE 500 μm. Bottom row: m) SCUBA 850 μm; n) Bolocam 1.1 mm, where only one source has been fitted.

Fig. 2 Positions (J2000) of the sources at all wavelengths; see text for references. Two black open circles are centred on the position of B1-bS and B1-bN after applying a positional offset correction found by comparing the PACS and Spitzer positions of S295, see text for details. The radius of the circles is the 1σ of the mean of the PACS coordinates, for B1-bS, and of the mean of PACS 100 μm, 160 μm, and SPIRE 250 μm coordinates for B1-bN. The positions at 350 μm, 500 μm, and 850 μm are very uncertain, because B1-bS and B1-bN are not resolved at those wavelengths. The blue open circle centred on Bolo 81 has a radius of 7′′, the pointing accuracy of Bolocam observations (Enoch et al. 2006).

The best-fit blackbody has a temperature of T ~ 30 K in both sources. Its size is about 130 AU at 235 pc for B1-bS, and is much smaller in B1-bN. For the latter, Ω_b is poorly constrained since we do not have a 70 μm detection. The parameters of the greybody component, i.e., the dusty envelope, are similar for both sources. The temperature is 8 K, not too much different from the 11.7 K temperature found by Rosolowsky et al. (2008) for the ammonia core NH3SRC 123 detected at the position of Bolo 81.

From the relation $\mathit{ṁ}\mathrm{=}\mathit{\alpha }{\mathit{c}}_{\mathrm{s}}^{\mathrm{3}}\mathit{/}\mathit{G}$, where c_s is the sound speed, we can estimate the mass accretion rate. There are two limiting cases for the gravitational collapse of an isothermal sphere (see McKee & Ostriker 2007, for a detailed discussion and references): in the fast collapse (the so-called Larson-Penston-Hunter solution) the infall is supersonic and α = 47, while in the slow collapse (or Shu’s solution) the infall is sonic and α = 0.975. For T = 8 K we derive in the former case ṁ = 5.3 × 10^-5 M_⊙/yr, and ṁ = 1.1 × 10^-6 M_⊙/yr for the latter case. Taking into account that c_s could also include contributions from a magnetic field or from turbulence, we can safely conclude that ṁ < 1 × 10^-4 M_⊙/yr.

Turning these mass accretion rates into an accretion luminosity is not straightforward because we need to know the mass and the radius of the accreting source. Because we modelled the emission of the central object with a blackbody, we cannot derive its mass. The value we found for the radius from the fitting procedure appears to be higher than the few AU expected for an FHSC: 130 AU for B1-bS and 48 AU, but with a 50% error, for B1-bN.

Fig. 3 Best-fit two-component models for B1-bS and B1-bN. Filled circles are Herschel data (for B1-bN the 70 μm flux is an upper limit), open circles are SCUBA data, and upper limits at λ ≤ 24 μm are from Spitzer IRAC and MIPS. The solid line is the best-fit model, a sum of a blackbody (dashed line dominating at shorter wavelengths) plus a greybody (dashed line more significant at longer wavelengths). Herschel fluxes are those reported in Table 2 multiplied by the colour correction factors, see Appendix A.1.

Such a large radius, however, can have a physical meaning. Bate (2011) found that in a rotating envelope a disc with a size of  ~100 AU can form before the stellar core; in the Saigo & Tomisaka (2011) model, again for a rotating core, 100 AU correspond to the radius where the envelope’s density profile changes from a r^-2 power law to a flatter profile that ends on the core’s surface, located at 20 AU when the FHSC first forms. The inner component of our model, then, can be physically interpreted as describing the thicker and inner regions of the envelope that eventually merge with the compact central objects. In other words, the parameters T_b and R_b of the blackbody give the average properties of this hybrid component, without being able to distinguish where the envelope ends and where the core starts.

The best-fit value for β, the exponent of the dust opacity, is 2, even if good fits for B1-bN can be obtained also with β = 1.5; see Appendix A. The wavelength at which τ = 1, λ₀, is 100 μm for B1-bS and a value slightly higher for B1-bN.

The mass of the dusty envelope in B1-bS is 9    ±    5 M_⊙, with the large uncertainty mainly due to the uncertainty in β. The mass of B1-bN corresponding to the best-fit model is also 9 M_⊙, but the large uncertainty in Ω_b and Ω_g causes a wide range of variation for the mass: 1 ≤ M/M_⊙ ≤ 23. Previous estimates of the masses were made by Hirano et al. (1999) who, with data at λ ≥ 350 μm, derived smaller masses for both sources, M ~ 1.7 M_⊙, as well as higher temperatures, T ~ 18 K. Hatchell et al. (2007a) derived an upper limit for B1-bS M_235pc < 23.1 M_⊙, but they also used the millimetre flux of Bolo 81 for the SED.

The radius $\sqrt{{\mathrm{\Omega }}_{\mathrm{g}}\mathit{/}\mathit{\pi }}$ of the envelope is  ~20′′ for B1-bS and  ~17′′ for B1-bN. Such a radius, however, can not be immediately compared with the sizes given as FWHM in Table 2, which result from fitting the sources in the map with a 2D-Gaussian profile. Indeed, when fitting the SED, we assumed that the source has a finite radius, R, while a Gaussian profile does not have any radius. To compare R with the observed FWHMs, we proceeded as follows: the integral I of a normalized 2D circular Gaussian over a circle of radius R_I is (Wörz 2006) $\mathit{I}\mathrm{=}\mathrm{1}\mathrm{-}{\exp }^{\mathrm{-}\frac{{\mathit{R}}_{\mathit{I}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\mathit{.}$Clearly, I = 1 only for R_I → ∞, but already when R_I = 3σI ~ 0.99¹. We then defined the radius R_G of a 2D Gaussian as R_G = 3σ, or, using the relation between FWHM and σ, R_G = 1.29 × FWHM. For B1-bS, the best-fit radius R is 20$\stackrel{\mathrm{″}}{.}$4, a circular 2D-Gaussian with R_G = R has FWHM = 15$\stackrel{\mathrm{″}}{.}$8, a value in between the measured size at 250 μm and at 350 μm. For B1-bN, the best-fit radius is 17′′, which corresponds to FWHM = 13′′, very close to the measured size at 250 μm.

It is evident from Table 2, however, that the size derived from the observations is poorly defined, since the FWHM increases with wavelength. This trend is to some extent related to the instrumental beam size. For instance, the FWHM of B1-bS is smaller at 100 μm than at 70 μm, which is likely related to the size of the PSF of the intrument, that at 100 μm (scanning speed of 20′′/s, nominal PACS compression mode) is smaller than at 70 μm (scanning speed of 60′′/s, double PACS compression mode). This trend is also likely to be related to the increase with wavelength of the background level, which makes it more difficult to distinguish the genuine source emission from the extended emission. It is physically plausible, however, to assume that this trend is, at least in part, intrinsic to the sources, because this phenomenon has been already observed elsewhere.

In high-mass star-forming regions, for instance, an increase of the size of the dense cores with wavelength is known and is taken into account by scaling the flux proportionally to the ratio of measured radii, assuming the size at 160 μm as the fiducial radius (see, e.g., Motte et al. 2010; Giannini et al. 2012). This approach is justified by theoretical models showing that the emission can be described as coming from an almost isothermal envelope whose mass increases linearly with the radius.

For cores in nearby star-forming regions, where the spatial resolution is sufficiently high to resolve the structure of the envelope, the wavelength-radius relationship could be a consequence of temperature stratification in the envelope. Namely, at shorter wavelengths, we see the inner and warmer part of the envelope. At longer wavelengths, however, the emission comes from the outer and generally colder part of the core. In our model, part of the blackbody radiation is absorbed by the envelope, which should cause a temperature stratification. We did not treat the radiative transfer of the system, however, assuming instead an isothermal envelope for simplicity. The temperature T_g so derived, and in general the whole set of parameters T_g,β,λ₀,Ω_g, along with the mass M from Eq. (A.9), can then be considered to be the averaged properties of the outer regions of the envelope, in the same way as T_b and R_b describe average properties of the inner regions, as discussed above. In other words, while our objects likely have a temperature stratification that causes an inrease of the observed radius with wavelength, our model shows that the SEDs can be described in a simple way as having one radius and one temperature at all wavelengths.

The consistency of our approach can be seen from the fact that 1) the derived size is comparable to the observed sizes at λ ~ 300 μm, and 2) λ₀ ≳ 100 μm. Indeed, in a greybody model, the dependence on the solid angle disappears in the optically thin regime, so that the derived Ω_g values are most likely to represent the sizes corresponding to the longest wavelength where the envelope is not yet completely thin. The robustness of our result can be verified by noting that T_g and M are very similar to the values found with the isothermal greybody model discussed at the beginning of this section.

After modelling the observed SED with two components we can separately derive the internal luminosity L_int of the sources and the bolometric luminosity L_bol.

In Table 4 we report L_b, the luminosity of the blackbody component. This is found with the standard equation ${\mathit{L}}_{\mathrm{b}}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{R}}_{\mathrm{b}}^{\mathrm{2}}\mathit{\sigma }{\mathit{T}}_{\mathrm{b}}^{\mathrm{4}}$ (the parameters are those of Table 3 with ${\mathit{R}}_{\mathrm{b}}\mathrm{=}\mathrm{\left(}\sqrt{{\mathrm{\Omega }}_{\mathrm{b}}\mathit{/}\mathit{\pi }}\mathrm{\right)}$). At any frequency ν, the envelope absorbs a part of the blackbody emission, leaving the amount B_ν(T_b)e^−τ_ν free to escape from the envelope. The integral over frequency of this radiation is reported in the table as L_em. Clearly, the fraction (1−e^−τ_ν) is absorbed inside the envelope. The column labelled L_g reports the luminosity of the greybody component found by numerical integration of Eq. (A.1) in the range 1   μm ≤ λ ≤ 10 mm. The observable luminosity predicted by the model is then L_em + L_g, which we report in the column L_bol. It can be compared with the observed luminosity that we report in the column L_SED, found by integrating the measured fluxes.

If there is no external source of energy, then by the conservation of energy L_g = L_b − L_em. The contribution of the interstellar radiation field to the sources luminosity could be, then, estimated as L_ISRF = L_g − (L_b − L_em). For the two sources, however, we find quite different results: L_ISRF = 0.03   L_⊙ for B1-bS and 0.15 L_⊙ for B1-bN. While a small difference between L_ISRF in the two sources could be reasonable, this large discrepancy is more likely caused by the uncertanties in the best-fit parameters, which are in turn due to the uncertanties in the measured fluxes, or to the consequence of the non-proper treatment of the radiation transfer in our model.

Finally, as an example of a model that treats the radiative transfer more rigorously, we fitted the SEDs using the on-line fitting tool by Robitaille et al. (2007), which compares the observed fluxes with a grid of theoretical SEDs. This comparison is important also because in this way we can test the hypothesis that our sources are more evolved than an FHSC. The models by Robitaille et al. (2007) are indeed appropriate to describe the emission from an evolved YSO, where a central star already formed. On the contrary, their grid does not contain models that are appropriate for very young sources like an FHSC. We found that no model can account for both the upper limits in the Spitzer bands and the far-IR/submm fluxes. The SEDs of our sources are not compatible with any of the 200 000 models of YSOs in the grid of Robitaille et al. (2007).

4. Discussion

The following features can be extracted from the observed SEDs of our sources, prior to any modelling: the colours in the PACS bands are not compatible with those of a greybody, while the SEDs of the starless or pre-stellar cores can be modelled with a greybody. An additional blackbody is required, which hints at the presence of an inner and warmer compact component. This second component is not visible at λ ≲ 70 μm, therefore it is not as evolved as other nearby protostellar objects. Another observational signature, although indirect, of the very young age of B1-bS comes from the statistics of the YSOs detected in Perseus. The preliminary analysis with CuTEx of the western region of Perseus gives a list of  ~40 tentative sources detected in all six Herschel bands. They are visible in the 24 μm Spitzer map, the only exception being B1-bS. This fact alone points to some peculiarity with this source.

In Table 5, we report the photometric data available in literature, or derived by us, for all FHSC candidates in the fundamental range between 24 μm and 160 μm. As already stated in Sect. 1, the SED only is not sufficient to identify an FHSC. The other sources reported in Table 5 have been proposed as FHSC candidates based on other observational properties, e.g., the u-v visibilities from interferometric observations. On the other hand, it is evident that as far as the SED is concerned, B1-bS remains an exception also when compared with the other candidates. It is the only source that shows clear evidence of an SED more evolved than a pre-stellar core, but less evolved than a Class 0 source.

Other age indicators directly derived from the SEDs are the ratio L₃₅₀/L_bol, i.e., the luminosity at λ ≥ 350 μm over the bolometric luminosity, and T_bol, i.e., the temperature of a blackbody with the same mean frequency as the observed SED. For a Class 0 source, L₃₅₀/L_bol > 5% (André et al. 1993) and T_bol < 70 K (Chen et al. 1995). Our sources have quite extreme values with respect to these limits, which confirms that they are indeed very young. B1-bS has L₃₅₀/L_bol = 25% and T_bol = 18 K; B1-bN has L₃₅₀/L_bol = 35% and T_bol = 14 K. From these signatures alone, however, it is not possible to distinguish an FHSC from a pre-stellar core, or a starless core.

In summary, B1-bS shows all features expected in the SED of an FHSC. B1-bN, the companion 20′′ north of B1-bS, is not detected at 70 μm, but otherwise it shares all characteristics of B1-bS.

For both sources, only the bolometric luminosities, found by integrating the observed fluxes, do not fit with FHSC predictions: L_bol is 0.49 L_⊙ and 0.28 L_⊙ for B1-bS and B1-bN, respectively, in excess of the maximum predicted for an FHSC, L ≲ 0.1 L_⊙ (Omukai 2007). This high luminosity could be ascribed to imperfect background removal, but even if we limit the integration to the PACS bands, where the diffuse emission is less prominent, we find L_PACS = 0.14    L_⊙ for B1-bS. It then seems reasonable to expect L ≳ 0.2 L_⊙, a value that also excludes the possibility that B1-bS is a very low luminosity object (VeLLOs; e.g., Di Francesco et al. 2007).

The observed bolometric luminosity, however, can be higher than the internal luminosity because of external heating of the envelope by the interstellar radiation field. We have seen in the previous section that it is not easy to estimate L_ISRF from our data. Instead, we can make use of the relationship between L_int and F₇₀ found by Dunham et al. (2008). For B1-bS, their relationship gives L_int ~ 0.04 L_⊙, a value that seems difficult to reconcile with the observed luminosities. On the other hand, Dunham et al. (2008) derived the relation by modelling the protostar with a fixed temperature of 3 000 K and a luminosity in the range of 0.03–10 L_⊙, which implies a radius in the range 0.5–11 R_⊙. These values are reasonable for a Class 0 source, but could be inappropriate for a still-younger object, like an FHSC. Dunham et al. (2008)’s result, however, clearly implies that for our sources L_int is surely lower than the measured L_bol. The precise factor, however, remains unknown.

Recently, Commerçon et al. (2012) found that the luminosity of an FHSC becomes higher than 0.1 L_⊙ only at late times, so that the vast majority of FHSCs appear as VeLLOs for almost all their lifetime. Nonetheless, these authors also found cases where only a few hundred years after the formation of the FHSC, the predicted luminosity is as high as 0.3 L_⊙, similar to what we derive for B1-bN. This increase in luminosity is caused by the increase of mass of the inner core which in turn is caused by the accretion from the envelope.

Another possible contribution to the luminosities of our sources is through contamination of the SEDs by an outflow, as recently found by Maury et al. (2012) for the prototypical Class 0 protostar VLA 1623–243. In our cases, however, the existence of an outflow in B1-bS/bN is not completely clear, therefore it is not possible to explore the reliability of this hypothesis. Drabek et al. (2012) found that ¹²CO line emission can contribute to the dust continuum at 850 μm. They estimated, however, that in NGC 1333, at positions far from outflows, the contamination is less than 20%. For the B1 region, Sadavoy et al. (in prep.) found that CO (3−2) emission appears to be relatively minor towards B1-bS and B1-bN, therefore any contamination of the dust continuum from the gas is likely insignificant at 850 μm.

The results of the fits we performed must be interpreted with caution because of the uncertanties in the fluxes. We tentatively suggest, however, that B1-bN is still slightly younger than B1-bS. It has an envelope slightly more optically thick than that of B1-bS, explaining the lack of detection in the shortest PACS band. For both sources, the envelope is cold and massive, and its internal thermal pressure alone is insufficient to prevent its gravitational collapse. This conclusion comes from comparing the mass of the envelope and the corresponding critical Bonnor-Ebert mass m_BE, the largest mass that an isothermal sphere of gas bounded by pressure can have without collapsing. When M_env > m_BE, the internal thermal pressure may not be high enough to support the core against internal gravity and external pressure, so that the envelope is undergoing, or is about to collapse. We found for m_BE, see Eq. (A.3), 0.4 M_⊙ (B1-bS) and 0.3 M_⊙ (B1-bN). In both cases, M_env ≫ m_BE so that the envelopes cannot be supported by internal thermal pressure alone, even if we cannot exclude that other physical mechanisms such as the amount of turbulence or the strength of the magnetic field (see, e.g., Basu et al. 2009) may still contribute to support the cloud against the gravitational collapse.

As far as the magnetic field is concerned, however, the largest mass that a magnetic field of 31 μG (as estimated for the B1 region by Matthews & Wilson 2002) can support is only 0.11 M_⊙ (from Stutz et al. 2007), assuming a radius of 0.023 pc (20′′at 235 pc, see Ω_g in Table 3).

To estimate the stability against turbulence we computed the virial parameter α (Bertoldi & McKee 1992) and considered a core to be gravitationally bound if $\mathit{\alpha }\mathrm{=}\frac{\mathrm{5}{\mathit{\sigma }}^{\mathrm{2}}\mathit{R}}{\mathit{GM}}\lesssim \mathrm{2}\mathit{,}$where σ is the velocity dispersion found by adding in quadrature the non-thermal component of the NH₃ lines measured by Rosolowsky et al. (2008) for NH3SRC 123, and the thermal component of a mean particle of molecular weight μ = 2.33: σ_NT = 0.325 km s^-1 and σ_T = 0.204 km s^-1. From Table 3R = 7.19 × 10¹⁶ cm and R = 5.99 × 10¹⁶ cm for B1-bS and B1-bN, respectively, for a distance of 235 pc. Then, α = 0.44 and 0.37 for the two sources. Clearly, the exact values of α are poorly constrained given the large uncertainties, especially in the mass of the sources. As long as M ≳ 2 M_⊙, however, α ≲ 2.

Important information on the evolutionary stage of a source can be obtained from the observations of its outflow. The conclusions drawn from spectral observations (e.g., Hatchell et al. 2007b; Öberg et al. 2010), however, have been affected by the lack of an adequate knowledge of the spatial distribution of the sources. We recall that Hirano et al. (1999) detected only B1-bS/bN, while in the Bolocam survey at 1.1 mm (Enoch et al. 2006) only Bolo 81 was detected. We also recall that the Spitzer maps (Rebull et al. 2007) show only S295, while in the SCUBA map at 850 μm B1-bS and B1-bN are blended. Thanks to Herschel data, we can clearly see B1-bS/bN and S295 for the first time in the same map, in particular in the PACS bands. This new definition will allow a better understanding of the spectroscopic observations. For instance, by comparing our Fig. 2 with the CH₃OH maps by Öberg et al. (2010), it appears possible that B1-bS coincides with the peak in the CH₃OH map, in between the two identified outflows (the authors give the coordinates of only one outflow, which we reported in Fig. 2).

The outflows render it likely that one of the two sources is a local density enhancement (knot). Indeed, Gueth et al. (2003) concluded that emission knots generated by the shocked outflow in the vicinity of a protostellar object can be misinterpreted, e.g., as starless clumps. This conclusion, however, has been derived from mm and submm maps; to what extent it also holds in the FIR is unknown.

We finally note that the number of FHSC candidates found in Perseus is quite high. Pineda et al. (2011) derived that the number of expected FHSCs in this region should be  ≤ 0.2, 30 times fewer than the six proposed candidates reported in Table 5. In addition to the obvious argument that eventually none of the six candidates may be confirmed to be an FHSC, there are other ways to explain this discrepancy. The expected number of FHSCs is found by assuming that the ratio of the number of FHSCs over the number of Class 0 sources is equal to the ratio between their respective lifetimes. Since the two ratios are equal only for continuous star formation, a first hypothesis is that we are observing a burst of star formation in Perseus. The other possibility, which like the previous one has already been suggested by Pineda et al. (2011), is that we assume incorrect lifetimes for the FSHC phase, the Class 0 phase, or both. Moreover the number of Class 0 sources may be underestimated (Schnee et al. 2012), which would also impact the estimated lifetime. Though the catalogue of detected sources in Perseus is not yet released (Pezzuto et al., in prep.), it is reasonable to suppose that Herschel observations will update the number of starless cores, pre-stellar cores, and Class 0 sources for all nearby star-forming regions. With this information, it is possible that the expected number of FHSCs will change.

Finally, we note that our sources are quite close to each other, with a projected separation of 20′′. If the physical distance is not very different from the projected one, then B1-bS and B1-bN are  ≳4700 AU apart, corresponding to a few Jeans lengths at 10 K. It is then possible that these sources formed more or less at the same time from the fragmentation of a larger structure. This possibility, if confirmed, would explain why we found two FHSC candidates in close proximity.

5. Conclusions

We have presented the results of fitting the SEDs of two sources in the Perseus star-forming region. Data were derived from Herschel photometry observations collected within the Gould Belt Survey programme (André et al. 2010). The two sources are B1-bS and B1-bN, which were discovered by Hirano et al. (1999) from interferometric observations at millimetre wavelengths. B1-bS is the only source in the western part of the star-forming region in Perseus, detected in all six Herschel bands, and not visible in the Spitzer 24 μm map. These two criteria are important for the SED-based detection of FHSC candidates. The SED alone is not sufficient to determine the evolutionary status of a source, and other proposed FHSC candidates have been observed with interferometry to better assess their status. But considering only to the photometric criteria, B1-bS is the only candidate that satisfies the detection at 70 μm and the non detection at 24 μm.

Herschel data were complemented with Spitzer, JCMT-SCUBA, and CSO-Bolocam data. The resulting SED was fitted with a two-component model that describes the emission in terms of a blackbody, which roughly corresponds to the compact central object, and a greybody. i.e., the surrounding dusty envelope. We also tried to model the SED with a simple greybody alone, which is adequate for a pre-stellar core, and with a proper radiative transfer model (Robitaille et al. 2007), suited for more evolved sources. Both latter models failed to reproduce the SED over the whole observed spectral range. We conclude that B1-bS shows almost all expected characteristics expected in an FHSC. Only its luminosity is too high with respect to the theoretical predictions. Indeed, the high luminosity is at present the strongest argument against our candidate being an FHSC.

The SED of B1-bN is similar to that of B1-bS, with the important difference that the former is not detected at 70 μm. For the rest, this source seems similar to B1-bS. The two sources are situated a few 10³ AU apart, corresponding to a few Jeans lengths. It is then possible that these two sources formed at almost the same time from the fragmentation of a larger structure.

Additional observations at long wavelengths with high spatial resolution, such as those that ALMA will undertake, are clearly needed to better understand the nature of B1-bS and B1-bN. It is reasonable, however, to conclude that Herschel has observed two young objects that could be in a very early stage of protostellar formation, and we propose them as two new FHSC candidates.

Online material

Appendix A: SED fitting and photometry

For our two-component models, we modelled each set of data as the sum of a blackbody and a greybody envelope: F_ν = B(T_b)e^−τΩ_b + G(T_g,λ₀,β)Ω_g, where $\mathit{G}\mathrm{\left(}{\mathit{T}}_{\mathrm{g}}\mathit{,}{\mathit{\lambda }}_{\mathrm{0}}\mathit{,\beta }\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}\mathrm{)}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{g}}\mathrm{\right)}{\mathrm{\Omega }}_{\mathrm{g}}\mathit{,}$(A.1)with τ, the optical depth, parametrized as a power-law, i.e., $\mathit{\tau }\mathrm{=}{\left(\frac{{\mathit{\lambda }}_{\mathrm{0}}}{\mathit{\lambda }}\right)}^{\mathit{\beta }}$(A.2)where λ₀ is the wavelength at which τ = 1, and Ω_g and Ω_b are the solid angles of each respective component. The independent, unknown parameters are then T_b,T_g,λ₀, and β, while the solid angles are computed by scaling the model fluxes to the data. This model describes our sources in terms of two components: a central source, the blackbody, embedded in a dusty envelope whose emission is modelled with a greybody. Since we do not impose that the envelope is optically thin at all wavelengths, the blackbody radiation is attenuated by a factor that describes the absorption caused by the dust opacity. The absorption term, strictly speaking, implies that the envelope cannot be isothermal, therefore, as explained in Sect. 3, the set of parameters T_g,λ₀,β,Ω_g should be considered only as describing the average physical conditions of the envelope. The density profile in the cores usually follows a power law, ρ ∝ r^−q, with typically q ≲ 2, so that the mass, and then the flux, increases with the radius. T_g should then be indicative of the temperature in the outer region of the envelope.

The best-fit model was found inside a grid prepared by varying the parameters in the following intervals: 5 ≤ T(K) ≤ 50 in steps of 1 K; 0 ≤ β ≤ 5 in steps of 0.5; and 10 ≤ λ₀(μm) ≤ 600 in steps of 10 μm². These long intervals, longer than physically plausible, were chosen to ensure that the best-fit parameters would not fall on the border of the grid. The best fit was then chosen as the one with the smallest χ². During the search, we applied a few constraints: a) T_b > T_g, i.e., the greybody envelope must be colder than the blackbody component; b) Ω_b < Ω_g, the blackbody must be smaller than the greybody; c) the model must predict a 24 μm flux lower than the upper limit; d) the predicted flux at 1.1 mm must be lower than the measured value; and e) the envelope mass must be in the range of 0.1–10 M_⊙.

Once the minimum χ² was found, the 1σ uncertainties in the parameters were found by considering all models with a χ² in the range ${\mathit{\chi }}_{\min }^{\mathrm{2}}\mathrm{\le }{\mathit{\chi }}^{\mathrm{2}}\mathrm{\le }{\mathit{\chi }}_{\min }^{\mathrm{2}}\mathrm{+}\mathrm{1}$ (Andrae 2010). For B1-bS, only eight models were found in this χ² range. For B1-bN, not having a 70 μm flux had the consequence that 242 models, out of a total of 27646, had ${\mathit{\chi }}^{\mathrm{2}}\mathit{<}{\mathit{\chi }}_{\min }^{\mathrm{2}}\mathrm{+}\mathrm{1}$. Among these 242 models, the best fit was the one reported in Table 3, but with very large uncertainties. To decrease the uncertainties, we added another constraint and selected only the 107 models with λ₀ ≤ 200 μm, in analogy with the results of the selection for B1-bS. Finally, assuming that the dust properties are the same in the two sources, we considered only the models with β = 2, discarding 70 models with β = 1.5.

For the Bonnor-Ebert mass, we used the formula reported in Könyves et al. (2010)${\mathit{m}}_{\mathrm{BE}}\mathrm{=}\mathrm{2.4}\frac{\mathit{R}{\mathit{a}}^{\mathrm{2}}}{\mathit{G}}\mathit{,}$(A.3)where a is the sound speed corresponding to the temperature found from the fit, assuming a molecular weight μ = 2.33μ_H = 3.90 × 10^-24   g. For R we used the radius found from the fit, i.e., $\mathit{R}\mathrm{=}\mathit{D}\sqrt{{\mathrm{\Omega }}_{\mathrm{g}}\mathit{/}\mathit{\pi }}$.

A.1. Photometry corrections

Before comparing the model fluxes with the observed data, two steps were performed. For the first step, a colour correction was made; the flux calibration of PACS and SPIRE was performed under the usual assumption that the SED of a source displays a flat νF_ν spectrum. For all other kinds of SED, the derived fluxes must be colour-corrected according to the intrinsic source spectrum. To correct the fluxes properly, however, we need to know the spectrum a priori, but that is what we aim to derive from the data themselves. To overcome this circular problem, we computed a correction the other way round: a set of greybody models were computed over a wide frequency range and the fluxes were derived as ${\mathit{F}}_{\mathrm{cc}}\mathrm{=}\mathrm{\int }{\mathit{F}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T,}{\mathit{\lambda }}_{\mathrm{0}}\mathit{,\beta }\mathrm{\right)}{\mathrm{RSRF}}_{\mathit{\nu }}\mathrm{d}\mathit{\nu ,}$(A.4)where RSRFν stands for the relative spectral response function of each of the PACS or SPIRE filters, computed by the instruments control teams. To obtain the corresponding flux density, F_cc was divided by the appropriate filter width that we computed by imposing the condition that the colour corrections at the effective wavelengths is 1 for an SED of constant νF_ν. The derived Δλ and Δν are reported in Table A.1. As a consistency check, we compared the colour corrections we obtained for a blackbody with the values found by the instrument teams³. For PACS, the agreement is better than 1% in all the bands for T ≥ 7 K. The agreement at lower temperatures is the same at 100 μm and 160 μm, while it starts to diverge at 70 μm, in particular our correction for T = 6 K is 1.3% higher, at T = 5 K is 37% higher. At such extremly low temperatures, however, the colour corrections are quite difficult to be estimated since the intrinsic spectrum differs notably from the calibration spectrum. For SPIRE the comparison is less obvious because the colour corrections for a blackbody are not tabulated as a function of the temperature, but only in the limiting case of a blackbody in the Rayleigh-Jeans regime. For a blackbody at 100 K, we found that our colour correction is less than 1.3% higher at 250 μm, and less than 1% for the other two bands.

The colour corrections for the attenuated blackbody were found by inverting the definition of the greybody, i.e., B_ν(T)e^−τ = B_ν(T) − G_ν(T), from which (B_ν(T)e^−τ)_cc = B_cc(T) − G_cc(T). In this way, the same grid of models could be used.

Once the best fit was found, we computed the true colour correction factors, i.e., the amount by which the observed fluxes must be multiplied to derive the instrumental corrected fluxes. For the best-fit models reported in Table 3, the colour corrections f_cc are given in Table A.2.

For the second step we made before comparing model fluxes with observed fluxes, we took into account the Gaussian fit used to derive the photometric flux. For example, it is known that the PACS PSF is not a Gaussian, therefore an error is introduced in the photometry because of the Gaussian fit. To derive these errors, we reduced a set of PACS observations of flux calibrators and the results of aperture photometry and of synthetic photometry (Gaussian fitting) were compared. This exercise was repeated for a set of isolated compact sources in the Perseus field with different integrated fluxes. For 70 μm, 100 μm, and 160 μm, we found correction factors of 1.6, 1.5, and 1.4, respectively. We are making additional tests to improve our knowledge of Gaussian fit errors. For SPIRE, the PSF is much closer to a Gaussian profile, therefore we did not need correction factors for SPIRE fluxes.

In summary, Table 2 reports the measured fluxes multiplied by the factors that correct for the Gaussian fit. In Fig. 3, the Herschel data are the fluxes from Table 2 multiplied by the colour-correction. For instance, the flux at 70 μm of B1-bS resulting from the Gaussian fit is 0.138 Jy, which becomes 0.22 Jy after the multiplication by 1.6. This is the value reported in Table 2. After being multiplied by 0.79, see Table A.2, the flux becomes 0.17 Jy, which is the value used in Fig. 3.

We do not yet have estimates of the uncertainties in the maps, therefore it is not possible to give reliable uncertainties to each flux. The uncertainties given in Table 2 were found after a comparison of CuTEx fluxes with those found with the getsources algorithm (Men’shchikov et al. 2012).

A.2. Photometry in SCUBA maps

The effective beam of SCUBA at 450 μm is a Gaussian of FWHM 17$\stackrel{\mathrm{″}}{.}$3 (Di Francesco et al. 2008), but the actual beam consists of two different spatial components: an inner one with an FWHM of 11′′, and an outer one with an FWHM of 40′′. The sizes we derived for B1-bS and B1-bN are smaller than 17$\stackrel{\mathrm{″}}{.}$3. Since CuTEx is more sensitive to compact sources than to extended sources, we conclude that what CuTEx fitted was just the inner component of the sources, while the extended component, larger than the sources’ separation, was seen by CuTEx as background. To estimate the total flux of the individual sources, we used the following approach. If a source has an intrinsic size θ, then ${\mathit{F}}_{\mathrm{tot}}\mathrm{=}{\mathit{P}}_{\mathrm{1}}{\mathrm{(}{\mathit{\theta }}_{\mathrm{1}}\mathit{/}{\mathrm{11}}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathrm{2}}{\mathrm{(}{\mathit{\theta }}_{\mathrm{4}}\mathit{/}{\mathrm{40}}^{\mathrm{)}}}^{\mathrm{2}}\mathit{,}$where ${\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\sqrt{{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}{\mathrm{11}}^{\mathrm{2}}}$ and similarly for θ₄. From the analysis of Di Francesco et al. (2008), we know that the peak fluxes are P₁ = 0.88P and P₂ = 0.12P, where P is the peak flux of the 17$\stackrel{\mathrm{″}}{.}$3 FWHM Gaussian, from which P₂ = P₁ × 0.12/0.88. Finally, expressing θ₄ as a function of θ₁, we have ${\mathit{F}}_{\mathrm{tot}}\mathrm{=}{\mathit{P}}_{\mathrm{1}}{\mathrm{(}{\mathit{\theta }}_{\mathrm{1}}\mathit{/}{\mathrm{11}}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{0.12}}{\mathrm{0.88}}{\mathit{P}}_{\mathrm{1}}\left(\frac{{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}{\mathrm{11}}^{\mathrm{2}}\mathrm{+}{\mathrm{40}}^{\mathrm{2}}}{{\mathrm{40}}^{\mathrm{2}}}\right)\mathrm{·}$(A.5)At 850 μm, CuTEx finds sizes that are larger than the effective beam of 22$\stackrel{\mathrm{″}}{.}$9, but since this beam is again the combination of two Gaussians (Di Francesco et al. 2008) of 19$\stackrel{\mathrm{″}}{.}$5 and 40′′FWHM, it is likely that also in this case CuTEx fits just the inner Gaussian and sees the second one as a background. The formula we used to find the total flux is the same as Eq. (A.5), with 19$\stackrel{\mathrm{″}}{.}$5 instead of 11′′. Since our sources are slightly larger than the effective beam, however, the applicability of Eq. (A.5) to the 850 μm flux is less robust.

A.3. The mass of a greybody

Without Herschel data to define the SED in the FIR, the emitting mass is usually derived from the flux measured at long wavelengths where the envelope becomes optically thin. In this limit Eq. (A.1) then becomes ${\mathit{F}}_{\mathit{\nu }}\mathrm{\approx }\mathit{\tau }{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{\Omega }\mathit{.}$(A.6)If the envelope is optically thin, we see the whole mass distribution M, so that $\mathit{\tau }\mathrm{=}{\mathit{\kappa }}_{\mathit{\nu }}\mathrm{\int }\mathit{\rho }\mathrm{d}\mathit{s}\mathrm{\approx }{\mathit{\kappa }}_{\mathrm{ref}}{\left(\frac{{\mathit{\lambda }}_{\mathrm{ref}}}{\mathit{\lambda }}\right)}^{\mathit{\beta }}\frac{\mathit{M}}{\mathit{A}}\mathit{,}$(A.7)where A is the projected area and κ_ref is the opacity at the reference wavelength λ_ref. As in other papers of the Gould Belt consortium, e.g., Könyves et al. (2010), we adopted an opacity of κ_ref = 0.1   cm²   g^-1 at λ_ref = 300 μm (Beckwith et al. 1990). With this choice, the GBS opacity law is very similar to the one by Hildebrand (1983), who adopted κ_ref = 0.1   cm²   g^-1 at λ_ref = 250 μm and β = 2 for λ ≥ 250 μm. For this work, only κ_ref is important because β was derived directly from the fit. λ_ref was shifted to 300 μm for consistency with older results obtained with ground-based facilities, such as IRAM or JCMT. By comparing the column density map of IC 5146 obtained from Herschel observations and the GBS dust opacity law with that obtained from SCUBA (sub)mm data and with a different opacity law, Arzoumanian et al. (2011) found that the GBS dust opacity law works well for the Herschel data. The exact value for κ_ref is in any case uncertain and likely dependent on the dust temperature. An uncertainty in the derived mass of a factor  ~2 cannot be excluded.

Substituting Eq. (A.7) into (A.6) and writing the solid angle as A/D² = πR²/D², with R equal to the outer radius of the envelope, assumed spherical, and D equal to the distance to the source, we arrive at the well-known formula $\mathit{M}\mathrm{\approx }\frac{{\mathit{F}}_{\mathit{\nu }}{\mathit{D}}^{\mathrm{2}}}{{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\mathrm{·}$(A.8)Since from our best-fit models λ₀ and β are also found, we derived a different formula to estimate the mass of a greybody envelope.

At wavelengths where the envelope is optically thin, we find by combining Eqs. (A.2) and (A.7) that ${\left(\frac{{\mathit{\lambda }}_{\mathrm{0}}}{\mathit{\lambda }}\right)}^{\mathit{\beta }}\mathrm{\approx }{\mathit{\kappa }}_{\mathrm{ref}}{\left(\frac{{\mathit{\lambda }}_{\mathrm{ref}}}{\mathit{\lambda }}\right)}^{\mathit{\beta }}\frac{\mathit{M}}{\mathit{A}}\mathit{,}$from which $\mathit{M}\mathrm{=}\frac{{\mathit{D}}^{\mathrm{2}}\mathrm{\Omega }}{{\mathit{\kappa }}_{\mathrm{ref}}}{\left(\frac{{\mathit{\lambda }}_{\mathrm{0}}}{{\mathit{\lambda }}_{\mathrm{ref}}}\right)}^{\mathit{\beta }}\mathrm{·}$(A.9)This new equation, equivalent for a greybody to the unnumbered equation on page 274 in Hildebrand (1983), gives the mass inside a sphere of radius R that, for the given opacity law, has optical depth 1 at λ₀. Of course, Eq. (A.9) can be applied only if we know β and λ₀, but these are what Herschel data allow us to find. The mass of the envelopes, discussed in Sect. 3, were found using the above formula where the greybody parameters are those derived from the fit that are reported in Table 3.

Appendix B: Sigma clipping with a threshold dependent on the sample size

Before generating the final maps, the bolometers’ time series have to be corrected for glitches caused by the impact of high-energy particles falling onto the detectors. To achieve this correction, we exploited the spatial redundancy provided by the telescope movement, which ensures that each sky pixel of the final map has been observed several times with different bolometers. Outlier detection was then made with a sigma-clipping algorithm. Namely, given a sample of N values, estimates for the mean and standard deviations are derived. All values that differ from the mean by more than α standard deviations are considered outliers and removed from the sample.

One problem with the sigma-clipping algorithm is the choice of the number of standard deviations above which a datum is considered an outlier. To avoid making an arbitrary choice of α, a formula has been derived that makes α dependent on the size of the sample. Below, we give some details of the method.

For a Gaussian distribution with mean m and standard deviation σ, the probability to derive a value x such that m − ασ < x < m + ασ, is given by erf($\mathit{\alpha }\mathit{/}\sqrt{\mathrm{2}}\mathrm{)}$ where erf(y) is the error function $\mathrm{erf}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}}{\sqrt{\mathit{\pi }}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{y}}\mathrm{e}\mathrm{-}{\mathit{t}}^{\mathrm{2}}\mathrm{d}\mathit{t}\mathit{.}$For example, if α = 3, the probability is erf($\mathrm{3}\mathit{/}\sqrt{\mathrm{2}}\mathrm{)}\mathrm{=}\mathrm{0.9973}$. Obviously, 1 − erf$\mathrm{(}\mathit{\alpha }\mathit{/}\sqrt{\mathrm{2}}\mathrm{)}$ gives the probability to derive a value outside the same interval. This probability tends to zero when α → ∞ (by definition erf(+ ∞) = 1).

In a sample of N data, the probability given by the error function, multiplied by N, can be interpreted as a the expected number of points inside or outside a certain interval around the mean. For any value of α, the number p of points that differ more than ασ from the mean is then $\mathit{p}\mathrm{\left(}\mathit{\alpha }\mathrm{;}\mathit{N}\mathrm{\right)}\mathrm{=}\left[\mathrm{1}\mathrm{-}\mathrm{erf}\left(\frac{\mathit{\alpha }}{\sqrt{\mathrm{2}}}\right)\right]\mathrm{·}\mathit{N}\mathit{.}$For an ideal Gaussian, any value of α is allowed. In a real experiment, however, p is an integer number and data differing from the mean by ασ, such that p(α;N) ≪ 1, are suspicious. To implement this condition in the sigma-clipping algorithm, we have to define a precise value of α, say $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{\alpha }\end{array}$, above which observed data are considered outliers and removed from the sample. To this aim, we define $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{\alpha }\end{array}$ as $\mathit{p}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{\alpha }\end{array}\mathrm{;}\mathit{N}\mathrm{\right)}\mathrm{=}\mathrm{1}$or $\mathrm{erf}\left(\frac{\begin{array}{}\mathrm{˜}\\ \mathit{\alpha }\end{array}}{\sqrt{\mathrm{2}}}\right)\mathrm{=}\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{\mathit{N}}\mathrm{·}$(B.1)All the points, if any, that are outside the interval $\mathit{m}\mathrm{-}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathit{\sigma }\mathrm{\le }\mathit{x}\mathrm{\le }\mathit{m}\mathrm{+}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathit{\sigma }$ are considered outliers and removed from the sample. This is equivalent to assume that $\mathit{p}\mathrm{(}\mathit{\alpha }\mathit{>}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathrm{;}\mathit{N}\mathrm{)}\mathrm{\equiv }\mathrm{0}$ instead of $\mathit{p}\mathrm{(}\mathit{\alpha }\mathit{>}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathrm{;}\mathit{N}\mathrm{)}\mathit{<}\mathrm{1}$.

When N → ∞ then $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathrm{\to }\mathrm{\infty }$ too. The mathematical property of the Gaussian distribution that for an infinitely large sample any value is allowed is preserved.

It is not possible to analytically invert Eq. (B.1) for a given N, but it is possible to derive an approximate analytical expression. Indeed, we found that in the range $\mathrm{1}\mathrm{\le }\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathrm{\le }\mathrm{5}$ Eq. (B.1), written in terms of log (N), can be approximated with the parabola $\log \mathrm{\left(}\mathit{N}\mathrm{\right)}\mathrm{=}\mathrm{0.0794}\mathrm{+}\mathrm{0.2282}\begin{array}{}\mathrm{˜}\\ \mathit{\alpha }\end{array}\mathrm{+}\mathrm{0.2005}\stackrel{}{{{\mathit{\alpha }}^{\mathrm{2}}}_{\mathrm{˜}}}\mathit{.}$(B.2)In Fig. B.1 we show the ratio between N as given by Eqs. (B.1) and (B.2). The ratio is always lower than 3%.

Fig. B.1 The ratio between $\mathit{N}\mathrm{\equiv }\mathit{N}\mathrm{\left(}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathrm{\right)}$ as given by Eq. (B.1), and $\mathit{N}\mathrm{\equiv }\mathit{N}\mathrm{\left(}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathrm{\right)}$ approximated with the parabola given in Eq. (B.2), as a function of the number of standard deviations. The difference between the two curves is always less than 3%. Fitting a cubic instead of a parabola does not improve the approximation significantly and makes the relation $\mathit{N}\mathrm{\equiv }\mathit{N}\mathrm{\left(}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathrm{\right)}$ harder to invert.

Inverting Eq. (B.2) is trivial, the result is $\begin{array}{}\mathrm{˜}\\ \mathit{\alpha }\end{array}\mathrm{=}\mathrm{-}\mathrm{0.569}\mathrm{+}\sqrt{\mathrm{-}\mathrm{0.072}\mathrm{+}\mathrm{4.99}\log \mathrm{\left(}\mathit{N}\mathrm{\right)}}\mathit{.}$(B.3)For a given N, the value of $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{\alpha }\end{array}$ is found from the above equation. The number of expected points distant from the mean by more than $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathit{\sigma }$ is zero. If one or more outliers are found instead, they are removed. New values of m and σ are recomputed and the new N′ < N is used in Eq. (B.3). This process is repeated until no other outliers are present in the sample. The procedure may not converge, for instance, if the noise is not Gaussian. For this reason, we did not iterate the formula. Instead, for each point of the map the detection of outliers was performed only once. There can be more than 1 outliers found at the first iteration, however.

In the central region of the map, where the coverage (and thus N) is high, the threshold for clipping is higher than in the outskirts of the map where the coverage is low. For instance, if a sky pixel has been observed with 40 bolometers, the above formula gives $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathrm{=}\mathrm{2.25}$. Accordingly, once we have estimated the mean m and the standard deviations σ, all values x_i such that ABS(x_i − m) > 2.25σ are flagged as outliers. If instead a pixel has been observed with 20 bolometers, the threshold is lowered to 1.96σ.

Acknowledgments

We thank the anonymous referee for the many comments and suggestions that made the paper better by far than the original. We also thank B. Commerçon for helping us to answer one of the points raised by the referee; P. Martin and J. Kirk for helpful clarifications on the dust opacity; and A. Men’shchikov for running getsources on our map to independently verify our extraction. S.P. spent one week at the Herzberg Institute of Astrophysics in Victoria, during which this work was finalized; he thanks all people at the institute for the kind and warm hospitality. K.L.J.R., D.P., M.P., D.E., and E.S. are funded under ASI contracts I/038/08/0 and I/005/11/0. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (UK) and including: Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). HIPE is a joint development by the Herschel SGSC, consisting of ESA, the NASA HSC, and the HIFI, PACS and SPIRE consortia.

References

All Tables

All Figures

Fig. 1 Maps of 1$\stackrel{\mathrm{\prime }}{.}$5  ×  1$\stackrel{\mathrm{\prime }}{.}$5 of B1-bS and B1-bN at all wavelengths, with sizes and positions derived from the Gaussian fits. The FWHM of the respective instrumental beam is shown as a circle in the bottom-left corner of each respective panel. At λ ≤ 100 μm, the source S295 is also visible. Colour bars are in Jy/beam for the SCUBA maps, and MJy/sr in all the others. In the Spitzer images, positions and sizes of the sources are copied from the 70 μm map. IRAC maps are in the top row, from left to right: a) 3.8 μm; b) 4.5 μm; c) 5.8 μm; d) 8.0 μm. Second row: e) Spitzer 24 μm; f) PACS 70 μm, B1-bN is not really detected and its flux corresponds to a 1σ rms point source; g) PACS 100 μm; h) PACS 160 μm, S295 is no longer visible. Third row: i) SPIRE 250 μm; j) SPIRE 350 μm; k) SCUBA 450 μm, where the beam shown is 11′′, see Appendix A.2; l) SPIRE 500 μm. Bottom row: m) SCUBA 850 μm; n) Bolocam 1.1 mm, where only one source has been fitted.

In the text

Fig. 2 Positions (J2000) of the sources at all wavelengths; see text for references. Two black open circles are centred on the position of B1-bS and B1-bN after applying a positional offset correction found by comparing the PACS and Spitzer positions of S295, see text for details. The radius of the circles is the 1σ of the mean of the PACS coordinates, for B1-bS, and of the mean of PACS 100 μm, 160 μm, and SPIRE 250 μm coordinates for B1-bN. The positions at 350 μm, 500 μm, and 850 μm are very uncertain, because B1-bS and B1-bN are not resolved at those wavelengths. The blue open circle centred on Bolo 81 has a radius of 7′′, the pointing accuracy of Bolocam observations (Enoch et al. 2006).

In the text

Fig. 3 Best-fit two-component models for B1-bS and B1-bN. Filled circles are Herschel data (for B1-bN the 70 μm flux is an upper limit), open circles are SCUBA data, and upper limits at λ ≤ 24 μm are from Spitzer IRAC and MIPS. The solid line is the best-fit model, a sum of a blackbody (dashed line dominating at shorter wavelengths) plus a greybody (dashed line more significant at longer wavelengths). Herschel fluxes are those reported in Table 2 multiplied by the colour correction factors, see Appendix A.1.

In the text

Fig. B.1 The ratio between $\mathit{N}\mathrm{\equiv }\mathit{N}\mathrm{\left(}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathrm{\right)}$ as given by Eq. (B.1), and $\mathit{N}\mathrm{\equiv }\mathit{N}\mathrm{\left(}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathrm{\right)}$ approximated with the parabola given in Eq. (B.2), as a function of the number of standard deviations. The difference between the two curves is always less than 3%. Fitting a cubic instead of a parabola does not improve the approximation significantly and makes the relation $\mathit{N}\mathrm{\equiv }\mathit{N}\mathrm{\left(}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\alpha }\end{array}\mathrm{\right)}$ harder to invert.

In the text