THE HERSCHEL ORION PROTOSTAR SURVEY: SPECTRAL ENERGY DISTRIBUTIONS AND FITS USING A GRID OF PROTOSTELLAR MODELS

In order to construct SEDs for our sample of 330 YSOs, we combined our own observations with data from the literature and existing catalogs. For the near-infrared photometry, we used J, H, and K_s data from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006). For the mid-infrared spectral region, we used Spitzer data from Kryukova et al. (2012) and Megeath et al. (2012): the Infrared Array Camera (IRAC; Fazio et al. 2004) provided 3.6, 4.5, 5.8, and 8.0 μm photometry, while the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004) provided 24 μm photometry. In addition, most of the YSOs in the HOPS sample were also observed with the Infrared Spectrograph (IRS; Houck et al. 2004) on Spitzer using the Short-Low (SL; 5.2–14 μm) and Long-Low (LL; 14–38 μm) modules, both with a spectral resolution of about 90 (see, e.g., Kim et al. [2016] for a description of IRS data reduction). Herschel PACS data at 70, 100, and 160 μm yielded far-infrared photometric data points (B. Ali et al. 2016, in preparation; the 100 μm data are from the Gould Belt Survey; e.g., André et al. 2010). Most YSOs were also observed at 350 and 870 μm (see Stutz et al. 2013) by the APEX telescope using the SABOCA and LABOCA instruments (Siringo et al. 2009, 2010, respectively). Thus, our SEDs have well-sampled wavelength coverage from 1.2 to 870 μm; we did not include additional data from the literature in order to preserve a homogeneous data set for all the objects in our sample.

The aperture radius used for the photometry varies depending on the instrument and wave band. The photometry in the 2MASS catalog was derived from point-spread function (PSF) fits using data from 4'' apertures around each object (see the Explanatory Supplement to the 2MASS All Sky Data Release and Extended Mission Products). Megeath et al. (2012) used an aperture radius of 24 for IRAC and PSF photometry for MIPS 24 μm data. We used aperture radii of 96 and sky annuli of 96–192 for PACS 70 and 100 μm images; we then applied aperture correction factors of 0.7331 and 0.6944 to the 70 and 100 μm fluxes, respectively. For PACS 160 μm, we used an aperture radius of 128, a sky annulus of 128–256, and an aperture correction factor of 0.6602. In some cases (background contamination, close companions) we used PSF photometry at 70 and 160 μm instead (see B. Ali et al. 2016, in preparation, for details). Finally, we adopted beam fluxes at 350 and 870 μm (with FWHMs of 734 and 19'', respectively). The IRS SL module has a slit width of 36, while the LL module is wider, with a slit width of 105. Sometimes the flux level of the two segments did not match at 14 μm (owing to slight mispointings or more extended emission from surrounding material measured in LL), and in these cases usually the SL spectrum was scaled by at most a factor of ∼1.4 (typically 1.1–1.2). In a few cases, especially when the LL spectrum included substantial amounts of extended emission or flux from a nearby object, the LL spectrum was scaled down to match the flux level of the SL spectrum at 14 μm, typically by a factor of 0.8–0.9. We discuss how the different aperture sizes are accounted for in the model fluxes in Section 4.2.

The SEDs of our HOPS sample are shown in Figure 1 together with their best-fit models from our model grid (see sections below); the data are listed in Table 2. Many objects display a deep silicate absorption feature at 10 μm and ice features in the 5–8 μm region, as expected for protostars. Those objects with very deep 10 μm features and steeply rising SEDs are likely deeply embedded protostars, often seen at high inclination angles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (22 panels)

Tarball of figure set images

A large fraction (203 out of 330) of the young stars in our sample have at least one Spitzer-detected source within a radius of 15''; in most cases, this "companion" is faint in the infrared and likely a background star or galaxy. Thus, the emission at far-IR and submillimeter wavelengths is expected to be dominated by the protostar or pre-main-sequence star, and we can assume that the SEDs are representative of the YSOs even if the nearby sources cannot be separated at these wavelengths. There are a few YSOs that have objects separated by just 1''–3'' and are only resolved in one or two IRAC bands (HOPS 22, 78, 108, 184, 203, 247, 293, 364); in these cases we used the flux at the IRAC position that most closely matched those at longer wavelengths. We note that some of these very close "companions" are likely outflow knots. There are also unresolved binaries, which appear as single sources even in the IRAC observations (Kounkel et al. 2016); in these cases our SEDs show the combined flux in all wave bands. If two point sources are not fully resolved and the resulting blended source is elongated, no IRAC photometry was extracted. In such cases, a protostar may not have IRAC fluxes even though it was detected in the Spitzer images.

There are also several protostars that lie close to other protostars: HOPS 66 and 370 (d = 149), HOPS 76 and 78 (d = 141), HOPS 86 and 87 (d = 121), HOPS 117 and 118 (d = 137), HOPS 121 and 123 (d = 76), HOPS 124 and 125 (d = 98), HOPS 165 and 203 (d = 133), HOPS 175 and 176 (d = 80), HOPS 181 and 182 (d = 102), HOPS 225 and 226 (d = 92), HOPS 239 and 241 (d = 124), HOPS 262 and 263 (d = 63), HOPS 316 and 358 (d = 69), HOPS 332 and 390 (d = 112), HOPS 340 and 341 (d = 47), and HOPS 386 and 387 (d = 99). HOPS 105 lies 87 to the north of an infrared-bright source, identified by Megeath et al. (2012) as a young star with a protoplanetary disk. This source is brighter than HOPS 105 in all Spitzer bands and at 70 μm, but it is well separated at all wavelengths. A similar situation applies to HOPS 128, which has a disk-dominated source 63 to the southeast. HOPS 108 is 66 from HOPS 64, which is brighter than HOPS 108 out to 8 μm, but not detected in the far-IR and submillimeter. HOPS 108 also lies 166 from HOPS 369 and 282 from HOPS 370. HOPS 140 has two neighboring sources, at 96 and 139, that are likely surrounded by protoplanetary disks; they are both brighter than HOPS 140 out to 8 μm, but at 70 μm and beyond HOPS 140 dominates. HOPS 144 lies 79 from HOPS 377; there is also a somewhat fainter, red source 117 to the northeast, which is not detected beyond 24 μm. This source also lies 97 to the southwest of HOPS 143. HOPS 173 forms a small cluster with HOPS 174 (at 71) and HOPS 380 (at 114); HOPS 174 is the brightest source out to 24 μm, but at 70 μm HOPS 173 takes on this role. Also HOPS 322, 323, and 389 form a group of protostars; HOPS 322 lies 134 from HOPS 389 and 201 from HOPS 323, while HOPS 323 and 389 are 102 apart. HOPS 323 is the brightest source.

Thus, there are 45 targets in our sample that have an object within 15'' that is bright in the mid- or far-IR and that is resolved with IRAC and MIPS. Given that Megeath et al. (2012) used PSF photometry for the MIPS 24 μm observations, they obtained reliable fluxes even for companions separated by less than 6'', the typical PSF FWHM. For fluxes at 70 and 160 μm, we also used PSF photometry for objects that were point sources, but too close for aperture photometry. In cases where the fluxes could not be determined even with PSF photometry, we had to adopt upper limits instead. Similarly, we performed PSF photometry on protostars without companions, but contaminated by extended or filamentary emission; if the PSF photometry did not return a good fit, we used the flux value from aperture photometry as an upper limit.

Since most of our targets have an IRS spectrum, in addition to data points from IRAC at 5.8 and 8 μm and from MIPS at 24 μm, we can detect discrepancies if flux values at similar wavelengths, but from different instruments, do not agree. They might be due to calibration or extraction problems in the IRS spectrum (for example, some extended emission around the target or a close companion), but also to variability. We assumed the former scenario if the flux deviations between IRS and IRAC and between IRS and MIPS were similar (and more than 10%, a conservative estimate for the typical calibration uncertainty), and in such cases scaled the IRS spectrum to the MIPS 24 μm flux. Even though this scaling could mask actual variability, it creates a representative SED for the YSO and yields an estimate of the protostellar parameters from model fits of the SED.

In Appendix A we identify potentially variable HOPS targets based on their mid-IR fluxes and find that about 5% of the protostars with IRS, IRAC, and MIPS data could be variable. The Young Stellar Object Variability (YSOVAR) program, which monitored large samples of protostars and pre-main-sequence stars in nearby star-forming regions with Spitzer at 3.6 and 4.5 μm (Morales-Calderón et al. 2011; Cody et al. 2014; Günther et al. 2014; Rebull et al. 2014, 2015; Poppenhaeger et al. 2015; Wolk et al. 2015), found that up to ∼90% of flat-spectrum and Class I YSOs are variable on a timescale of days, with typical changes in brightness of 10%–20%. On longer timescales (years as opposed to days), 20%–40% of members of young clusters show long-term variability, with the highest fraction for those clusters with more Class I protostars (Rebull et al. 2014). In Orion, the fraction of variable Class I protostars is ∼85% (Morales-Calderón et al. 2011). Using a larger sample of protostars in Orion and IRAC data at 3.6, 4.5, 5.8, and 8.0 μm, Megeath et al. (2012) found that, on a timescale of about 6 months, 60%–70% of Orion protostars show brightness variabilities of ∼20%, with some as high as a factor of four. Thus, given that our SEDs consist of noncontemporaneous data sets, small flux discrepancies should be common, but we also expect some protostars with large mismatches.

One protostar with a large discrepancy between various data sets is HOPS 223. It is an outbursting protostar (also known as V2775 Ori; Caratti o Garatti et al. 2011), and for its SED we had 2MASS, IRAC, and MIPS data from the pre-outburst phase available, while the IRS spectrum, PACS, and APEX data are from the post-outburst period. Thus, its SED does not represent an actual state of the object, and the derived T_bol and L_bol values are unreliable. Pre- and post-outburst SEDs and model fits for this protostar can be found in Fischer et al. (2012). HOPS 223 is the only protostar with an SED affected by extreme variability. A few more protostars, HOPS 71, 132, 143, 228, 274, and 299, show notable discrepancy between the IRAC and IRS fluxes, and to a minor extent between MIPS 24 μm and IRS, and thus have somewhat unreliable SEDs and SED-derived parameters. HOPS 383, which was identified as an outbursting Class 0 protostar by Safron et al. (2015), does not appear variable in the SED presented here, since we adopted post-outburst IRAC 3.6 and 4.5 μm fluxes obtained by the YSOVAR program (Morales-Calderón et al. 2011; Rebull et al. 2014) to construct a representative post-outburst SED for this object.

The SEDs provide the means to determine L_bol, T_bol, and the 4.5–24 μm spectral indices for our sample of protostars. For measuring the near- to mid-IR SED slope ($n=d\mathrm{log}(\lambda {F}_{\lambda })/d\mathrm{log}(\lambda )$), we chose a spectral index between 4.5 and 24 μm to minimize the effect of extinction on the short-wavelength data point; also, the IRAC 4.5 μm fluxes for our HOPS targets are more complete than the IRAC 3.6 μm fluxes owing to the lower extinction at this wavelength. For calculating L_bol and T_bol, we used all available fluxes for each object, including the IRS spectrum, assumed a distance of 420 pc, and used trapezoidal summation; for T_bol, we applied the equation from Myers & Ladd (1993). Figure 2 shows the distribution of n_4.5–24, T_bol, and L_bol values for our targets. There is a peak in the distribution of spectral indices around 0, while the distribution of T_bol values is relatively uniform from about 30 to 800 K. The bolometric luminosities cover a wide range, with a broad peak around 1 L_⊙. The median L_bol, T_bol, and n_4.5–24 values are 1.1 L_⊙, 146 K, and 0.68, respectively.

Zoom In Zoom Out Reset image size

Our distribution of L_bol values is very similar to the observed luminosity function of Orion protostars presented in Kryukova et al. (2012); both distributions peak around 1 L_⊙ and include values from ∼0.02 L_⊙ up to several hundred L_⊙. Some differences between the two distributions are expected, given that Kryukova et al. (2012) only had Spitzer 3.6–24 μm data available and thus had to extrapolate L_bol from the measured near- to mid-infrared luminosity. The main difference is a somewhat larger number of protostars with L_bol ≲ 0.5 L_⊙ for the Kryukova et al. (2012) Orion sample; our median L_bol value amounts to 1.1 L_⊙, while their value is 0.8 L_⊙. However, with the contaminating sources removed from their sample (which tend to have lower luminosities; see Kryukova et al. 2012 for details), their median bolometric luminosity and our value match. Overall, given that Orion is considered a region of high-mass star formation, its luminosity function is similar to that of other regions where massive star forms (Kryukova et al. 2012, 2014), and it is different from that of low-mass star-forming regions such as Taurus and Ophiuchus (Kryukova et al. 2012; Dunham et al. 2013). Compared to the sample of 230 protostars in 18 different molecular clouds studied by Dunham et al. (2013), the observed (i.e., not extinction-corrected) L_bol values of those protostars span from 0.01 to 69 L_⊙, with a median value of 0.9 L_⊙. However, given that almost half the protostars in the Dunham et al. (2013) sample lack far-IR and submillimeter data, the true luminosities are likely higher, which would bring the median closer to the Orion value. Finally, we note that the distribution of observed T_bol values from Dunham et al. (2013) is similar to our distribution for Orion protostars; the median T_bol of their and our sample is 160 and 146 K, respectively, and the bulk of their protostars also has T_bol values between 30 and 1000 K, with a tail down to temperatures of ∼10 K and another tail up to T_bol = 2700 K.

To separate our targets into Class 0, Class I, Class II, and flat-spectrum sources, we used the 4.5–24 μm spectral index (n_4.5–24) and/or bolometric temperature (T_bol): Class 0 protostars have n_4.5–24 > 0.3 and T_bol < 70 K, Class I protostars have n_4.5–24 > 0.3 and T_bol > 70 K, flat-spectrum sources have −0.3 < n_4.5–24 < 0.3, and Class II pre-main-sequence stars have n_4.5–24 < −0.3. Based on this, we identify 92 targets as Class 0 protostars, 125 as Class I protostars, 102 as flat-spectrum sources, and 11 as Class II pre-main-sequence stars (see Table 1 and Figure 3). There are nine protostars with T_bol values between 66.5 and 73.5 K (which corresponds to a ±5% range around the Class 0−I boundary of 70 K); six of them have T_bol > 70 K (HOPS 1, 18, 186, 256, 322, 370), and the other three have T_bol values just below 70 K (HOPS 75, 250, 361). These protostars' classification is less firm than for the other HOPS targets. There are also a few flat-spectrum sources whose classification is more uncertain: HOPS 45, 183, 192, 194, 210, 264, and 281 should be Class I protostars based on their 4.5–24 μm spectral index, but when considering the IRS spectrum (specifically, the 5–25 μm spectral index), they fall into the flat-spectrum regime (n_5–25 < 0.3). Also, for HOPS 45 and 194 the T_bol values are relatively high (>500 K). Similarly, HOPS 33, 134, 242, 255, and 284 should be Class II pre-main-sequence stars based on their 4.5–24 μm spectral index, but the spectral slope over the IRS wavelength range suggests that they are flat-spectrum sources. In these cases where the n_4.5–24 and n_5–25 spectral indices were somewhat discrepant, we adopted the latter, and thus these objects were classified as flat-spectrum sources.

Zoom In Zoom Out Reset image size

There are five objects with T_bol < 70 K and n_4.5–24 < 0 (HOPS 164, 340, 341, 373, 405); despite their negative 4.5–24 μm SED slopes, their SEDs either show or imply a deep silicate absorption feature at 10 μm, rise steeply in the mid- to far-IR, and their long-wavelength emission is strong. Thus, their T_bol values are low, and we identify them as Class 0 protostars, even though they have 4.5–24 μm spectral indices not typical of embedded protostars. In particular, HOPS 341, 373, and 405 are likely young protostars with dense envelopes (Stutz et al. 2013; see also Section 7.2.1). In the case of HOPS 373, the 4.5 μm flux may be contaminated by bright H₂ emission from an outflow shock, rendering the n_4.5–24 value more unreliable. This might also explain the negative 4.5–24 μm spectral index for the other four protostars.

Finally, the few Class II objects in our sample were thought to be potential protostars prior to their observations with Herschel. Their 4.5–24 μm SED slopes are usually just slightly more negative than the cutoff for a flat-spectrum source (−0.3); three Class II pre-main-sequence stars (HOPS 22, 184, 201) have SEDs that are typical of disks with inner holes, displaying a 10 μm silicate emission feature and a rising SED from 12 to about 20 μm (e.g., Kim et al. 2013). The SEDs of the other Class II objects are similar to those of flat-spectrum sources; thus, they could have (remnant) envelopes that contribute to their long-wavelength emission.

Our HOPS sample is mostly complete in the number of Class 0, Class I, and flat-spectrum sources in the areas of Orion surveyed by Spitzer, excluding the Orion Nebula (see Megeath et al. 2012; Stutz et al. 2013). Of the 357 unique YSOs originally identified in Spitzer data that were included in the HOPS sample and observed with PACS, 322 were detected at least at 70 μm, which amounts to a fraction of 90%. We removed likely contaminants and added 16 new protostars discovered in PACS data to get to our sample of 330 YSOs, most of which are protostars. Our lowest L_bol source is HOPS 208, with L_bol = 0.017 L_⊙. This protostar also has the lowest PACS 70 μm flux in our sample (8.2 mJy). Overall, our sample has 27 protostars with L_bol < 0.1 L_⊙, which places them in the luminosity range of very low luminosity objects (VeLLOs; di Francesco et al. 2007; Dunham et al. 2008). The number of VeLLOs in our sample is likely larger, given that VeLLOs are defined as having internal luminosities less than 0.1 L_⊙, and the bolometric luminosity has contributions from both the internal luminosity and that due to external heating (see Dunham et al. 2008). In addition, our sample could miss fainter flat-spectrum sources and Class 0 and Class I protostars. In fact, there are several faint YSOs without PACS data that were excluded from our sample, but do have Spitzer detections (see Appendix C.2.1).

There are 3040 models in the grid; they cover eight values for the total (i.e., intrinsic) luminosity, four disk radii, 19 envelope infall rates (which correspond to envelope densities), and five cavity opening angles. Each model is calculated for 10 different inclination angles, from 182 to 872, in equal steps in cos(i) (starting at 0.95 and ending at 0.05), resulting in 30,400 different model SEDs. The values for the various model parameters are listed in Table 3. Since there are a large number of parameters that can be set in the Whitney radiative transfer models, we focused on varying those parameters that affect the SED the most, leaving the other parameters at some typical values. For example, we assumed a stellar mass of 0.5 M_⊙, a disk mass of 0.05 M_⊙, and an envelope outer radius of 10,000 au. The stellar mass enters the code in two ways. First, it is needed to relate the density of the envelope to the infall rate (see Equation (1) below). Since we fit the density of the envelope, the infall rate plays no role in the best-fit envelope parameters; any stellar mass can be chosen to determine the infall rate for a given best-fit envelope density. Second, the stellar mass is combined with the stellar radius and disk accretion rate to set the disk accretion luminosity. Given that the accretion luminosity is the actual parameter that influences the SED, it does not matter which of the three factors is varied. For simplicity and reasons described below, we varied the disk accretion rate and the stellar radius, but left the stellar mass constant, to achieve different values for this component of the luminosity.

Table 3.  Model Parameters

Parameter	Description	Values	Units
Stellar Properties
M*	Stellar mass	0.5	M⊙
T*	Stellar effective temperature	4000	K
R*	Stellar radius	0.67, 2.09, 6.61	R⊙
Disk Properties
Mdisk	Disk mass	0.05	M⊙
Rdisk	Disk outer radius	5, 50, 100, 500	au
A	Radial exponent in disk density law	2.25	⋯
B	Vertical exponent in disk density law	1.25	⋯
${\dot{M}}_{\mathrm{disk},1}$	Disk-to-star accretion rate for Rstar = 0.67 R⊙	0, 1.14 × 10−8, 5.17 × 10−8	M⊙ yr−1
${\dot{M}}_{\mathrm{disk},2}$	Disk-to-star accretion rate for Rstar = 2.09 R⊙	3.67 × 10−7, 1.63 × 10−6	M⊙ yr−1
${\dot{M}}_{\mathrm{disk},3}$	Disk-to-star accretion rate for Rstar = 6.61 R⊙	1.14 × 10−5, 5.15 × 10−5, 1.66 × 10−4	M⊙ yr−1
Rtrunc	Magnetospheric truncation radiusa	3	R*
fspot	Fractional area of the hot spots on the starb	0.01	⋯
Envelope Properties
Renv	Envelope outer radiusc	10,000	au
ρ1000	Envelope density at 1000 aud	0.0, 1.19 × 10−20, 1.78 × 10−20, 2.38 × 10−20,	g cm−3
		5.95 × 10−20, 1.19 × 10−19, 1.78 × 10−19,
		2.38 × 10−19, 5.95 × 10−19, 1.19 × 10−18,
		1.78 × 10−18, 2.38 × 10−18, 5.95 × 10−18,
		1.19 × 10−17, 1.78 × 10−17, 2.38 × 10−17,
		5.95 × 10−17, 1.19 × 10−16, 1.78 × 10−16
Rc	Centrifugal radius of TSC envelope	= Rdisk	au
θ	Cavity opening angle	5, 15, 25, 35, 45	degrees
bcav	Exponent for cavity shapee (polynomial)	1.5	⋯
zcav	Vertical offset of cavity wall	0	au
Derived Parameters
L*	Stellar luminosityf	0.1, 1, 10	L⊙
Ltot	Total luminosity (star + accretion)g	0.1, 0.3, 1.0, 3.1, 10.1, 30.2, 101, 303	L⊙
Parameters for Model SEDs
i	Inclination angle	18.2, 31.8, 41.4, 49.5, 56.7,	degrees
		63.3, 69.5, 75.6, 81.4, 87.2
	Aperture radii for model fluxesh	420, 840, 1260, 1680, ..., 10080	au

Notes. The dust opacities used for these models are those called "icsgra3" from Ormel et al. (2011).

^aThis radius applies to the gas. The inner disk radius for the dust is equal to the dust destruction radius. The scale height of the disk at the dust sublimation radius is set to the hydrostatic equilibrium solution. ^bThe hot spots are caused by the accretion columns that reach from the magnetospheric truncation radius to the star. ^cThe inner envelope radius is set to the dust destruction radius. ^dThe actual input parameter for the Whitney code is the envelope infall rate, which can be derived from ρ₁₀₀₀ using Equation (2). The first six ρ₁₀₀₀ values correspond to envelope infall rates of 0, 5.0 × 10⁻⁸, 7.5 × 10⁻⁸, 1.0 × 10⁻⁷, 2.5 × 10⁻⁷, and 5.0 × 10⁻⁷ M_⊙ yr⁻¹; the other values can be similarly deduced. ^eThe cavity walls are assumed to have a polynomial shape; no material is assumed to lie inside the cavity. Also, the ambient density (outside the envelope) is 0. ^fThe three values of L_* correspond to the three different stellar radii. ^gThe total luminosities combine the stellar luminosities and the accretion luminosities (which depend on ${\dot{M}}_{\mathrm{disk}}$). ^hFor each model, the emitted fluxes are calculated for 24 apertures ranging from 420 to 10,080 au, in steps of 420 au.

Download table as:  ASCIITypeset image

The total luminosity for each system consists of the stellar luminosity (derived from a 4000 K stellar atmosphere model), the accretion luminosity resulting from material accreting through the disk down to the disk truncation radius, and the accretion luminosity from the hot spots on the stellar surface, where the accretion columns, which start at the magnetospheric truncation radius, land (these columns are not included in the modeled density distribution, since they do not contain dust and do not have a source of opacity in the radiative transfer models). Typically, the accretion luminosity from the hot spots is much larger than the disk accretion luminosity; in our models, the former is about a factor of 9 larger than the latter. We chose three different stellar radii, 0.67, 2.09, and 6.61 R_⊙ (with the same stellar temperature), resulting in three different stellar luminosities. Since both components of the accretion luminosity depend on the disk accretion rate, choosing a total of eight different disk accretion rates (three for the 0.67 R_⊙ star, two for the 2.09 R_⊙ star, and three for the 6.61 R_⊙ star) results in eight values for the total luminosity used in the grid (see Table 3). The input spectrum produced by the central protostar depends on the relative contributions from the intrinsic stellar luminosity (which peaks at 0.7 μm) and the accretion luminosity (which is radiated primarily in the UV). In the models, it can be altered to some degree by choosing different combinations of the disk accretion rate and stellar radius (the former affects only the accretion luminosity, while the latter affects both the stellar and accretion luminosity). However, the effect of the input spectrum on the output SED is negligible. Consequently, we cannot reliably measure the relative contributions of stellar and accretion luminosity through our SED fits. Instead, we adjusted the particular values for the stellar radius and disk accretion rate to set the values of the total luminosity.

For our model grid, we chose four values for the disk outer radius, which we set equal to the centrifugal radius (R_c). In a TSC model, the centrifugal radius is the position in the disk where material falling in from the envelope accumulates; owing to envelope rotation, material from the envelope's equatorial plane lands at R_c, while material from higher latitudes falls closer to the star. The disk could extend beyond R_c, but in our models it ends at R_c. In this work, we use the terms "disk (outer) radius" and "centrifugal radius" interchangeably. The primary effect of R_c is to set the rotation rate of the infalling gas and thereby determine the density structure of the envelope (Kenyon et al. 1993).

The largest number of parameter values in our grid is for the envelope infall rate. The envelope infall rate used as an input in the Whitney code sets the density of the envelope for a given mass of the protostar. Since the SED depends on the density of the envelope (and not directly on the infall rate, which is only inferred from the density and the acceleration due to gravity from the central protostar), in this work we report a reference envelope density instead of the envelope infall rate as one of our model parameters. For the TSC model, the envelope infall rate ${\dot{M}}_{\mathrm{env}}$ and the reference density at 1 au in the limit of no rotation (R_c = 0) are related as follows (see Kenyon et al. 1993):

$$\begin{eqnarray}&&{\rho }_{1}=5.318\times {10}^{-14}\left(\displaystyle \frac{{\dot{M}}_{\mathrm{env}}}{{10}^{-5}{M}_{\odot }\;{\mathrm{yr}}^{-1}}\right){\left(\displaystyle \frac{{M}_{*}}{1{M}_{\odot }}\right)}^{-1/2}{\rm{g}}\;{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 1 }$$

where M_* is the mass of the central protostar, which is assumed to be 0.5 M_⊙ in our model grid. The density distribution in the envelope follows a power law, ρ ∝ r^−3/2, at radii larger than the centrifugal radius, R_c, but then flattens as a result of the rotation of the envelope. The density reported by ρ₁ assumes a spherically symmetric envelope with a −3/2 power-law exponent valid down to the smallest radii, and it is higher than the angle-averaged density of a rotating envelope at 1 au. To quote densities that are closer to actual values found in the modeled rotating envelopes (which have R_c values ranging from 5 to 500 au), we report ρ₁₀₀₀, the density at 1000 au for a ρ ∝ r^−3/2 envelope with a 0.5 M_⊙ protostar:

$$\begin{eqnarray}{\rho }_{1000} & = & {\rho }_{1}{\left(\displaystyle \frac{1}{1000}\right)}^{3/2}\\ & = & 2.378\times {10}^{-18}\left(\displaystyle \frac{{\dot{M}}_{\mathrm{env}}}{{10}^{-5}{M}_{\odot }\;{\mathrm{yr}}^{-1}}\right)\;{\rm{g}}\;{\mathrm{cm}}^{-3}.\end{eqnarray} \tag{ 2 }$$

Thus, the range of reference densities probed in our model grid, from 1.2 × 10⁻²⁰ to 1.8 × 10⁻¹⁶ g cm⁻³ (see Table 3), would correspond to envelope infall rates from 5.0 × 10⁻⁸ to 7.5 × 10⁻⁴ M_⊙ yr⁻¹, assuming M_* = 0.5 M_⊙ (this does not account for a reduction of the infalling mass due to clearing by outflow cavities). In Figure 5, we show the radial density profiles for two TSC models with 5 and 500 au centrifugal radii. The density profiles are azimuthally symmetric and show the flattening of the density distribution inside R_c due to envelope rotation. These plots demonstrate that the density ρ₁ is much higher than the angle-averaged density at 1 au; ρ₁₀₀₀ seems to yield more physical values for the density in the envelope at 1000 au, even for R_c values of 500 au.

Zoom In Zoom Out Reset image size

As can be seen from the values of the envelope density in Table 3, there is one set of models with an envelope density of 0. These are models that do not contain an envelope component; the entire excess emission is caused by the circumstellar disk. If an object is best fit by such a model, it would indicate that it is more evolved, having already dispersed its envelope.

The cavities in our models range from 5° to 45° and are defined such that $z\propto {\tilde{r}}^{1.5}$, where $\tilde{r}$ and z are the cylindrical coordinates for the radial and vertical direction, respectively, and ${\tilde{r}}_{\mathrm{max}}={z}_{\mathrm{max}}\;\mathrm{tan}\theta$, with θ defined as the cavity opening angle that is specified in the parameter file of the Whitney radiative transfer code. In this code, z_max is set to the envelope outer radius. Thus, a polynomial-shaped cavity, which is wider at smaller $\tilde{r}$ values and then converges toward the specified opening angle, is somewhat larger than this opening angle at the outer envelope radius (see Figure 6). This effect is most noticeable at larger cavity opening angles, but negligible for small cavities. A different definition of the cavity size, where ${\tilde{r}}_{\mathrm{max}}={R}_{\mathrm{env}}\mathrm{sin}\theta$ and ${z}_{\mathrm{max}}={R}_{\mathrm{env}}\mathrm{cos}\theta$ (with R_env as the envelope outer radius), results in z values that are a factor of $1/\mathrm{cos}\theta$ larger, and thus the cavity reaches the specified opening angle at the outer envelope radius. For this work, the adopted definition of the cavity opening angle is inconsequential, but it becomes relevant when comparing the results of SED modeling to scattered light images that reveal the actual cavity shape and size. We also note that in our models the cavities are evacuated of material, so there is no dust and gas inside the cavity; in reality, there might be some low-density material left that would add to the scattered light (see Fischer et al. 2014).

Zoom In Zoom Out Reset image size

Figures 7–11 display a few examples of model SEDs from our grid to show the effect of changing those model parameters that influence the resulting SED the most. The inclination angle has a strong effect on the near- and mid-infrared SED (Figure 7). While a low inclination angle results in an overall flat SED in this wavelength region, increasing the inclination angle causes a deeper silicate absorption feature at 10 μm and a steep slope beyond it. The far-infrared to millimeter SED is not affected by the inclination angle, since emission at these wavelengths does not suffer from extinction through the envelope.

Zoom In Zoom Out Reset image size

The cavity opening angle affects the SED shape at all wavelengths (Figure 8). A small cavity only minimally alters the SED compared to a case without a cavity; there is still a deep silicate absorption at 10 μm and steep SED slope, but the cavity allows some scattered light to escape in the near-IR. A larger cavity results in higher emission at near- and mid-infrared wavelengths and reduced emission in the far-infrared. The effect of the cavity on the SED would change if a different shape for the cavity walls were adopted. For example, cavities where the outer wall follows the streamlines of the infalling gas and dust evacuate less inner envelope material than our polynomial-shaped cavities, resulting in deeper silicate absorption features and steeper mid-infrared SED slopes for the same cavity opening angle (see Furlan et al. 2014). Thus, our cavity opening angles are tied to our assumed cavity shape.

Zoom In Zoom Out Reset image size

The effect of the centrifugal radius is somewhat similar to those of the cavity opening angle and inclination angle, but less pronounced (Figure 9). Small disk radii imply more slowly rotating, less flattened envelopes and depress the near- and mid-infrared fluxes more than larger disk radii, but even with large disk radii (and more flattened envelopes) there is still sufficient envelope material along the line of sight to cause a pronounced 10 μm absorption feature. Overall, our models do not directly constrain the size of the disk; the opacity is dominated by the envelope. Furthermore, the flattening of the envelope that is determined by R_c has a similar effect on the SED as changing the outflow cavity opening angle.

Zoom In Zoom Out Reset image size

Changing the envelope density causes shifts in the SED in terms of both wavelength and flux level: the higher the envelope density, the less flux is emitted at shorter wavelengths, and the more the peak of the SED shifts to longer wavelengths (Figure 10). Deeply embedded protostars have SEDs that peak at λ > 100 μm, steep mid-IR SED slopes, and deep silicate absorption features. The effect of the envelope density on the SED is different from that of the inclination angle, especially in the far-IR: while the SED is not very sensitive to the inclination angle in this wavelength region, the ratio of, e.g., 70 and 160 μm fluxes changes considerably depending on the envelope density.

Zoom In Zoom Out Reset image size

The total luminosity of the source has an effect on the overall emission level of the protostar, but does not strongly affect the SED shape. The main effect is that the peak of the SED shifts to longer wavelengths as the luminosity decreases (λ_peak ∝ L^−1/12; Kenyon et al. 1993). Especially when comparing models with L_tot values that differ by a factor of a few, the SED shapes are similar (Figure 11). Thus, one could scale a particular model by a factor between ∼0.5 and 2 and get a good representation of a protostar that is somewhat fainter or brighter, without having to rerun the model calculation with the different input luminosity.

Zoom In Zoom Out Reset image size

The model fluxes are computed for 24 different apertures, ranging from 420 to 10,080 au in steps of 420 au (which corresponds to 1'' at the assumed distance of 420 pc to the Orion star-forming complex). For these SED fluxes, no convolution with a PSF is done, and therefore the spatial distribution of the flux is solely due to the extended nature of protostars. Since the envelope outer radius is chosen to be 10,000 au, the largest aperture encompasses the entire flux emitted by each protostellar system. However, most of the near- and mid-infrared emission comes from smaller spatial scales, so an aperture of about 5000 au will already capture most of the flux emitted at these wavelengths.

For a more accurate comparison of observed and model fluxes, in each infrared photometric band where we have data available, we interpolate model fluxes from the two apertures that bracket the aperture used in measuring the observed fluxes (4'' for 2MASS, 24 for IRAC, PSF photometry for MIPS 24 μm, with a typical FWHM of 6'', 96 for PACS 70 and 100 μm, 128 for PACS 160 μm). For the IRS data points, we use fluxes interpolated for a 53 aperture, since the spectra are composed of two segments, SL (5.2–14 μm; slit width of 36) and LL (14–38 μm, slit width of 105), and, if any flux mismatches were present, the SL segment was typically scaled to match the LL flux level at 14 μm (see, e.g., Furlan et al. 2008). So, fluxes measured in an aperture with a radius of 53 roughly correspond to fluxes from a 106-wide slit.

Given that our targets are typically extended and that the near- to mid-infrared data have relatively high spatial resolution, measuring fluxes in small apertures (a few arcseconds in radius) will truncate some of the object's flux, so it is important to choose similar apertures for the model fluxes. From about 30 to 100 μm, the model fluxes calculated for smaller apertures are not very different from the total flux (i.e., the flux from the largest aperture), which is a result of the emission profile in the envelope and the lower spatial resolution at longer wavelengths. To check whether extended source emission in the far-infrared might affect the flux we measure in our models, we calculated a small set of model images at 160 μm, convolved them with the PACS 160 μm PSF, and compared the fluxes from the model images to those written out for the model SEDs (which we refer to as "SED fluxes"; these are the fluxes from the models in the grid). Model images would be the most observationally consistent way to measure the flux densities, but they are too computationally expensive and would not represent a significant gain.

In Figure 12 we show the fluxes derived for a particular model at 160 μm using different methods. The fluxes measured in the convolved model image are lower than the SED fluxes; this is caused by the wide PACS 160 μm PSF, which spreads flux to very large radii. Since the shape of the PSF is known, we can correct for these PSF losses (assuming a point source and using standard aperture corrections). The fluxes corrected for these PSF losses are very similar to the SED fluxes, typically within ∼5%–10% at apertures larger than 5''. Since our observed fluxes correspond to these PSF-corrected fluxes (we apply aperture corrections to our fluxes measured in a 128 aperture to account for PSF losses), adopting the SED fluxes from the largest aperture would yield model fluxes that are somewhat too high. Thus, we chose to adopt the SED flux measured in a 128 aperture as a good approximation for the model flux we would get if we had model images available for all models in the grid and measured aperture-corrected fluxes in these images. We note that in our PACS data, the 160 μm sky annulus, which extends from 128 to 256 (see B. Ali et al. 2016, in preparation), can include extended emission from surrounding material and also some envelope emission. In these cases, we often used PSF photometry to minimize contamination from nearby sources and nebulosity; however, PSF fitting was not used for more isolated sources since the envelopes can be marginally resolved at 160 μm and thus deviate slightly from the adopted PSF shape.

Zoom In Zoom Out Reset image size

For the SABOCA and LABOCA data, beam fluxes were adopted; the FWHM of the SABOCA beam is 73, while for the LABOCA beam it is 19''. In order to determine which aperture radius corresponds best to beam fluxes, we created a similar set of model images as above at 350 and 870 μm, convolved them with Gaussian PSFs, and measured fluxes in the model images using different apertures (see Figures 13 and 14, where we show the results for one model). Fluxes measured in the convolved model image are smaller than the SED fluxes, especially at aperture radii smaller than the FWHM of the beam. We find that the beam fluxes for SABOCA and LABOCA are best matched by SED fluxes from apertures with radii half the size of the FWHM of the beam, i.e., 365 for SABOCA and 95 for LABOCA (thus, the aperture sizes are the same as the beam FWHM). This is again an idealized situation, since the measured SABOCA and LABOCA beam fluxes also include extended emission (if the source lies on top of background emission), and thus they could be higher than those from the model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In our models, the luminosity is determined by the central protostar and the accretion; no external heating is included. The interstellar radiation field (ISRF) could increase the temperature in the outer envelope regions, thus causing an increase in the longer-wavelength fluxes (e.g., Evans et al. 2001; Shirley et al. 2002; Young et al. 2003). It is expected that external heating has a noticeable effect only on low-luminosity sources (≲1 L_⊙), while objects with strong internal heating are not affected by the ISRF. Moreover, the strength of the ISRF varies spatially (Mathis et al. 1983), and thus its effect on each individual protostar is uncertain. Nonetheless, in the following we estimate the effect of external heating on model fluxes by using a different set of models.

For this model calculation, we used the 2012 version of the Whitney radiative transfer code (Whitney et al. 2013), which allows for the inclusion of external illumination by using the ISRF value in the solar neighborhood from Mathis et al. (1983); to vary the ISRF strength, the adopted value can be scaled by a multiplicative factor and extinguished by a certain amount of foreground extinction. We calculated a small number of models with and without external heating and then compared their far-infrared and submillimeter fluxes. One set of models has L_tot = 0.1 L_⊙, R_c = 100 au, θ = 15°, and four different reference densities ρ₁₀₀₀, ranging from 2.4 × 10⁻¹⁷ to 2.4 × 10⁻²⁰ g cm⁻³. The other set has the same parameters except for L_tot, which is 1.0 L_⊙. We calculated models without external heating, with heating from an ISRF equal to that in the solar neighborhood, and with ISRF heating 10 times the solar neighborhood value. For these models, we did not include any foreground extinction for the ISRF; thus, the ISRF heating in these models can be considered an upper limit—especially the 10-fold increase over the ISRF in the solar neighborhood represents an extreme value. Figure 15 shows a few examples of model SEDs with and without external heating. External heating results in flux increases in the far-IR and submillimeter; as expected, it affects low-luminosity sources more, and its effects are also more noticeable for higher-density envelopes.

Zoom In Zoom Out Reset image size

For a more quantitative comparison of model fluxes in the far-IR and submillimeter, we computed the fluxes for each model in six different bands, those of MIPS 24 μm, PACS 70, 100, and 160 μm, and SABOCA (350 μm) and LABOCA (870 μm), using apertures as described in Section 4.2. The model fluxes are affected by poorer signal-to-noise ratios at the longest wavelengths, so the 870 μm fluxes are less reliable. We subtracted the fluxes of the models without external heating (F_{no.ext.heating}) from those with external heating (F_ext.heating) to determine the flux excess due to external heating. The ratios of these excess fluxes and the model fluxes with external heating ($({F}_{{\rm{ext.heating}}}-{F}_{{\rm{no.ext.heating}}})/{F}_{{\rm{ext.heating}}}$) are shown in Figure 16. Given that these ratios depend on the inclination angle to the line of sight, we show them as average values for all 10 inclination angles, as well as the range subtended by all inclination angles. We note overall smaller flux ratios at 350 μm owing to the smaller aperture size chosen in this wave band (see Section 4.2).

Zoom In Zoom Out Reset image size

Our analysis shows that heating by the ISRF results in flux increases in the far-IR and submillimeter that are about a factor of 2–3 higher for envelopes of low-luminosity sources (L_tot = 0.1 L_⊙) than for those with higher luminosity. Also, the effect of external heating is more noticeable at longer wavelengths (where apertures/beams are also larger) than at shorter ones; given our chosen apertures, the largest effect occurs at 160 and 870 μm. We also note that the flux increases due to heating by the ISRF are smallest for the lowest ρ₁₀₀₀ value probed, 2.4 × 10⁻²⁰ g cm⁻³; at 160 μm, the flux increase is largest for intermediate envelope densities. Finally, the flux increases in the far-IR and submillimeter are far larger for a solar-neighborhood ISRF scaled by factor of 10 than for an unscaled ISRF; for the L_tot = 0.1 L_⊙ models, an unscaled ISRF increases the fluxes from a few percent (at ≲100 μm) to 50% (at 870 μm), while an ISRF scaled by a factor of 10 increases these fluxes by 30%–75%. Thus, for low-luminosity protostars, up to ∼75% of a protostar's 870 μm flux could be due to external heating, if the environment is dominated by an extremely strong ISRF.

To estimate how the contribution of external heating would modify derived model parameters, in Figures 17 and 18 we compare model SEDs that include external heating by an ISRF 10 times stronger than in the solar neighborhood and model SEDs without this additional heating. For the latter, we used models from our model grid and tried to reproduce the SEDs with external heating. For the models with L_tot = 0.1 L_⊙, the effect of external heating can be reproduced by increasing the luminosity by factors of a few, increasing ρ₁₀₀₀ by up to an order of magnitude, and increasing the cavity opening angle and inclination angle by a small amount. For the L_tot = 1.0 L_⊙ models, just increasing the reference density by a factor of 2.5 results in a good match to the long-wavelength emission of our externally heated models; however, the shorter-wavelength flux is either under- or overestimated. A better match is achieved with models having the same reference density as the externally heated models, but with slightly larger cavity opening angles and inclination angles, and luminosities about a factor of 2 larger. Thus, if the far-IR and submillimeter fluxes were contaminated by emission resulting from extremely strong external heating, a model fit using models from our grid (which does not include external heating) could overestimate the envelope density by up to an order of magnitude and the luminosity by a factor of 2–5. The cavity opening and inclination angles would also be more uncertain, but not by much. For a more realistic scenario with more modest external heating (which would also include the effect of local extinction), the effect on model parameters would be smaller.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the latter point, we explored the effect of extinction on the ISRF by calculating a few more models with L_tot = 0.1 L_⊙, R_c = 100 au, θ = 15°, ρ₁₀₀₀ = 2.4 × 10⁻¹⁸ g cm⁻³, an ISRF 10 times stronger than that in the solar neighborhood, and A_V values for the ISRF of 2.5, 10, 20, and 50. The model SEDs are shown in Figure 19. Compared to ISRF heating without any foreground extinction, already A_V = 2.5 causes a decrease by a factor of 1.5–2 in the overall emission at far-IR wavelengths. With A_V of 10 and 20, the far-IR emission decreases by factors of up to ∼3.5 and 4, respectively, compared to a strong ISRF that is not extinguished. The fraction of excess emission due to external heating at 160 μm decreases from an average of 0.8 for A_V = 0 (see Figure 16) to 0.6, 0.3, and 0.2 for A_V = 2.5, 10, and 20, respectively. Therefore, considering that typical A_V values in Orion are ∼10–20 mag (Stutz & Kainulainen 2015), it is likely that the effect of external heating on model parameters of low-luminosity sources does not exceed a factor of ∼2 in luminosity and ∼5 in envelope density.

Zoom In Zoom Out Reset image size

Figure 21 displays the histogram of R values of the best model fits for the 330 objects in our HOPS sample that have Spitzer and Herschel data (more than two data points at different wavelengths) and are not contaminants (see Section 2). The median R value is 3.10, while the mean value is 3.29. Fitting a Gaussian to the histogram at R ≤ 7 yields 3.00 and 2.24 as the center and FWHM of the Gaussian, respectively. The distribution of R values implies that, on average, the model deviates by about three times the average fractional uncertainty from the data. This is not unexpected, given that we fit models from a grid to observed SEDs that span almost three orders of magnitude in wavelength range, with up to 29 data points. The fewer the data points, the easier it is to achieve a good fit; in fact, the eight protostars with R < 1, HOPS 371, 391, 398, 401, 402, 404, 406, and 409, have SEDs with measured flux values at only 4–5 points. Starting at R values of about 1, R can be used as an indicator of the goodness of fit. However, in some cases a noisy IRS spectrum can increase the R value of a fit that, judged by the photometry alone, does not deviate much from the observed data points. In other cases, mismatches between different data sets, like offsets between the IRAC fluxes and the IRS spectrum, can result in larger R values. These might be interesting protostars affected by variability and are thus ideal candidates for follow-up observations.

Zoom In Zoom Out Reset image size

When looking at the SED fits in Figure 1 (and the corresponding R values in Table 1), we estimate that an R value of up to ∼4 can identify a reliable fit (with some possible discrepancies between data and model in certain wavelength regions). When R gets larger than about 5, the discrepancy between the fit and the observed data points usually becomes noticeable; the fit might still reproduce the overall SED shape but deviate substantially from most measured flux values.

In Figure 22, we show the histogram of R values separately for the three main protostellar classes in our sample. The median R value decreases from 3.27 for the Class 0 protostars to 3.18 for the Class I protostars to 2.58 for the flat-spectrum sources. There are 4 Class 0 protostars and 4 Class I protostars with R values between 1.0 and 2.0, but 17 flat-spectrum sources in this R range. These numbers translate to 17% of the flat-spectrum sources in our sample, 4% of the Class 0 protostars, and 3% of the Class I protostars. When examining objects' R values between 2.0 and 4.0, there are 51 Class 0 protostars (55% of Class 0 protostars in the sample), 91 Class I protostars (73% of the Class I sample), and 74 flat-spectrum sources (73% of the flat-spectrum sample).

Zoom In Zoom Out Reset image size

Thus, close to 90% of flat-spectrum sources are fit reasonably well (R values <4), representing the largest fraction among the different classes of objects in our sample. This could be a result of their source properties being well represented in our model grid, but also lack of substantial wavelength-dependent variability (see, e.g., Günther et al. 2014), which, if present, would make their SEDs more difficult to fit. About three-quarters of Class I protostars also have best-fit models with R < 4; this fraction drops to about two-thirds for the Class 0 protostars. The latter group of objects often suffers from more uncertain SEDs owing to weak emission at shorter wavelengths (which, e.g., results in a noisy IRS spectrum); they might also be more embedded in extended emission, such as filaments, which can contaminate the far-IR to submillimeter fluxes. Another factor that could contribute to poor fits is their presumably high envelope density, which places them closer to the limit in parameter space probed by the model grid. Overall, 75% of the best-fit models of the protostars in our sample have R < 4.

When examining the SED fits of objects with R values larger than 5.0, several have very noisy IRS spectra (HOPS 19, 38, 40, 95, 164, 278, 316, 322, 335, 359). In a few cases the measured PACS 100 and 160 μm fluxes seem too high compared to the best-fit model (e.g., HOPS 189), which could be an indication of contamination by extended emission surrounding the protostar.

Of particular interest are objects where variability likely plays a role in a poor fit. As mentioned in Section 3, variability among protostars is common; we found in Appendix A that about 5% of our targets display noticeable (≳50%) mismatches between the IRS, IRAC, and MIPS fluxes that could be due to intrinsic variability. The SED fits of objects for which the flux mismatches between IRS and IRAC and between IRS and MIPS are different are particularly affected, since in that case we did not scale the IRS spectrum to match the MIPS 24 μm flux. HOPS 228 exemplifies such a case: there is a clear discrepancy between the IRAC and IRS fluxes (a factor of 2.1–2.7) and also between MIPS 24 μm and IRS (a factor of 0.8); even though the fit gives more weight to the IRS data, they are not fit well, especially the silicate absorption feature. The R value of 5.74 for the fit of HOPS 228 reflects the discrepant data sets and poor fit. HOPS 223 is another case where the IRS fluxes do not match the shorter-wavelength data (they are more than an order of magnitude larger); however, it is a known FU Ori source (see Fischer et al. 2012), and the SED presented here contains both pre- and post-outburst data. The model fit is very poor, which can also be gauged by the R value of 8.41.

There are also objects with overall good fits whose SEDs show discrepancies that may be signs of variability or contamination. For example, for the Class I protostar HOPS 71 the IRAC fluxes are a factor of 1.8–2.4 lower than the IRS fluxes in the 5–8 μm region, and also the MIPS flux is about 20% lower. The best-fit model (R = 3.63) fits the SED extremely well beyond about 6 μm, with some discrepancy at shorter wavelengths. There is a source just 11'' from HOPS 71 that is detected in 2MASS and Spitzer data, but not by PACS; this object, HOPS 72, is likely an extragalactic object (see Appendix C.2.2) that could contaminate the IRS fluxes. Thus, in this case, wavelength-dependent contamination by a companion could explain the discrepancies observed in the SED.

Another example is HOPS 124, which is a deeply embedded Class 0 protostar. For this object, the mismatch between IRS and IRAC and MIPS fluxes decreases with increasing wavelength (from a factor of 2.5 to a factor of 1.4); for the SED fit, the IRS spectrum was scaled by 0.7 to match the MIPS 24 μm flux. As with HOPS 71, there is a nearby source that could contaminate some of the fluxes, especially at shorter wavelengths: HOPS 125, a flat-spectrum source, lies 98 from HOPS 124 and is brighter than HOPS 124 out to ∼20 μm, but then much fainter at longer wavelengths. The best-fit model of HOPS 124 (R = 2.43) matches the mid- to far-IR photometry and also most of the IRS spectrum well.

As an example of a probably variable flat-spectrum source, HOPS 132 has IRAC fluxes that lie a factor of 1.3–1.7 above those of IRS and a MIPS 24 μm flux that is a factor of 0.6 lower. It does not have a close companion; the nearest HOPS source, HOPS 133, is 27'' away. The IRS spectrum was not scaled, and since the SED fitter gave more weight to the spectrum, it is fit well, but the IRAC photometry is underestimated and the MIPS photometry overestimated. Nonetheless, the R value of the best fit is 2.87.

Overall, the SED fits of objects that are likely variable or suffer from some contamination are less reliable, but it is not always clear from the R value of the best fit. The SED fitting procedures assume that the protostars are not variable, so when large mismatches between different data sets are present, the fit will appear discrepant with at least some of the observed data points, but the R value would not end up particularly high if, e.g., the IRS spectrum was fit exceptionally well. However, given the data sets we have for these protostars, our SED fits will still yield the best possible estimate for the protostellar parameters describing these systems.

The histogram of best-fit ρ₁₀₀₀ values (which is the density of the envelope at 1000 au; see Section 4.1) is shown in Figure 23. The median value of the distribution amounts to 5.9 × 10⁻¹⁹ g cm⁻³; this corresponds to a ρ₁ value of 1.9 × 10⁻¹⁴ g cm⁻³. There is a spread in values: 69 objects have densities ρ₁₀₀₀ smaller than 5.0 × 10⁻²⁰ g cm⁻³ (6 of them have actually no envelope), 89 fall in the 5.0 × 10⁻²⁰ to 5.0 × 10⁻¹⁹ g cm⁻³ range, 96 are between 5.0 × 10⁻¹⁹ and 5.0 × 10⁻¹⁸ g cm⁻³, 60 between 5.0 × 10⁻¹⁸ and 5.0 × 10⁻¹⁷ g cm⁻³, and 16 have ρ₁₀₀₀ values larger than 5.0 × 10⁻¹⁷ g cm⁻³.

Zoom In Zoom Out Reset image size

We also calculated the envelope mass (M_env) within 2500 au for the best-fit models (see Figure 24 for their distribution). The 2500 au radius is close to half the FWHM of the PACS 160 μm beam at the distance of Orion (i.e., ∼6'') and thus roughly represents the spatial extent over which we measure the SEDs. This envelope mass is determined from the integrated envelope density of our best-fit models, with allowances made for outflow cavities, and thus only valid in the context of our models. The median envelope mass within 2500 au amounts to 0.029 M_⊙. The majority of protostars have model-derived masses in the inner 2500 au of their envelopes around 0.1 M_⊙; just 22 objects have M_env (<2500 au) larger than 1.0 M_⊙. Of the 330 modeled objects, 291 have M_env (<2500 au) smaller than 0.5 M_⊙ (6 of these 291 objects have no envelope).

Zoom In Zoom Out Reset image size

Figure 25 contains the histogram of the total luminosities derived from the best-fit models. These luminosities consist of the stellar, disk accretion, and accretion shock components. The median total luminosity amounts to 3.02 L_⊙, while the values cover four orders of magnitude, from 0.06 L_⊙ (for HOPS 336) to 607 L_⊙ (for HOPS 288 and 361). Since the minimum and maximum values for the total luminosity in our grid amount to 0.1 and 303.5 L_⊙, respectively, and our scaling factors range from 0.5 to 2.0, our fitting procedure can return best-fit luminosities that range from 0.05 to 607 L_⊙. Thus, two protostars are reaching the upper limit allowed for total luminosities in our grid; it is possible that even better fits could be achieved by increasing the luminosity further.

Zoom In Zoom Out Reset image size

From the distribution of best-fit outer disk radii in Figure 26, it is apparent that most protostars are fit by small disks whose radius is only 5 au. Since the outer disk outer radius is the centrifugal radius in our models, infalling material from the envelope tends to accumulate close to the star for most sources. Thus, the disk radius is tied to the envelope structure; a small centrifugal radius implies higher envelope densities at smaller radii and a less flattened envelope structure. The median disk radius is 50 au, but the number of objects with disk radii ≥50 au is roughly evenly split among the values of 50, 100, and 500 au.

Zoom In Zoom Out Reset image size

The distribution of best-fit cavity opening angles is displayed in Figure 27. Most protostars seem to have either very small (5°) or very large (45°) cavities; the median value is 25°. When dividing the envelope densities by cavity opening angle (see Figure 28, left column), differences emerge: the distributions of ρ₁₀₀₀ values are significantly different when comparing objects with θ = 5° and θ ≥ 35°, objects with θ = 15° and θ ≥ 25°, and objects with θ = 25° and θ = 45°. The Kolmogorov–Smirnov (K-S) tests yield significance levels that these subsamples are drawn from the same parent population of ≲0.015. Thus, there seems to be a difference in the distribution of envelope densities among the best-fit models with smaller cavity opening angles and those with larger cavities. Protostars with larger cavities (≥35°) tend to have higher envelope densities (their median ρ₁₀₀₀ values are about an order of magnitude larger compared to objects with cavities ≤15°).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 28 (middle column) also shows the distribution of total luminosities for the different cavity opening angles. The only significant difference can be found for the θ = 5° histogram as compared to the histograms for larger θ values (K-S test significance level ≲0.03); the luminosities of models with θ = 5° have a different distribution, and also their median value is 1.45 L_⊙, as compared to ∼3–5 L_⊙ for the models with larger cavities. So, protostars with small cavities seem to have lower total luminosities.

The distribution of centrifugal radii for different cavity opening angles (right column in Figure 28) shows that, independent of cavity size, most objects have R_disk = 5 au. However, the distribution among the four different disk radii becomes flatter for the largest cavity opening angles; the histograms for θ = 35° and θ = 45° are very similar (K-S test significance level of 0.98). There is also no significant difference (K-S test values >0.075) between the θ = 15° and θ = 25° histograms and between the θ = 5° and θ ≥ 35° histograms. The distributions of disk radii for the other cavity opening angles are all different from one another (K-S test significance levels <0.015). Overall, Figure 28 shows that protostars best fit by models with large cavity opening angles are also fit by models with higher envelope densities and larger centrifugal radii.

In Figure 29, we show the distribution of the inclination angles for the best-fit models. There is a clear concentration of models in the 60°–70° range; the median inclination angle is 63°. This median value is close to 60°, which is where the probability for isotropically distributed inclination angles reaches 50% (i.e., the probability of observing an inclination angle less than 60° is the same as the probability of observing i > 60°). However, the details of the distributions differ. The cumulative probability of finding an inclination angle less than a certain value, i_c, is $1-\mathrm{cos}({i}_{c})$, assuming a random distribution of inclination angles. For inclination angles i₁ and i₂, the probability for i₁ < i < i₂ is cos(i₁) − cos(i₂). Thus, since the inclination angles in our model grid were chosen to be equally spaced in cos(i) (there are five values <60° and five values >60°), one would expect a flat distribution in Figure 29 if the best-fit inclination angles were randomly distributed (see the green dashed histogram). However, we find a distribution peaked at 63° and 70°. This can also be seen in Figure 30, where we compare our observed cumulative distribution of inclination angles to that of randomly distributed ones. Our distribution shows a deficit at inclination angles below 60° and is just slightly higher at large inclination angles. A K-S test of the two distributions yields a 5.6% chance that they are drawn from the same parent distribution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To examine whether the distribution of envelope parameters changes with inclination angle (which could imply a degeneracy), Figure 31 shows the reference envelope density ρ₁₀₀₀, the total luminosity, and the cavity opening angle binned by three ranges of inclination angles. None of the three model parameters show a significantly different distribution for any of the inclination bins (K-S test significance levels are ≳0.1, except for the cavity opening angles for the lowest and middle inclination range, for which the K-S test significance value is 0.02). The median ρ₁₀₀₀ values for the i = 18°–41°, 49°–63°, and 69°–87° inclination bins are all 5.9 × 10⁻¹⁹ g cm⁻³. Even though not shown in Figure 31, the objects whose best-fit model does not include an envelope are only found at i ≥ 49°. It is noteworthy that protostars with the highest envelope densities do not have inclination angles in the 69°–87° range; it is not clear whether this is an observational bias, whether our observed sample does not contain high-density, edge-on protostars, or whether this is due to biases in the fitting procedure and/or model grid. The median values for the total luminosity do not differ by much for the different bins of inclination angle, increasing from 2.9 to 4.1 L_⊙ from the lowest to the middle inclination range and then decreasing to 2.0 L_⊙ for the highest inclination angles. The few protostars with very high L_tot values have large inclination angles (i ≥ 49°). Finally, the distribution of cavity opening angles is quite similar for different ranges in inclination, except for a somewhat larger number of θ = 45° values at intermediate inclination angles. Half the objects in the i = 18°–41° and 69°–87° inclination bins have θ ≤ 15° (with the most common value 5°), while almost half the objects at intermediate inclination angles have θ ≥ 35° (the most common value is 45°).

Zoom In Zoom Out Reset image size

In Figure 32, we show ratios of the total and bolometric luminosities as a function of inclination angle and foreground extinction (i and A_V are adopted from the best model fits). The total luminosity is the intrinsic luminosity from the best-fit model of each object, while the bolometric luminosity is derived by integrating the fluxes of the observed SED. It is expected that L_tot is higher than L_bol for objects seen at higher inclination angles, since for these objects a large fraction of the emitted flux is not directed toward the observer (and thus deriving bolometric luminosities from observed fluxes will underestimate the intrinsic source luminosity). Conversely, objects seen more face-on should have lower L_tot values compared to L_bol. Our data and model fits yield L_tot values that are usually higher than the L_bol values measured from the SED; the discrepancy is larger for the more highly inclined sources. The median L_tot/L_bol ratio is 1.5 for protostars with inclination angles in the 18°–41° range, 2.5 for the i = 49°–63° range, and 3.5 for inclination angles ≥69°. The fact that L_tot > L_bol even for i = 18°–41° could be related to the typically smaller cavity opening angles for this range of inclination angles (see Figure 31); less flux, especially at shorter wavelengths, is detected since the opacity along the line of sight is still high owing to the small cavities.

Zoom In Zoom Out Reset image size

Foreground extinction also plays a role in increasing the L_tot/L_bol ratio. The median ratio of these luminosities increases from 1.8 for the A_V = 0–5 mag range to 5.0 for A_V = 25–30; it decreases somewhat for the next A_V bin, but reaches 5.9 at A_V = 40–50 (the 23 objects with A_V > 50, not shown in Figure 32, have a median L_tot/L_bol ratio of 8.2). Among the 22 objects with best-fit A_V values of 0–5 mag and inclination angles ≤50°, only four have L_tot/L_bol ratios that are larger than 1.5 (they are HOPS 57, 147, 199, and 201; in most cases the model overestimates the near-IR emission).

Figures 33–38 divide the histograms of the best-fit reference density ρ₁₀₀₀, inclination angle, cavity opening angle, total luminosity, disk radius, and foreground extinction, respectively, by protostar class. As explained in Section 3, we divided our targets into Class 0, Class I, flat-spectrum, and Class II objects based on their mid-infrared (4.5–24 μm) spectral index and bolometric temperature (see also Table 1). Thus, Class 0 and I protostars have a spectral index >0.3, and Class 0 protostars have T_bol values <70 K, but, as mentioned in Section 3, there are a few protostars whose spectral index or T_bol value places them very close to the transition region between Class 0 and I or between Class I and flat spectrum. Given that our sample contains just 11 Class II pre-main-sequence stars, we did not include them in the following histograms; they will be discussed in Section 7.2.3.

The distributions of reference densities (Figure 33) are different for all SED classes; none are consistent with being drawn from the same parent population (K–S test significance level <0.01). Overall, Class 0 protostars have higher envelope densities than Class I and flat-spectrum sources; the median ρ₁₀₀₀ values decrease from 5.9 × 10⁻¹⁸ g cm⁻³ to 2.4 × 10⁻¹⁹ g cm⁻³ to 1.2 × 10⁻¹⁹ g cm⁻³ for these three groups. The lower and upper quartiles for ρ₁₀₀₀ are 1.8 × 10⁻¹⁸ and 1.8 × 10⁻¹⁷ g cm⁻³ for the Class 0 protostars and 2.4 × 10⁻²⁰ and 1.2 × 10⁻¹⁸ g cm⁻³ for the Class I and flat-spectrum objects. We will discuss some implications of these differences in derived envelope densities in Section 7.2.

Zoom In Zoom Out Reset image size

For the inclination angles (Figure 34), the distributions are significantly different for all protostellar classes, too (K-S test significance level ≪0.01). As was shown in Figure 29, a random distribution of inclination angles would result in equal numbers of protostars at each value; there is a deficit of Class 0 and Class I protostars at i ≲ 60°, and there are also few Class I protostars and hardly any flat-spectrum sources at the highest inclination angles. The median inclination angle is highest for Class 0 protostars (70°) and then decreases somewhat for Class I protostars (63°) and even more for flat-spectrum sources (57°). Similar to the envelope density, the median inclination angle decreases as one progresses from Class 0 to flat-spectrum sources.

Zoom In Zoom Out Reset image size

In the distributions of cavity opening angles (Figure 35), significant differences can be found between Class 0 and Class I protostars and between Class I protostars and flat-spectrum sources (K-S test significance level ≪0.01). The median cavity opening angle is 15° for the Class I protostars, but 25° for the other two classes. About 40% of Class I protostars have θ = 5°, while the distribution among the different cavity opening angles is flatter for the other two object classes. The large fraction of Class I protostars with small cavities could be the result of degeneracy in model parameters (see Section 7.2) or our assumptions on envelope geometry (see Section 7.4). There are notably few flat-spectrum sources with a 5° cavity opening angle; most of them have cavity opening angles of 15° or 45°.

Zoom In Zoom Out Reset image size

When comparing the total luminosities for the different SED classes (Figure 36), the distribution of L_tot values is different for the Class 0 protostars when compared to the other two classes (K-S test significance level <0.015), but similar for Class I protostars and flat-spectrum sources. The median total luminosity for Class 0 protostars is 5.5 L_⊙, compared to 2.0 L_⊙ for Class I protostars and 3.0 L_⊙ for flat-spectrum sources. Both Class 0 and I protostars cover close to the whole range of L_tot values in the model grid (∼0.06–600 L_⊙), while flat-spectrum sources span a more limited range, from 0.1 to 316 L_⊙.

Zoom In Zoom Out Reset image size

The distribution of centrifugal radii for the whole sample showed a preference for 5 au (see Figure 26). When separating the best-fit disk radii by protostellar class (Figure 37), it is clear that the trend for small centrifugal radii is driven by the flat-spectrum sources and also Class I protostars. The fraction of Class 0 protostars with R_disk = 5 au is 17%; it increases to 46% and 73% for Class I protostars and flat-spectrum sources, respectively. The median disk radius decreases from 100 au for Class 0 protostars to 50 au for Class I protostars to 5 au for flat-spectrum sources. All three histograms are significantly different from one another (K-S test significance level ≪0.001). The unexpectedly small centrifugal radii for Class I protostars and flat-spectrum sources could point to parameter degeneracies (see Section 7.2) or the need to revise certain model assumptions (see Section 7.4).

Zoom In Zoom Out Reset image size

Finally, the distribution of best-fit foreground extinction values (Figure 38) is similar for all three object classes (K-S test significance level >0.03). Even the median values are close: A_V = 9.2 for Class 0 protostars, A_V = 8.9 for Class I protostars, and A_V = 10.1 for flat-spectrum sources. Most objects are fit with relatively low foreground extinction values. As can be seen from Figure 39, the majority of protostars have best-fit A_V values well below the maximum A_V values determined from column density maps, which were used as the largest allowed A_V values for the SED fitter. The ratio of model-derived A_V to observationally constrained maximum A_V is lower than 0.5 for about 60% of the sample.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 40, we plot the reference densities ρ₁₀₀₀ versus the foreground extinction for Class 0, Class I, and flat-spectrum sources. As was already seen in Figure 38, the extinction along the line of sight is similar for all three classes, with most objects in the A_V ∼ 0–30 regime. Class 0 protostars, which have higher envelope densities, tend to have lower A_V values from foreground extinction; the highest-density envelopes are spread among a wide range of A_V values. The result is similar for Class I protostars. Flat-spectrum sources display a range in envelope densities at various foreground extinction values; the lowest-density envelopes typically have A_V < 20. Thus, foreground extinction does not seem to affect the classification of protostars. This result is also supported by the statistical analysis of Stutz & Kainulainen (2015), who found that, for A_V values up to 35, the misclassification of a Class I protostar as a Class 0 protostar due to foreground extinction (which results in a lower T_bol) is low.

Zoom In Zoom Out Reset image size

We found differences in the best-fit envelope densities and inclination angles for the various protostellar classes. The result that Class 0 protostars tend to have larger inclination angles and envelope densities compared to Class I and flat-spectrum objects can also be seen in Figure 41. There are very few Class 0 protostars with low inclination angles; most have relatively high density and i > 60°. Class I protostars are best fit by somewhat lower inclination angles than Class 0 protostars and also lower ρ₁₀₀₀ values. The best-fit reference density for Class I protostars decreases as the inclination angle increases; thus, higher-density protostars are typically classified as Class I protostars only if they are not seen at close to edge-on orientations. Flat-spectrum sources are spread out in density–inclination space, but intermediate inclination angles and low envelope densities are common. There is a relatively large number of objects at i = 18° and a deficit of objects at high inclination angles. The highest-density flat-spectrum sources are seen at inclination angles <50°, while the lower-density objects cover almost the full range of inclination angles.

Zoom In Zoom Out Reset image size

The median parameter values we determined from the best fits for the Class 0, Class I, and flat-spectrum sources (see Table 4) can be used to show representative median SEDs for each protostellar class. In Figure 42, we show model SEDs whose parameter values are equal to the median values found for each of the three protostellar classes. It is apparent that the large envelope density and higher inclination angle for Class 0 protostars cause a deep absorption feature at 10 μm and a steeply rising SED in the mid- and far-IR, with a peak close to 100 μm. In Class I protostars, the SED is less steep and peaks at a shorter wavelength than the median SED of Class 0 protostars. Flat-spectrum sources show the strongest near-IR emission of the three protostellar classes; their median SED is very flat out to 70 μm, but at longer wavelengths it is very similar in shape and flux level to that of Class I protostars.

Zoom In Zoom Out Reset image size

Given that the R values are a measure of the goodness of fit in units of the fractional uncertainty, we can use models that lie within a certain range of the best-fit R value to estimate uncertainties for the various model parameters. For each modeled HOPS target, we tabulated the model parameters for all those models that lie within a difference of 0.5, 1.0, 1.5, and 2.0 of the best-fit R. We then computed the mode (i.e., the value with the highest frequency) for the inclination angle, total luminosity, ρ₁₀₀₀, cavity opening angle, outer disk radius, and A_V in each of the ΔR bins for each object. For any given protostar, the models in each ΔR bin span certain ranges in parameter values; while the modes do not capture the full extent of these ranges, they convey the most common value within each parameter range. The farther away a mode is from the best-fit value, the more poorly constrained the model parameter. Conversely, if a mode of a certain parameter is close to or matches the best-fit value, especially for ΔR = 1.5 or 2, that particular model parameter is well-constrained. In Figures 43–48 we plot the mode versus the best-fit value for six model parameters and four ΔR bins for all 330 targets in our sample. The larger ΔR, the larger the spread in modes is expected to be for each parameter value.

For example, Figure 43 shows that even when considering all models with an R value of up to 2 larger than the best-fit R (ΔR = 2), the inclination angle for objects with a best-fit i of 18° is well-constrained; most modes lie at i = 18°, too, and only a few modes can be found at larger inclination angles. However, objects with best-fit i values of 32° or 41° typically can also be fit by models with lower inclination angles (the majority of modes lies below the line where mode and best-fit value are equal). Inclination angles ≳63° are better constrained, since their modes mostly lie at high inclination angle values, but there are protostars with modes of i = 18°, too.

Zoom In Zoom Out Reset image size

The modes for the total luminosity (Figure 44) show a small spread for models within ΔR = 0.5, but the spread increases as R increases, with some objects displaying up to an order of magnitude in variation of L_tot. As illustrated in Figure 45, the reference density ρ₁₀₀₀ is usually well-constrained; however, as R increases, the modes of the ρ₁₀₀₀ values are often lower than the best-fit values. For the cavity opening angle (Figure 46), many models up to ΔR = 2 have modes of θ = 45°, independent of the best-fit value. Similarly for the centrifugal radius (Figure 47), R_disk = 500 au is a common value. For all four disk radii, the modes tend to be larger than the best-fit values; in particular, objects with a best-fit R_disk of 5 au have a very uncertain disk radius. In general, it looks like our models do not constrain the disk radius and cavity opening angle well. The foreground extinction (Figure 48) displays a certain range of modes for each best-fit value, but objects with A_V ≲ 20 typically have more reliable A_V values from their model fits.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 43–48 allow us to gauge general trends between best-fit values and modes for different model parameters. For results on individual objects, we refer to Appendix B, where we show plots of the difference between the modes and the best-fit values of the major model parameters for all modeled HOPS targets. In this way it is possible to estimate which models are better constrained and thus which objects have more reliable SED fits. In addition, in Appendix B we also include contour plots of R values for different pairs of model parameters for a few targets to illustrate typical parameter degeneracies, which also contribute to parameter uncertainties.

We compared a grid of TSC models to each target by ranking the models using a statistic, R, which measures the deviation between observed and model fluxes in logarithmic space. We did not model each source by further adjusting the model parameters, but instead identified the best-fit SED from our model grid. Thus, we are bound by the range and sampling of parameters chosen for the grid, and while we constructed the grid with the aim of covering the typical parameter space for protostars, it is limited to discrete values. It is likely that many protostars have best-fit parameters that would fall between those sampled by the grid, and a few objects could have parameter values that lie beyond the limits set by the grid. In addition, TSC models are axisymmetric and have mostly smooth density and temperature profiles, and they do not include external heating. They assume a rotating, infalling envelope with constant infall rate and with the gravitational force dominated by the central protostar, but the true envelope structure is likely more complex. The models would not apply to the collapse of a cloud in an initial filamentary or sheet-like geometry or to multiple systems with, e.g., more than one outflow cavity (e.g., Hartmann et al. 1996; Tobin et al. 2012).

Despite the relatively simple models that we use, many of the observed SEDs are fit remarkably well: 75% of the fits have R < 4. In those cases, the continuum traced by the IRS spectrum, the silicate absorption feature at 10 μm, and the PACS fluxes are all accurately reproduced by the model. Even many flat-spectrum sources, which often do not display any spectral features in the mid-infrared and have an overall flat SED out to 30 or 70 μm, often have models that fit them very well. About 75% of Class I protostars and ∼70% of Class 0 protostars have R < 4, while close to 90% of flat-spectrum sources have R values in this range. This validates the choice of parameter values for our model grid. Additional constraints, like limits on foreground extinction or information on the inclination and cavity opening angles from scattered light images or mapping of outflows, would allow us to further test and refine the models. We have used limits on the extinction in our analysis. Although Hubble Space Telescope (HST) scattered light images have been used to constrain models for one HOPS protostar (Fischer et al. 2014), scattered light images are not available for many of our targets. We therefore chose to focus on fitting the SEDs of all of our targets in a uniform way to a well-defined set of models. Future studies will incorporate scattered light images and compare the results to those from the SED fits (J. Booker et al. 2016, in preparation).

The best-fit models from our grid for the HOPS targets both reproduce the SEDs and yield estimates for their protostellar parameters, mostly envelope properties. However, these are not necessarily unique fits to the data for three reasons. First, there are degeneracies in the model parameters; increasing the envelope density or inclination angle, or decreasing the cavity opening angle or disk radius, results in a steeper mid-IR SED slope and deeper silicate features. Each of these parameters affects the SED differently (just the general trends are the same), and the best fit for each object tries optimizing them. The next best fit, however, could be a different combination of these parameters, especially if the SED is not well-constrained by observations (see Section 6.4). Second, although the TSC models reproduce the observed SEDs, other models with different envelope geometries may also be able to fit the same SEDs. The modeling presented here is only valid in the context of the TSC models of single stars, and the resulting derived properties are only valid within that framework. Last, neglecting external heating could result in overestimated envelope densities and luminosities, with the most noticeable effects (ρ₁₀₀₀ and L_tot too large by factors of a few) on low-luminosity sources exposed to strong radiation fields (see Section 4.3). From the distribution of best-fit L_tot values, we estimate that ∼20% of HOPS targets in our sample could be affected by external heating. Even though we do not know the strength of external heating for each protostar, it is likely that external heating would only result in relatively small changes in the derived envelope parameters for these protostars.

When comparing envelope parameters sorted by SED classes, we found that envelope densities and inclination angles decrease from the sample of Class 0 protostars through that of Class I protostars to that of flat-spectrum objects. The former is likely an evolutionary effect, while the latter confirms the results of previous work (e.g., Evans et al. 2009) that the inclination angle has an important effect on the SED and that the evolutionary state of an object cannot be derived from SED slopes alone. Thus, there is a difference between the "stage" and "class" of an object (Robitaille et al. 2006); Stage 0 and I objects are characterized by substantial envelopes, Stage II objects are surrounded by optically thick disks, with possibly some remnant infalling envelopes, and Stage III objects have optically thin disks.

In general, the trends we see among model parameters are a consequence of the definition of a protostar based on its SED: in order to be classified as a Class 0 or I object, a protostar is required to have a near- to mid-infrared SED slope larger than 0.3. A protostellar model with a small cavity opening angle, small centrifugal radius, and/or high inclination angle will generate such an SED, since it increases the optical depth along the line of sight. Models with a large cavity will only yield a rising SED in the 2–40 μm spectral range if their envelope density is large or the inclination angle is relatively high.

We find that Class 0 protostars can be best fit not only by very high envelope densities but also moderately high envelope densities and large inclination angles. The bolometric temperature, T_bol, which is used to separate Class 0 from Class I protostars, is inclination dependent; some Class I protostars are shifted to Class 0 protostars if they are viewed more edge-on. The higher-density Class I protostars tend to have lower inclination angles (but still >50°); thus, their evolutionary stage could be similar to more embedded protostars that are seen edge-on and classified as Class 0 protostars. Conversely, some Class 0 protostars with large inclination angles, but lower envelope densities, could be in a later evolutionary stage than typical Class 0 protostars. Similarly, Class I protostars with large i and low ρ₁₀₀₀ values could be edge-on Stage II objects (whose infrared emission is dominated by a disk). Finally, low-inclination Stage 0 and Stage I protostars can appear as flat-spectrum sources (Calvet et al. 1994).

Nevertheless, the observed trend in envelope densities suggests that the variations in the observed SEDs track, in great part, an evolution toward lower envelope densities and lower infall rates. Assuming a certain mass for the central star, the reference density in our models can be used to infer an envelope infall rate (${\dot{M}}_{\mathrm{env}}\propto {\rho }_{1000}\sqrt{{M}_{*}}$). As mentioned in Section 4.1, this infall rate is model dependent and therefore tied to the assumptions of the models. With this in mind, the median ρ₁₀₀₀ values for the Class 0, Class I, and flat-spectrum protostars in our sample correspond to envelope infall rates of 2.5 × 10⁻⁵, 1.0 × 10⁻⁶, and 5.0 × 10⁻⁷ M_⊙ yr⁻¹, respectively, for a 0.5 M_⊙ star. Using a more realistic assumption of larger stellar mass for more evolved protostars, the infall rates for Class I and flat-spectrum protostars would be larger than these values by a factor of a few. However, just larger stellar masses cannot explain the large decrease of a factor of 50 in the median envelope density from Class 0 to flat-spectrum protostars; to achieve such a decrease with a constant infall rate of 2.5 × 10⁻⁵ M_⊙ yr⁻¹, the stellar mass would have to increase by a factor of 2500. Thus, within the context of our model fits, we can conclude that, as envelopes become more tenuous, the infall rates also decrease.

Other trends are also apparent. Class 0 protostars and flat-spectrum sources show a relatively flat distribution of cavity opening angles. On the other hand, the best fit for a large fraction of Class I protostars (40%) results in θ = 5°. This could point to a degeneracy in the models, since protostars with small cavity opening angles tend to have lower envelope densities (and also lower total luminosities); thus, the smaller cavity partly compensates for the lower opacity resulting from the lower envelope density (see also Figure 49).

Zoom In Zoom Out Reset image size

Even though our models do not yield reliable disk properties, we can make a statement about the difference in the best-fit centrifugal radii (or R_disk), which are tied to the structure of the rotating envelope given by the model fits. It should be noted that the centrifugal radii set a lower limit to the disk radii, since disks may spread outward owing to viscous accretion. Most Class I protostars and flat-spectrum sources are fit with a centrifugal radius of just 5 au. Since the smallest centrifugal radius in our model grid is 5 au and the next value is 50 au, we can state that, except for Class 0 protostars, most protostars in our sample have R_disk < 50 au, and some might even have R_disk < 5 au.

Small disks of those sizes have been observed; radio interferometry of the multiple protostellar system L1551 IRS 5 shows disks whose semimajor axes are ≲20 au (Rodríguez et al. 1998; Lim & Takakuwa 2006). However, there is a degeneracy between the centrifugal radius and the envelope density; for protostars with low envelope densities, the small centrifugal radius can compensate the decrease in opacity by concentrating more material closer to the star. As can be seen in Figure 50, most protostars with R_disk = 5 au also have lower envelope densities. Inclination angle also plays a role; protostars seen at i > 80° typically have larger centrifugal radii. In addition, our envelope models include outflow cavities, which allow some of the mid-IR radiation to escape. In order to generate model SEDs with silicate absorption features and steep mid-IR slopes at low to intermediate inclination angles, a lower R_disk value is needed. We will discuss the potential implications of the small cavity sizes and R_c values for Class I protostars and flat-spectrum sources in Section 7.4.

Zoom In Zoom Out Reset image size

7.2.1. The Most Embedded Protostars

7.2.2. Flat-spectrum Sources

A particularly interesting group of protostars that are not easy to categorize are the flat-spectrum sources. They are thought to include objects in transition between Stages I and II, when the envelope is being dispersed (Greene et al. 1994). In particular those with low envelope densities could be more evolved protostars, or they could be protostars that started out with more tenuous envelopes. On the other hand, flat-spectrum sources could also be highly inclined disk sources (see Crapsi et al. 2008), or protostars surrounded by dense envelopes, but seen close to face-on (Calvet et al. 1994). This type of misclassification could have a large effect on the lifetimes of the earlier protostellar stages and thus on the timeline of envelope dispersal. Among the 330 HOPS targets in our sample, we identified 102 flat-spectrum sources based on their flat (−0.3 to +0.3) spectral index from 4.5 to 24 μm (or 5–25 μm in a few cases). Thus, they compose a fairly large fraction of our protostellar sample. Of these 102 objects, 47 have a negative spectral index and 55 have one between 0 and +0.3; 41 have a spectral index between −0.1 and 0.1, which results in a very flat mid-infrared SED.

Despite a flat SED slope between 4.5 and 24 μm, many flat-spectrum sources display a weak silicate emission or absorption feature at 10 μm, which may indicate the presence of a very tenuous infalling envelope or may be the result of the viewing geometry. Some SEDs are very flat out to 100 μm, others have negative spectral slopes beyond 40 μm, and again others a rising SED from the mid- to the far-IR. There are also objects with more pronounced absorption features due to not only silicates but also ices, as are typically found in Class 0 and I protostars, but also edge-on disks (see HOPS 82, 85, 89, 90, 92, 129, 150, 200, 210, 211, 281, 304, 331, and 363). Only two flat-spectrum sources have prominent silicate emission features, and their SEDs are reminiscent of protoplanetary disks (see HOPS 187 and 199). Thus, flat-spectrum sources likely include objects of a variety of evolutionary stages.

The latter conclusion can also be drawn when analyzing the distribution of envelope reference densities and inclination angles for flat-spectrum sources. In Figure 41, we showed that flat-spectrum sources typically have intermediate inclination angles and lower envelope densities. To compare their properties more directly to Class 0 and I protostars, in Figure 51 we show the median best-fit ρ₁₀₀₀ value at each best-fit inclination angle; it is larger for Class 0 and I protostars than for flat-spectrum sources at all inclination angles. For Class 0 and I protostars, the median ρ₁₀₀₀ value is highest at intermediate inclination angles, decreases at larger inclination angles, and then increases again for i > 80°. For flat-spectrum sources, the median ρ₁₀₀₀ value is relatively flat over the 18°–63° region but has its peak value at i = 41°; it decreases for larger inclination angles. The only flat-spectrum source with a best-fit inclination angle of 81°, HOPS 357, has a very low envelope density (the lowest value for this parameter in the model grid), and its spectrum displays a deep silicate absorption feature.

Zoom In Zoom Out Reset image size

Overall, this shows that, while a range of envelope densities and inclination angles can explain flat-spectrum sources, their envelope densities are typically lower than for Class 0 and I protostars. The higher-density objects are seen at low to intermediate inclination angles, while only the lowest-density objects are seen closer to edge-on. Some of the high-density flat-spectrum sources could actually be more embedded protostars (Stage 0 objects) seen face-on (which would be classified as Class 0 objects if seen at larger inclination angles). Thus, in terms of envelope evolution, they include a diverse group of objects.

We note that even though we find that flat-spectrum sources have in general lower envelope densities than Class 0 and Class I objects, their best fit does include an envelope in almost all cases; just 3 of the 102 flat-spectrum sources are best fit without an envelope. This seems to contradict recent findings by Heiderman & Evans (2015), who found that only about 50% of flat-spectrum sources were actually protostars surrounded by envelopes. This could be partly explained by different criteria used to select flat-spectrum sources; in the Heiderman & Evans (2015) sample, flat-spectrum sources are selected by their extinction-corrected 2–24 μm spectral index (see also Evans et al. 2009; Dunham et al. 2013), while our sample uses a flat 4.5–24 μm spectral index. Moreover, in their study Heiderman & Evans (2015) detected the presence of an envelope via HCO⁺ emission, and they found that almost all sources detected in the submillimeter are also detected in HCO⁺ (but the opposite does not always hold). For our sample of Orion protostars, we find that 75% of Class 0+I protostars observed with SABOCA (350 μm) are detected, while only 47% of flat-spectrum sources have detections. For LABOCA observations (870 μm), these two fractions amount to 41% and 21%, respectively. Thus, we find that flat-spectrum sources have a ∼50% lower submillimeter detection rate than Class 0+I protostars. Flat-spectrum sources without submillimeter detections would likely also not display HCO⁺ emission and thus would be considered as protostars without envelopes by Heiderman & Evans (2015).

To compare how our submillimeter detections correlate with the presence of an envelope, in Figures 52 and 53 we show the derived best-fit reference envelope densities as a function of 350 or 870 μm fluxes for the combined Class 0+I sample and the flat-spectrum sources. We also differentiate the distribution of envelope densities between measured flux values and upper limits; at 870 μm, the upper limits are often cases where the sources are not detected due to confusion with bright, spatially varying emission. We find that even protostars with upper limits at 350 and 870 μm are best fit with an envelope; however, the envelope density is lower for objects with upper limits in the submillimeter. This is especially evident for Class 0+I protostars; for flat-spectrum sources, the distributions of envelope densities for submillimeter detections and upper limits show significant overlap. Four times as many flat-spectrum sources have upper limits instead of detections at 870 μm, but their derived ρ₁₀₀₀ values span almost the full range of values. Furthermore, the median ρ₁₀₀₀ value of 1.19 × 10⁻¹⁹ g cm⁻³ for sources without detections is relatively close to the median value of 5.95 × 10⁻¹⁹ g cm⁻³ for the sources with 870 μm detections. Thus, our model fits do not rely on submillimeter detections to yield a best fit with an envelope; in most cases the near- to far-IR SED is sufficient to constrain the properties of the envelope.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

7.2.3. Sources without an Envelope and Class II Objects

The luminosity distribution of protostars is a significant constraint on protostellar evolution, and it is important to understand the effect of the envelope on the observed luminosity (e.g., Offner & McKee 2011). The bolometric luminosity distribution of the HOPS protostars is very similar to that determined for the Spitzer-identified protostars by Kryukova et al. (2012) with a peak near 1 L_⊙ (Figure 2). In contrast, the distribution of the total luminosities from the models shows a peak near 2.5 L_⊙ (Figure 36), indicating that the luminosities of protostars may be systematically underestimated by the bolometric luminosities, which do not take into account the inclination angle (and thus beaming of the radiation along the outflow cavities), as well as foreground extinction (see Figure 32 in Section 6.2).

Higher intrinsic luminosities for protostars could help address the "luminosity problem" first pointed out by Kenyon et al. (1990), who found that the luminosities of protostars are lower by about an order of magnitude than a simple estimate of the expected accretion luminosity. However, an increase in the luminosity by a factor of 2.5–3 would not solve the problem; solutions proposed by other authors, such as mass-dependent accretion rates (Offner & McKee 2011) or episodic accretion events (Dunham & Vorobyov 2012), are still needed.

Our best-fit models also suggest that Class 0 protostars have a different distribution of L_tot values compared to Class I protostars or flat-spectrum sources. Their median total luminosity is higher, which could be an indication of larger accretion luminosities for younger protostars. We must bear in mind the caveats and degeneracies mentioned above; in particular, in some cases the higher luminosity could be related to the adoption of an overly large inclination angle, which results in most of the emitted radiation not reaching the observer. Nevertheless, these differences have potentially important implications for protostellar evolution, which will be discussed in a future publication (W. Fischer et al. 2016, in preparation).

Although the TSC models provide impressive fits to the SEDs, some of the observed trends suggest problems with the models. First, the distribution of inclination angles (Figure 29) deviates from what we expect from a randomly oriented sample of protostars. Although this could result from unintentional selection biases in our sample of protostars, it may also be the effect of applying the wrong envelope model to the data.

Furthermore, our data show flat distributions in cavity open angles for Class 0 and flat-spectrum sources, but an excess of small cavities for the Class I protostars (Figure 35). We also find that protostars with large cavities often have high envelope densities (Figure 49). For example, models with high envelope densities viewed more edge-on require large cavity opening angles and high L_tot values to generate sufficient mid-IR flux; this is the case for a few of our highest-luminosity objects (HOPS 87, 108, and 178). These trends do not support the notion of increasing cavity size with later evolutionary stage, which would be expected if outflows play a major role in dispersing envelopes (Arce & Sargent 2006). This may suggest that cavity sizes are not growing with time; however, this may also imply a deviation from spherical symmetry for the initial configuration of the collapsing envelopes. Such a deviation may result if the envelope collapses from the fragmentation of a flattened sheet or elongated filament.

Finally, we find an excess of small values of R_disk, and therefore small centrifugal radii, for Class I and flat-spectrum protostars (Figure 37). This is contrary to the expectation from the TSC model, in which the late stages of protostellar evolution are characterized by the infall of high angular momentum material from large radii and hence larger values of R_c. This may imply that disks sizes are small, but it may also be the result of incorrect assumptions about the distribution of angular momentum in the TSC model.

In total, these "conundrums" that arise from our model fits hint that the current models do not realistically reproduce the structure of collapsing envelopes. Future high-resolution observations at submillimeter and longer wavelengths that resolve the structure and motions of envelopes may provide the means to develop more refined models that can fit the SEDs with more realistic envelope configurations.

Table 6 lists the newest HOPS targets discovered in Herschel PACS data (HOPS 394 to 409; Stutz et al. 2013; Tobin et al. 2015) with their coordinates and original identifiers from Stutz et al. (2013).

C.2.1. Young Stellar Objects

Among our HOPS sample, there are 41 targets that are likely YSOs, but they lack PACS measurements at 70 and 160 μm and were therefore not included in the modeling sample. There are four additional targets with HOPS numbers that were not modeled; they are HOPS 109, 111, 212, and 362, and they are duplicates of HOPS 40, 60, 211, and 169, respectively. Table 7 lists the 41 likely YSOs in the HOPS catalog that were not part of the modeling sample; their SEDs are shown in Figure Set 57. Among them, 17 were not observed by PACS at 70 μm, while 24 were observed, but not detected at 70 μm. The majority of these targets are Class I protostars or flat-spectrum sources; only one is a Class 0 protostar, and five are Class II pre-main-sequence stars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (3 panels)

Tarball of figure set images

Table 7.  YSOs in the HOPS Sample with No PACS Data

Object	R.A.	Decl.	Class	Lbol	Tbol	n4.5–24	PACS Flag
	(°)	(°)		(L⊙)	(K)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
HOPS 0	88.6171	1.6264	I	0.011	652.2	0.514	−1
HOPS 8	83.8880	−5.9851	I	0.013	329.8	0.419	0
HOPS 9	83.9550	−5.9843	I	0.006	281.5	0.858	−1
HOPS 14	84.0799	−5.9251	flat	0.042	464.0	0.246	0
HOPS 23	84.0745	−5.7818	I	0.012	346.8	0.539	0
HOPS 25	83.8443	−5.7415	flat	0.045	646.6	0.165	0
HOPS 31	83.8219	−5.6741	flat	0.024	634.7	0.304	0
HOPS 34	83.7954	−5.6585	I	0.013	235.5	0.762	0
HOPS 35	83.8331	−5.6503	I	0.044	305.2	0.884	0
HOPS 37	83.6986	−5.6237	flat	0.016	913.4	0.230	0
HOPS 51	83.8160	−5.5015	II	0.518	130.2	⋯	0
HOPS 52	83.8180	−5.4924	flat	0.641	610.8	−0.163	0
HOPS 54	83.3437	−5.3841	II	0.097	1879.3	−0.37	−1
HOPS 62	83.8524	−5.1916	flat	0.660	1154.1	0.043	0
HOPS 63	83.8538	−5.1671	flat	0.516	544.5	0.004	0
HOPS 64	83.8625	−5.1650	I	15.347	29.7	0.503	0
HOPS 69	83.8551	−5.1400	flat	2.778	31.3	−0.189	0
HOPS 79	83.8662	−5.0934	flat	0.086	666.2	−0.137	0
HOPS 103	83.5508	−4.8353	flat	0.142	1484.3	−0.032	0
HOPS 104	83.7782	−4.8338	I	0.044	337.3	0.837	0
HOPS 110	84.0093	−5.0472	I	0.014	244.0	⋯	−1
HOPS 126	85.0408	−7.1650	flat	0.132	1865.3	−0.136	−1
HOPS 151	84.6787	−6.9447	II	0.061	799.4	−0.505	0
HOPS 155	84.3160	−7.2972	flat	0.013	393.8	0.133	−1
HOPS 162	84.1291	−6.8780	II	0.015	909.9	0.352	−1
HOPS 180	84.2475	−6.1710	II	0.011	1493.5	0.578	−1
HOPS 195	84.0002	−6.1206	flat	0.032	659.7	0.399	0
HOPS 217	85.7965	−8.4056	I	0.008	323.8	0.773	−1
HOPS 230	85.6283	−8.1515	flat	0.267	1260.2	−0.104	−1
HOPS 231	85.1189	−8.5486	flat	0.024	386.0	−0.239	−1
HOPS 269	85.3625	−7.7094	flat	0.025	230.2	0.023	−1
HOPS 289	84.9865	−7.5017	I	0.095	331.1	0.868	0
HOPS 296	85.3215	−2.3021	I	0.022	326.0	0.931	0
HOPS 302	85.0934	−2.2610	flat	0.383	1367.2	−0.032	−1
HOPS 307	85.3077	−1.7844	0	0.748	57.1	1.506	−1
HOPS 314	86.6505	−0.3414	I	0.015	276.2	1.112	−1
HOPS 327	86.6139	0.1477	flat	0.020	991.0	0.145	−1
HOPS 328	86.5561	0.1759	I	0.012	326.3	0.868	−1
HOPS 330	86.7140	0.3298	flat	0.121	385.2	0.285	0
HOPS 332	86.8821	0.3391	flat	0.249	145.5	0.045	0
HOPS 360	86.8629	0.3425	I	1.017	43.2	⋯	0

Note. Column (1) lists the HOPS number of the object, columns (2) and (3) its J2000 coordinates in degrees, column (4) the type based on SED classification, column (5) the bolometric luminosity, column (6) the bolometric temperature, column (7) the 4.5–24 μm SED slope, and column (8) a flag identifying whether the object was not observed by PACS (flag value of −1) or not detected by PACS at 70 μm (flag value of 0).

Download table as:  ASCIITypeset image

Most of these YSOs have very faint fluxes in the near- to mid-IR. They could be deeply embedded protostars, like HOPS 307, or just very low mass protostars with weak envelope emission. Objects with little excess emission out to about 8 μm, a 10 μm silicate emission feature, and a more or less steeply rising SED beyond 12 μm are likely transitional disks (see Kim et al. 2013); good examples are HOPS 51 and 54. It is possible that some of the YSOs in this sample are actually extragalactic contaminants, in particular objects with flat SEDs (see more about this subset of our HOPS sample in the next subsection).

C.2.2. Contaminants

Our HOPS sample contains 29 targets that turned out to be likely extragalactic contaminants. These objects are listed in Table 8, and their SEDs are shown in Figure Set 58. Most galaxies were identified based on the presence of PAH features or emission lines in their IRS spectra (in particular, the 5–15 μm region), the absence of a silicate absorption feature at 10 μm, and an overall flat or slightly rising mid-infrared continuum. Clear examples are HOPS 21, 27, 46, 48, 55, 61, 72, 106, 161, and 301.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (2 panels)

Tarball of figure set images

Table 8.  Likely Extragalactic Contaminants in the HOPS Sample

Object	R.A.	decl.	Lbol	Tbol	n4.5–24
	(°)	(°)	(L⊙)	(K)
(1)	(2)	(3)	(4)	(5)	(6)
HOPS 21	84.0421	−5.8356	0.074	584.5	0.824
HOPS 27	84.0905	−5.6995	0.102	118.2	0.458
HOPS 39	84.0934	−5.6069	0.100	159.4	1.010
HOPS 46	83.6758	−5.5509	0.141	1081.9	−0.083
HOPS 48	83.7773	−5.5477	0.608	611.0	0.252
HOPS 55	83.4754	−5.3638	0.316	101.5	1.732
HOPS 61	83.3579	−5.2007	0.045	721.4	−0.031
HOPS 67	83.8445	−5.1428	0.044	278.7	0.900
HOPS 72	83.8571	−5.1296	0.545	693.0	0.068
HOPS 83	83.9822	−5.0771	0.131	293.9	−0.119
HOPS 97	83.8704	−4.9608	0.190	403.8	0.125
HOPS 101	83.7843	−4.9027	3.355	481.2	−0.207
HOPS 106	84.0518	−4.7544	0.016	359.7	0.943
HOPS 112	85.1833	−7.3786	0.014	390.2	0.093
HOPS 146	84.6840	−7.0112	0.053	519.7	−0.178
HOPS 161	84.1448	−7.1871	0.062	179.1	0.458
HOPS 196	83.8371	−6.3062	0.065	165.5	−0.042
HOPS 202	83.4330	−6.2295	0.018	736.3	0.534
HOPS 205	85.7620	−8.7971	0.067	427.8	0.549
HOPS 218	85.7912	−8.2232	0.007	332.5	0.653
HOPS 292	84.4787	−7.6890	0.068	280.5	0.421
HOPS 301	85.4366	−2.2654	2.955	518.8	0.271
HOPS 306	85.7630	−1.8013	0.038	310.8	0.060
HOPS 308	85.8082	−1.7195	0.042	156.0	0.791
HOPS 309	85.6973	−1.4131	0.041	111.9	−0.168
HOPS 313	85.2532	−1.1529	0.030	154.1	1.039
HOPS 339	86.4733	0.4243	0.128	398.3	0.246
HOPS 348	86.7511	0.3438	0.286	84.0	⋯
HOPS 351	83.8809	−5.0797	0.016	217.1	⋯

Note. Column (1) lists the HOPS number of the object, columns (2) and (3) its J2000 coordinates in degrees, column (4) the bolometric luminosity, column (5) the bolometric temperature, and column (6) the 4.5–24 μm SED slope.

Download table as:  ASCIITypeset image

Two objects, HOPS 202 and HOPS 205, have a tentative 10 μm silicate emission feature and a steep rise of their SED beyond 12 μm, but also PAH emission features in their IRS spectrum; they could be transitional disks and not galaxies.

In some cases, e.g., HOPS 308, targets classified as extragalactic contaminants are very faint in the near- to mid-infrared; instead of galaxies, they could be very low mass or deeply embedded protostars.

The mid-IR SED of HOPS 339 is mostly flat but displays a sharp 10 μm absorption feature; based on the SED alone, it would not necessarily be classified as a galaxy, but high-resolution near-IR HST images resolve its extended emission and reveal a spiral galaxy (J. Booker et al. 2016, in preparation).

Finally, HOPS 349, 350, 352, 353, 356, and 381 have poorly sampled SEDs (just one or two flux measurements), so their nature is quite uncertain. We list their coordinates in Table 9.

Table 9.  Targets in the HOPS Sample with Uncertain Nature

Object	R.A.	Decl.
(1)	(2)	(3)
HOPS 349	83.8592	−5.1426
HOPS 350	83.8758	−5.1386
HOPS 352	83.8617	−5.0675
HOPS 353	88.5556	1.7175
HOPS 356	85.5341	−1.4438
HOPS 381	83.7816	−5.6986

Note. Column (1) lists the HOPS number of the object, columns (2) and (3) its J2000 coordinates in degrees.

Download table as:  ASCIITypeset image