An HST Survey of Protostellar Outflow Cavities: Does Feedback Clear Envelopes?

Low-mass protostars are characterized by a rapid evolution, with the accretion of the stellar mass, the formation of disks, and potentially the initiation of planet formation occurring within 0.5 Myr (Cassen & Moosman 1981; Arce & Sargent 2006; Dunham et al. 2014; ALMA Partnership et al. 2015; Dipierro et al. 2015). The defining characteristic of the protostellar phase is the presence of a dusty, infalling envelope that absorbs and reprocesses most of the luminosity from the central protostar. In the initial phases of protostellar evolution, the envelope dominates the mass, while in the later phases, most of the mass is already accreted onto the star. Even in these later phases, the mass of the envelope typically exceeds that of the circumstellar disks surrounding the central protostar (e.g., Fischer et al. 2014); hence, infall in these phases shapes the properties of circumstellar disks and sets the stage for planet formation. Understanding the factors that govern the evolution of the envelopes—and thereby influence mass accretion and the properties of nascent disks—is a key problem in star and planet formation studies.

This evolution is accompanied by a rapid change in the shape of the spectral energy distributions (SEDs) produced by the reprocessing and scattering of radiative energy in the evolving disks and envelopes (Furlan et al. 2016). Since the central protostar is deeply embedded in its envelope, the effective temperatures and photospheric luminosities of protostars typically cannot be measured directly. In most cases, unlike pre-main-sequence stars, they cannot be reliably placed on HR diagrams and compared to evolutionary tracks to estimate masses and ages. Instead, the evolution of protostars is largely inferred from the shape of their SEDs. This evolution is typically measured by sorting protostars into bulk evolutionary classes based on the percentage of luminosity radiated in the submillimeter, their near- to mid-infrared spectral index or T_bol, the bolometric temperature (e.g., Adams & Shu 1985; Andre et al. 1993; Myers & Ladd 1993; Stutz et al. 2013; Dunham et al. 2014; Furlan et al. 2016). The observed sequence of evolutionary classes, Class 0, Class I, and flat-spectrum, shows the peak of the SED shifting from the far-infrared to the mid-infrared and the SED flattening as the protostars evolve and the envelopes dissipate (e.g., Furlan et al. 2016). Class II objects can be identified by their decreasing near- to mid-IR SED slopes and are primarily pre-main-sequence stars with disks that have exited the protostellar phase.

Due to the flattening of envelopes by rotation and the clearing of cavities in the envelopes by outflows, the luminosity of the protostars is not radiated isotropically, but is preferentially beamed along the rotation axis of the protostars. The resulting SEDs depend on the inclination of the protostars, and the effects of inclination on the SEDs are difficult to disentangle from those due to evolution (Kenyon et al. 1993; Whitney et al. 2003). To circumvent this degeneracy, Whitney & Hartmann (1993) and Robitaille et al. (2007) proposed a set of evolutionary stages that are dependent on the physical properties of the envelopes and not inclination; however, it is often difficult to reliably infer the stage of a protostar from the observed SED alone. Nevertheless, taking into account the uncertainties due to inclination, Furlan et al. (2016) demonstrate that the observed SEDs of the distinct evolutionary classes require the dissipation of the envelope, with the density of the envelope gas (as inferred by model fits to the SEDs) dropping by a factor of 50 between the Class 0 and flat-spectrum phases. This shows that the envelopes decrease dramatically in density during the Class I phase.

Although SEDs are currently the primary information we have on large samples of protostars, imaging at millimeter, submillimeter, and near-infrared wavelengths can be used to study protostellar evolution by resolving structures in the envelope that may change as protostars evolve (e.g., Arce & Sargent 2006). Hubble Space Telescope (HST) near-infrared images of protostars resolve structures seen directly in light scattered by dust grains in an envelope or in silhouette against the scattered light, placing constraints on the envelopes and disks that are complementary to those inferred from SEDs. HST imaging of protostars by Padgett et al. (1999), Allen et al. (2002), Terebey et al. (2006), and Fischer et al. (2014) show outflow cavities illuminated in scattered light, edge-on disks seen in absorption, and shadows cast into the envelopes by flared disks.

Of particular interest is the role of feedback from outflows in driving the evolution of protostars by clearing the envelope and halting infall. SED-based measurements cannot reliably constrain outflow cavity sizes (Furlan et al. 2016); hence, studies of the growth of outflow cavities must rely on observations that spatially resolve structures in envelopes. The CO observations of nine Class 0, I, and II sources by Arce & Sargent (2006) showed a widening in outflow size with evolutionary class. Bolstering their sample by nine sources in the literature, they found evidence that outflow cavity sizes increase progressively as protostars evolve. Tobin et al. (2007) and Seale & Looney (2008) used Spitzer IRAC images of protostar outflow cavities illuminated in scattered light to study the growth of cavities, and the latter authors found some evidence of outflow cavity growth with evolution, although with significant scatter.

These studies suggest that feedback from outflows play a significant role in the decrease or halting of infall and accretion. Although accretion from the disk can continue after infall stops, the resulting increase in mass is small compared to the stellar mass. By reducing or halting infall, feedback can also play an important role in the star formation efficiency inferred from the core mass function. In particular, the mass function of cores identified in submm measurements can reproduce the initial mass function if each core forms a star with a star formation efficiency (defined by the stellar to initial core mass) of ∼30%–40% (Alves et al. 2007; Könyves et al. 2015). Furthermore, simulations of protostars including feedback can produce star formation efficiencies of 50% or lower (Hansen et al. 2012; Machida & Matsumoto 2012; Machida & Hosokawa 2013; Offner & Arce 2014; Offner & Chaban 2017).

There are difficulties, however, in explaining the low star formation efficiency with feedback alone. Single-dish radio observations suggest that outflows may carry too little mass to clear out the envelope in 0.5 Myr (Hatchell et al. 2007; Curtis et al. 2010). Furthermore, even large cavities clear less than half of the envelope mass (Frank et al. 2014). These studies, however, likely underestimate the amount of entrained, low-velocity gas in the outflow (Dunham et al. 2014). ALMA observations can now map these lower-velocity flows and show whether they transport a significant fraction of the envelope gas over the lifetime of the protostar (Zhang et al. 2016).

To further investigate the impact of outflows on protostellar envelopes, we use in this work the largest HST survey of protostars to date. This survey focuses on the sample of protostars targeted by the Herschel Orion Protostar Survey, or HOPS. The protostars were identified using combined 2MASS and Spitzer photometry from the Spitzer Orion Survey (Megeath et al. 2012, 2016) and observed with Herschel and APEX to obtain well-sampled 1.2 μm to 870 μm SEDs. This sample was supplemented by very red protostars discovered with Herschel (Stutz et al. 2013; Tobin et al. 2016). Furlan et al. (2016) published the SEDs of the entire sample and then presented model fits to 319 of the protostars and 11 pre-main-sequence stars after rejecting likely extragalactic contamination and sources without Herschel detections. The HST survey examined 304 of these sources, (enumerated in Table 4), using initially NICMOS at 1.60 and 2.05 μm, and then after the failure of NICMOS, WFC3 at 1.60 μm. A search for binary systems using these data was published by Kounkel et al. (2016).

The morphologies of outflow cavities carved by the outflows can be seen by mapping the location of the cavity wall in scattered light. The volumes of the cavities carved by outflows can then be directly measured. The mechanism for creating these cavities, whether by jet precession, wide-angled winds, or jet entrainment (Raga & Cabrit 1993; Matzner & McKee 1999; Lee et al. 2001; Ybarra et al. 2006) is still debated. Independent of the underlying mechanism, the scattered-light cavities provide a direct measurement of the cleared gas with the 80 au (018) resolution of HST. These are used in this work to estimate the fraction of the volume cleared, which provides an estimate of the fraction of mass cleared.

In Section 2, we discuss the observations used in this paper. We make use of radiative transfer modeling, described in Section 3, to understand the morphologies of the observed cavities and to calibrate the relationship between the scattered-light distributions and the cavity properties. In Section 4, we present the morphologies of the observed protostars and our analysis of the cavity sizes. Finally, we discuss the implications for protostellar evolution in Section 5. Images of the protostars in our sample are shown in Appendix A.

The HST observations were assembled from two GO programs and a snapshot program. The bulk of the sample was observed in program GO 11548. The Near Infrared Camera and Multi-Object Spectrometer's (NICMOS) F205W and F160W filters were used for a total of 87 orbits in August and September of 2008 to image 92 objects in the HOPS catalog, before the failure of the cryocooler of NICMOS. After the 2009 June deployment of the Wide Field Camera 3 (WFC3), 126 orbits were used between 2009 August and 2010 December to observe 237 HOPS objects with the F160W filter. The observation and reduction of these data is described in Kounkel et al. (2016). A subsequent program using WFC3, SNAP 14181, was designed to target multiple star-forming regions within 500 pc. It completed observations during 114 orbits between 2015 December and 2017 September, 10 of which imaged 13 objects in the Orion molecular clouds. A final WFC3 study, program GO 14695, targeted four objects in Orion with weak 24 μm fluxes atypical of protostars. These observations were conducted in 2016 September with four orbits. For these final two programs, we used the standard data products produced from the calwf3 data reduction pipeline, which were then combined with AstroDrizzle from the DrizzlePac package using a drop size of 1 onto a pixel scale of 013.

The NICMOS observations used the NIC2 camera, which has a 0075 pixel size and resolution of 015. Integration times were 1215.4 s and 767.6 s for F160W and F205W filters, respectively. The WFC3 integration times were 2496.2 s for GO 11548, 1596.9 s for SNAP 14181, and 2396.9 s for GO 14695. All have a 018 angular resolution and a pixel size of 013. In this work, we adopt a distance to Orion of 420 parsecs for consistency with Furlan et al. (2016). This is within the range of distances found in Kounkel et al. (2018) and Großschedl et al. (2018) through APOGEE and Gaia measurements. At this distance, both NICMOS and WFC3 resolve structures down to 80 au scales.

Nine images taken with NICMOS are excluded from this analysis, due to the lack of guide star tracking; these contain HOPS 46, 47, 134, 139, 149, 227, 250, 271, and 276. Three WFC3 images, containing HOPS 293, 330, and 336, were also excluded due to what appear to be tracking failures. One additional WFC3 observation was excluded due to an apparent pointing error, with its target object, HOPS 100, only partially appearing on the edge of the frame. Three images where only one guide star was used, those containing HOPS 10, 177, 316, and 358, may suffer from a small amount of rotation during the exposure, although this is not apparent in the data. These are included in our program. Twenty-seven of the HOPS targets were imaged by both NICMOS and WFC3 due to their proximity to other protostars. Of these sources, only HOPS 250 showed a clear difference between the two observations due to the tracking failure.

Some of the HOPS targets were classified as potential extragalactic contamination by Furlan et al. (2016) based on the presence of PAH features in their Spitzer IRS spectrum, lack of silicate absorption at 10 μm, or the shape of the mid-infrared continuum (see appendix of Furlan et al. (2016)). The HST observations provide an independent means for separating galaxies from protostars. Only one source, HOPS 339, is conclusively determined by its morphology to be a galaxy and is omitted from the table in the appendix of this work (Table 4). The WFC3 image of this source is shown in Appendix D. Conversely, we add back into our sample and assign a class to HOPS 48, 67, and 301. These have morphologies in WFC3 imaging indicative of protostellar cavities. The nature of the remaining potentially extragalactic sources could not be clarified through WFC3 imaging. In program GO 14695, two of the four targeted sources were found not to be protostars; one was a galaxy and one an outflow knot; neither of these has a HOPS number (see Appendix D). In total, we imaged 304 objects in our sample. We note that seven of these were determined to be Class II objects by their SEDs in Furlan et al. (2016). Since these sources are in the HOPS sample and may have residual envelopes, we keep them in the analysis. We typically use the term "protostars" to refer to this entire sample. In addition, we serendipitously observed two Class II sources with nebulosity in our images (Kounkel et al. 2016). We describe these objects in Appendix C.

In order to quantify the shape of the observed cavities, we used a monochromatic Monte Carlo radiative transfer code, HO-CHUNK.ttsscat.20090521 (based on Whitney & Hartmann 1992, 1993). With this code, we simulated 1.60 μm images of a half-solar-mass star surrounded by a flared disk with a power-law radial density and scale height, and an envelope following the Terebey, Shu, and Cassen (TSC) model described in Terebey et al. (1984) (see also Ulrich 1976; Cassen & Moosman 1981). We examined six envelope densities (each corresponding to a different mass infall rate in Table 1), five cavity half-opening angles (see Figure 1 for the definition), five disk sizes, four disk masses, two variations on disk flaring, and ten inclinations. Table 1 shows the parameters used in our model grid. All models adopt an identical photon flux from the central star and assume fully evacuated cavities containing no material. These model images were convolved with the HST WFC3IR point-spread function (PSF) for the F160W filter. In this paper, we are primarily interested in variations in the observed near-IR morphology due to changes in envelope density, cavity half-opening angle, and inclination.

Zoom In Zoom Out Reset image size

In these models, the mass infall rate is used as a parameter to control the densities of the envelopes. The infall rate is combined with an adopted central stellar mass of 0.5M_⊙ to scale the envelope density using Equation (3) from Kenyon et al. (1993). See Furlan et al. (2016) for further discussion on this scaling.

The disk and envelope dust opacities are from a model by Ormel et al. (2011) that adopts a 2:1 mixture of ice-coated silicates and bare graphite grains, where the depth of the ice coating is 10% of the particle radius. The particles are subjected to time-dependent coagulation; we choose a coagulation time of 0.3 Myr. This is identical to the dust model used in Fischer et al. (2014) and Furlan et al. (2016). In the near-infrared, the opacities predicted by this model are slightly smaller than those of the often-cited OH5 opacities (Ossenkopf & Henning 1994). The reasons for adopting this model are described in Furlan et al. (2016). In Appendix E, we assess the dependence of the cavity appearance on the dust law.

Motivated by the shape of the observed outflow cavities, we used a parabolic model (n = 2) shown in Figure 1 for the outflow/envelope boundary in our models. In Section 4.2, we relax this constraint and use the power-law fit

$$\begin{eqnarray}&&z=A{\left|\rho \right|}^{n},\end{eqnarray} \tag{ 1 }$$

where the resulting power-law index, n, may be 1 or greater. For a given power law, the cavity half-opening angle depends on the adopted outer radius of the envelope; only for the case of a conical cavity (n = 1) is the half-opening angle independent of the adopted outer radius.

Examples of models from the grid are shown in Figure 2, which displays the effect of differing inclinations and cavity half-opening angles. Several model parameters, such as the radius and temperature of the central protostar or the presence of hot spots, are not constrained by either the SEDs or the near-infrared images. The surface brightness found in an image depends on the monochromatic luminosity of the protostar (which in turn depends on the temperature, radius and presence of hot spots), but the morphology of the image depends primarily on envelope density, outflow cavity shape, and inclination. The rest of our model parameters are chosen to cover a range of physical parameters observed in the fitting done by Furlan et al. (2016). This allows us to compare (in Appendix F) the values for the parameters determined by the fits to the SEDs and those determined from the near-infrared images.

Zoom In Zoom Out Reset image size

As shown by the models, the observed morphologies of the cavities trace the light scattered at a discontinuity in the dust density; in this case, the discontinuity is the boundary of a cleared cavity. If the protostar is seen edge-on, both cavities carved by the bipolar outflow are apparent. For these edge-on cases, a dust disk obscures the scattered light, creating a dust lane (Figure 2). If the system is inclined such that the extinction toward the far cavity is significantly higher than that toward the nearer one, a bowl-shaped unipolar structure is seen due to the obscuration of the more distant cavity. The envelope itself can be directly illuminated if the density is low enough for near-infrared photons to penetrate past the cavity walls and scatter off grains deep in the envelope. In these cases, the disk can cast shadows in the envelope that are also apparent for edge-on inclinations.

To compare our cavities to those measured in other analyses that adopt different models for their shapes (e.g., Arce & Sargent 2006; Furlan et al. 2016), we will determine the fraction of the envelope volume within the cavities. This is a measure of the amount of gas cleared by the outflow. The volume of the cavities in these models depends only on the power-law exponent n, the half-opening angle θ, and the outer envelope radius ${R}_{\max }$ (Figure 1). In Figure 3, we show the dependence of the fraction of the envelope volume cleared by the cavity on the cavity half-opening angle and the cavity exponent.

Zoom In Zoom Out Reset image size

An alternative metric for characterizing cavity sizes is the fraction of the envelope mass cleared by the outflows, i.e., the fraction of mass that would be found in an initially spherical symmetric core with a ρ^−1.5 density law and an outer radius of 8000 au. We compare the volume and mass fraction cleared in Figure 4. We find the fraction of the mass cleared can be up to 9% more than the volume cleared, and that the volume cleared is a lower limit to the mass cleared. We note that this is an instantaneous mass fraction of the current envelope, and it will differ from the total fraction of the envelope mass entrained and ejected by the outflow over the history of a protostellar collapse. Furthermore, it does not include the mass launched and ejected from the system by disk winds, X-winds, or accretion-driven stellar winds (e.g., Watson et al. 2016).

Zoom In Zoom Out Reset image size

In this section, we classify the protostars on the basis of their 1.60 μm morphologies in the HST images. We then examine how the morphologies depend on the properties derived from the model fits to their SEDs. For protostars with detected outflow cavities, we develop an algorithm to measure the shape of the outflow cavity, and we calibrate this approach using the radiative transfer models in our grid.

We exclude the images taken with the F205W filter from this analysis, as only the 83 objects successfully observed with NICMOS have these data. Furthermore, with a small number of exceptions, the morphologies are identical in the two NICMOS bands.

4.1. Protostellar Morphologies

The HST images resolve protostars at various stages of evolution, different inclinations, and differing amounts of envelope material. In these images, light primarily from the photospheres of the central protostars is scattered by dust grains in the envelopes, delineating structures present in the envelopes. In many of the images, the structures are similar to those caused by the outflow cavities in our model grid.

As a first step in our analysis, we divide all protostars into five morphological categories (Figure 5). The presence of a bipolar nebula, such as two scattered-light lobes separated by a dark lane or two outflow cavities, define the bipolar category. Sources with only one cavity visible make up the unipolar category. Unresolved protostars without detectable nebulosity are defined as point sources. Sources too deeply embedded to detect in the F160W band are considered nondetections (not shown in Figure 5). The final category comprises irregular protostars; these may result from background contamination (e.g., coincidence with a more extended reflection nebula) or from true inhomogeneities in the structure of the protostellar envelope. For the bipolar, unipolar, and irregular categories, the presence of an unresolved point source in the nebula is noted; these are likely to be the central protostar or light scattering off of structures within 80 au of the protostar, which is the smallest scale we can resolve in our images.

Zoom In Zoom Out Reset image size

In total, 141 HOPS objects exhibit extended structures in scattered light. The classification of all protostars are found in Table 4, and their breakdown is summarized in Table 2. Of these, 31 show a bipolar structure indicative of an edge-on inclination, although some cases show the point source of the central protostar near or offset from the midplane of the dark lane, implying that they are not exactly edge-on. One bipolar source was serendipitously observed in the same field as HOPS 334. This source was first identified as a candidate protostar by Stutz et al. (2013); based on their values for T_bol and L_bol, it is determined by the criteria in Furlan et al. (2016) to be a Class 0 protostar. In this paper, we introduce this source into the HOPS catalog as HOPS 410 (Table 4). Fifty-nine objects show nebulosity appearing to be a cavity on one side, with 36 of those having detected point sources near the base of the cavity. Fifty-one remaining protostars are classified as irregular. Images of sources with unipolar, bipolar, irregular, and point-like morphology are shown in Appendix A. Two additional Class II sources with nebulosity that were serendipitously discovered in our observation are shown in Appendix C.

Table 2. Breakdown of F160W Morphologies

	Point	No Point	Total
	Source	Source
Nondetections	⋯	60 + 3	60 + 3 a
Point Source b	93 + 7	⋯	93 + 7 c
Irregular	39	12	51
Unipolar	36	23	59
Bipolar	16	15	31
Total	184 + 7	110 + 3	294 + 10 d

Notes.

^a One of these three sources is likely an extragalactic source, and two are of uncertain nature (Furlan et al. 2016). ^b Sources without associated nebulosity ^c Six of these seven sources are likely extragalactic sources, and one is of uncertain nature. ^d Includes seven likely extragalactic sources and three of uncertain nature.

Download table as:  ASCIITypeset image

Approximately half of our sample, 163 objects, have no resolvable nebulosity in these observations. ¹¹ One hundred of these are detected as one or more isolated point sources; these have been analyzed to determine the companion fractions throughout the Orion molecular clouds (Kounkel et al. 2016). We refer to these as point sources without associated nebulosity. In these cases, any nebulosity around the source appears to be part of an extended nebula that is illuminated by other stars in the region or is very faint and tenuous and does not delineate a clear structure around the point source. As we will discuss in Section 4.4, the scattered light from cavities and envelopes illuminated by these sources are likely too faint to detect against the PSF of the point source. The remainder of the sources are nondetections.

Emission along jets, most likely dominated by the [Fe ii] line at 1.66 μm, is observed in 13 protostars, with three additional tentative detections. These are the bipolar protostars HOPS 133, 150, 186, and 216; the unipolar sources HOPS 29 (shown in Figure 5), 164, and 310; the irregular protostars HOPS 98, 188, 234, and 386; the point source 279; and the protostar HOPS 152, which, although not detected directly at 1.60 μm, is situated at a location that is an apparent source of jet emission. Tentative detections of jets are found toward the point-source protostars HOPS 3, 344, and 345.

In Figure 6, we plot the number of protostars versus bolometric temperature for four morphological groups: nondetections, point sources, unipolar or bipolar sources, and protostars with irregular morphologies. The bolometric temperature is a measure of the evolutionary stage of the protostar, although it also has some dependence on inclination (Ladd et al. 1998; Furlan et al. 2016). We also include the standard evolutionary classes, as determined with the criteria in Furlan et al. (2016). These figures demonstrate the strong dependence of detectability and morphology, in the near-IR, on the class of a protostar. The least evolved protostars (Class 0) are predominantly not detected, due to the greater optical depths in their envelopes. In comparison, the most evolved sources (i.e., flat-spectrum protostars and Class II pre-main-sequence stars) are dominated by unresolved point sources due to the low density of dust (and therefore low scattering probability) in their sparsely filled or nonexistent envelopes. ¹² Protostars with unipolar and bipolar cavities show a broad range of T_bol, but peak in the Class I phase (T_bol ∼ 100 K) and contain a significant fraction of Class 0 objects. Finally, the irregular protostars consist largely of Class I and flat-spectrum sources.

Zoom In Zoom Out Reset image size

We show the distribution of bolometric luminosities for each morphological class in Figure 7. The luminosity distributions for the three non-irregular classes display a shift in median luminosity, with point sources, unipolar/bipolar protostars, and nondetections having median L_bol of 0.5, 1.4, and 2.7 L_⊙, respectively. This change is small compared to the full range of bolometric luminosities probed, from 0.05 to 480 L_⊙. It is likely due to a decline in the luminosity with increasing age, as found by Fischer et al. (2017).

Zoom In Zoom Out Reset image size

4.2. Direct Measurements of Cavity Sizes

For protostars with unipolar or bipolar morphologies, we fit a power law to the shape of the cavities, in order to estimate the amount of the envelope that has been cleared by the outflows. This analysis relies on a custom edge detection routine developed to locate the outer contours of the cavities in the images. The methodology is illustrated in Figure 8. It is similar to the Sobel filter described in Danielsson & Seger (1990), constrained to the dimension perpendicular to the axis of the cavity. The image is first rotated such that the cavity is aligned with the positive y axis in an x − y Cartesian plane; this defines our adopted axis for the cavity. In three bipolar cases, those of HOPS 136, 280, and 333, we were able to measure the shape of both cavities. For each image, a 1D Gaussian smoothing kernel was chosen by eye to account for noise and applied to every slice of constant y. The width of the smoothing kernel is between 2 and 4 pixels, approximately 03.

Zoom In Zoom Out Reset image size

We then calculate the second-order finite difference along the slice (i.e., parallel to the x-axis) using the equation

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}I}{{{dx}}^{2}}=I(x+2)-2I(x+1)+I(x),\end{eqnarray} \tag{ 2 }$$

where I(x) is the F160W intensity at pixel x. The second-order finite difference, as an approximation to the second derivative, is zero at the inflection points of the slice. The inflection points allow us to define an "edge" of the cavity. The width of the smoothing kernel is increased to obtain a consistent edge, as a small smoothing kernel can produce a discontinuous edge. Inflection points are inspected to ensure that only those tracing the cavity (as opposed to structure within or outside the outflow cavity) are retained. We use this definition of an edge as it is bounded by the peak of the intensity and the background. More sophisticated techniques (e.g., Canny 1986) have a limited advantage, due to the presence of unrelated structures in the line of sight that cannot be treated as random noise.

To determine the half-width of the cavity, x, at a given position along the cavity axis, y, we measure the full width of the cavity between the two walls and then divide by two. Thus, the central axis of the outflow is defined by a curve tracing the midpoint of the two walls. Note that the y position is the distance along a straight line that starts at the base of cavity and extends along the adopted cavity axis (Figure 8).

In order to relate the detected edge to the physical cavity in the envelope, we ran the edge detection routine on our model grid. We compared the edges measured for the models as a function of the observed inclination to the shape of the projected cavity wall for the same model as observed from an edge-on inclination; this allows us to correct for the effect of inclination of the shape of the outflow. The location of the projected wall is given by the analytic equation y = A∣x∣ⁿ (Figure 1) , where A and n are determined by the parameters of our model that are described in more detail below.

For most models, the detected edges of the cavities differ systematically from those of the projected wall (Figure 8); this is due to the combined effects of inclination, the penetration of the light from the central protostars past the cavity wall into the envelope, and systematic biases of the edge detection routine. The inclination alone will broaden the cavity by 7%–25% for a 60°–40° (to correspond to 7°–25°, respectively) inclination, assuming a parabolic cavity.

Figure 8 shows our edge fitting routine applied to a model image of a protostar with an inclination of 60° and a cavity half-opening angle of 15°. The black solid line indicate the projected cavity wall of this model, as observed from this inclination. The black dashed line indicates where the cavity wall would be for the same analytical shape, but observed at an edge-on inclination—almost negligible, even for a 60° inclination. The detected edges (in blue) are characteristically wider than the known cavity wall.

To quantify the difference between the observed and actual edge, we determined the ratio of the distance to the "detected" edge in the model to that of the known, projected distance to the wall. For a given model, this ratio was found to be approximately constant as a function of the distance along the outflow axis. Thus, a single ratio can describe the difference between the observed and actual outflow cavity for a given source. Using the grid of models described in Section 3, the ratio was measured as a function of the cavity half-opening angle and inclination.

At lower inclinations, the line of sight toward the central protostar is more likely to be directly into the cavity or to pass through less envelope material, thus significantly increasing the probability to observe the protostar as a point source (see Section 4.4). In these cases, the cavity walls are difficult to detect against the PSF. Additionally, because there exist more possible lines of sight toward edge-on or nearly edge-on orientations than pole-on or nearly pole-on orientations, the probability that a protostar will be observed at a given inclination decreases as inclination decreases, assuming that the cavity may face any direction randomly. For these reasons, we averaged the ratios determined from all models for each cavity size, considering only inclination angles from 90° to 50°. The ratios are displayed as dashed lines in Figure 9, which shows that they are predominantly constant as a function of half-opening angle, except at the smallest opening angles, and that they have a weak dependence on inclination. The standard deviation over all parameters aside from cavity size and inclination are shown as the error bars in this figure.

Zoom In Zoom Out Reset image size

These ratios shown in Figure 9 were applied to the measured half-width of the cavities from the HST images. Generally, we initially divided the half-width by 1.3, which is the approximate average ratio for 50°–80° inclination cavities with a half-opening angle greater than 15°. For cavities of these sizes that were also bipolar and thus presumably near 90° in inclination, we restricted our initial ratios to 1.1. For cavities with narrower opening angles and unipolar and bipolar morphologies, we chose initial ratios of 2.0 and 2.3, respectively. We then iterated when necessary, modifying the correction ratio until the combination of the ratio and recovered half-opening angle was consistent with combinations observed in our models, shown in Figure 9.

We fit the cavity width as a function of the distance along the outflow axis to the function

$$\begin{eqnarray}&&y=A| x-{x}_{0}{| }^{n}+{y}_{0}.\end{eqnarray} \tag{ 3 }$$

to both the model grid discussed in Section 3 and observed images (Figure 10). Here, (x₀, y₀) identify the location of the protostar, and are fixed to the center of our model images. In the observed data, the parameters (x₀, y₀) are manually centered on the central protostar when apparent from a point source or an area of maximum flux along the profile of the cavity, or were placed along the disk absorption lane in the case of some bipolar sources. The midpoints of the cavity, as shown in Figure 11, were used to fit a center line, which in turn allowed us to perform a final small-angle rotation correction. Our two detected edges were then considered for fitting in three ways: both the left and right edges were independently fit with a power-law profile, and after folding over the now-vertical center line, both detected edges were simultaneously fitted. This allowed us to counter minor asymmetries in detected cavity edges as well as outlying points biasing our fitting regime on a single edge. The recorded parameters were those given by the single-edge fit with an exponent n greater than 1, or in cases where both edges met this criterion, the parameter from the folded fit was recorded. The exponent n, which is referred to as the cavity exponent, gives a measure of the collimation of the outflow cavity and may be indicative of the physical mechanism behind the outflow creation. For example, Shu et al. (1991) show how a shell of molecular gas composed of the outflow and swept-up material has a shape dependent on the angular distribution of the outflow.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By allowing n to be an unconstrained free parameter in our fitting (with the caveat that it be greater than 1), we allow for conical cavities (n = 1) as well as parabolic cavities (n = 2). The amplitude A parameterizes the size of the cavity. For the model used by the HO-CHUNK code, this relates the radius of the envelope (${R}_{\max }$) and the cavity half-opening angle θ by:

$$\begin{eqnarray}&&A={R}_{\max }^{1-n}{\cot }^{n}\theta .\end{eqnarray} \tag{ 4 }$$

The value of θ is only dependent on A for conical cavities (n = 1), but for other values of n, the value of θ depends on our choice of ${R}_{\max }$, which we set to 8000 au. Error analyses for functions of the fitted values are discussed in Appendix B.

For three protostars with bipolar morphologies, we were able to measure parameters for both cavities. In all other bipolar cases, the cavity appearing brighter was fitted. From our monochromatic model grid, we can see that inclination is responsible for variations in brightness between the two cavities. We expect the closer cavity to have a stronger signal due to lesser extinction along the line of sight, although inhomogeneous envelopes may also be responsible for differences in cavity brightness.

An example of our fitting technique applied to the protostar HOPS 136 can be seen in Figure 11. Fischer et al. (2014) determined that this protostar is a late stage object with a 10(±2)° half-opening angle for a ${R}_{\max }=10,000\,\,{au}$ envelope. The detected edges of its northern cavity are compared with power-law curves as given by Equation (3) in Figure 10, revealing that the northern cavity of this protostar is best fit by a 87 half-opening angle, ¹³ in close agreement with Fischer et al. (2014).

For 30 of the 90 protostars in our sample with unipolar or bipolar morphologies, we use this technique for measuring the cavity shape, and tabulate the values of n in Table 4 along with the half-opening angle. The median uncertainty for n, as obtained from the least-squares fitting of Equation (3), is δ n ∼ 0.14. We find from Equation (C1) that uncertainties in half-opening angle measurements are on average δ θ < (0.3°). We also include in Table 4 the volume fraction of the envelope cleared by the cavity as in Section 3. We calculate this assuming a spherical envelope and a cavity volume given by the profile in Equation (3). For simplicity, we calculate the uncertainty in this measurement for the case of a conical cavity of the same half-opening angle. Uncertainties in this metric are ∼5% of the measured volume fraction (Appendix B). We note that these uncertainties describe the accuracy of our fits after carefully determining the axis of the outflow cavity and selecting regions for fitting over which our edge detection was well-behaved and not confused by errant nebulosity, background stars, or the PSF of the central protostars. These determinations and selections may introduce systematic uncertainties that are not accounted for in the formal uncertainties.

The remaining 60 unipolar and bipolar protostars are excluded due to the inability to accurately and fully trace the cavity with the HST images. For several sources, we see a morphology indicative of an edge-on or nearly edge-on disk but do not see evidence for a cavity (e.g., HOPS 65 and HOPS 200); these may be pre-main-sequence stars with disks. Other protostars have cavities that are too faint to reliably trace (e.g., HOPS 220 and HOPS 235), show only one edge of a cavity wall—due to either a nonuniform extinction or an irregularly shaped envelope (e.g., HOPS 18 and HOPS 310), or are coincident with nebulosity—making it impossible to disentangle the cavity from larger-scale structures (e.g., HOPS 387 and HOPS 384). Finally, some cavities exhibit morphologies inconsistent with a power-law cavity (e.g., HOPS 8 and HOPS 232). In general, the factors that prevented automated fitting appeared incidental and not obviously correlated with apparent cavity size. Future efforts will focus on expanding the range of brightness levels and morphologies analyzed, as well as improving our understanding of the nature of cavities with only one apparent wall.

Figure 12 shows the range of fitted exponents n and cavity half-opening angles measured in this work. We find the mean and median of the cavity exponents are 1.9 and 1.5, respectively, indicating that parabolic cavities are a reasonable model assumption. Cavity exponents vary significantly between the protostellar cavities; we show examples of this variation in Figure 13. Finally, we show our distribution of cavity half-opening angles and volume fractions in Figure 14.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3. Cavity Sizes versus SED-derived Properties

The SEDs of the protostars provide information on both their evolutionary phase as well as their total luminosity (e.g., Whitney et al. 2003). Correlations between the SED-derived properties of protostars with the cavity sizes provide a means to probe the evolution of cavities as well as, potentially, their dependence on the final mass of the protostar (Fischer et al. 2017). Figure 15 shows two ways of parameterizing the cavity size, half-opening angle, and volume fraction cleared, against an assortment of evolutionary indicators and luminosities derived from the 1.6 μm to 870 μm SED (Furlan et al. 2016). We calculate the half-opening angles with Equation (4) using the values of A and n and an envelope radius of 8000 au.

Zoom In Zoom Out Reset image size

We quantify the degree of correlation by finding the Spearman Rank Correlation Coefficient r, a measure of the monotonic correlation between two variables in our thirty measured protostars. A correlation coefficient of 1 or -1 implies a strictly monotonic correlation. The Spearman coefficients and the two-sided p-value for a hypothesis test are given in Table 3 for each of the diagnostic indicators and methods of parameterizing the cavity size. The hypothesis test uses a null hypothesis of no correlation; therefore, low p-values indicate evidence of a correlation and evolutionary trend. For the three bipolar sources with both cavities measured, the found parameters of the two cavities were averaged before computing Spearman Coefficients and p-values.

Table 3. Spearman Coefficients and p-values

Evolutionary Diagnostic	Versus Half-opening Angle		Versus Volume Fraction
	r	p-value	r	p-value
Tbol	0.24	0.20	0.23	0.21
${\dot{M}}_{\mathrm{infall}}$	0.16	0.41	0.12	0.52
Lbol	0.26	0.16	0.29	0.12

Download table as:  ASCIITypeset image

We do not find statistically significant correlation between cavity size and T_bol or mass infall rate, (which should be considered a proxy for envelope density, as discussed in Section 3). As shown in Figure 14, the sample of protostars is dominated by Class I sources; at 1.60 μm, many Class 0 protostars are not detected, while flat-spectrum sources are often point sources or have irregular nebulosity (see Figure 6). Hence, these results can be primarily interpreted as a lack of evidence for an evolution in cavity properties across the Class I phase. The wide scatter in cavity sizes does not appear to be the result of evolution, but must depend on other environmental or intrinsic factors.

A higher correlation coefficient is found between cavity size and luminosity, with more luminous objects tending to have larger cavities; however, the p-values show that we cannot rule out the null hypothesis.

4.4. The Prevalence of Point Sources

We detect cavities toward 90 (30%) members of our sample, while 100 (33%) are observed as point sources without nebulosity. Since protostars are surrounded by envelopes that scatter light, the substantial number of point sources without detected scattered-light nebulosity is surprising. In this section, we examine why the point-source morphology is common and test whether the number of point sources implies an observational bias in our cavity size distribution.

Protostars may be observed as point sources without detected cavities in two primary cases. First, the central protostar is observed along a line of sight directly into the cavity. In this case, the brightness of the PSF from the central protostar, which will not be attenuated by the envelope, can be significantly stronger than scattered light from the surrounding cavity walls, which may only contribute a diffuse scattering around the bright protostar (Figure 2). Even if the line of sight grazes the cavity wall, the bright PSF can dominate over the nebulosity. Second, a low-density envelope leads to a more diffuse, lower surface brightness cavity wall and a brighter point source; once again, the cavity may not be visible against the PSF. In both of these cases, the extended nebulosity often found in the Orion clouds can also hide the scattered light from the cavities.

The first case may lead to a bias against detecting large cavities. For envelopes with large cavities, the central protostar can be directly observed over a larger range of inclinations. Furthermore, since the walls of the cavity are further from the protostar, they will have systematically lower densities than narrower cavities, and they will intercept less flux from the central star. Consequently, the walls will be fainter and harder to detect for large cavities (Figure 2).

To determine the combinations of inclinations, cavity sizes, and envelope densities that lead to the point-source morphology, and to ascertain potential biases in our observed cavity size distribution, we use a Monte Carlo simulation that combines the model grid in Section 3, the envelope densities from the SED model fitting of Furlan et al. (2016), and several adopted cavity size distributions. The steps of the simulation are as follows.

We first determined for each model in our grid whether a cavity would be detected by the WFC3 observations. To determine whether a cavity is detectable, two criteria were applied to each model. Nondetections of cavities were noted when no distinct edge that delineates a cavity is found in the image using the technique of Section 4.2 or when the signal in a cut taken across the cavity 2000 au from the central protostar has a peak value below the typical rms of an image. At 2000 au, every protostar with a detected cavity shows nebulosity; if the signal from the nebulosity in the models is below the typical rms values in the WFC3 images, then it is unlikely that the cavity would be detected. The typical rms was obtained from 30 × 30 off-source patches in our HST F160W images chosen to avoid point sources or outflow cavities. These patches commonly included extended, diffuse nebulosity that is common in the Orion molecular clouds.

We then performed a Monte Carlo simulation to predict the number of point sources without cavities we would detect for different assumed cavity half-opening angle distributions. We sampled the models drawing randomly for four parameters: infall rate (i.e., envelope density), inclination, inner F160W flux, and cavity half-opening angle. The distribution of infall rates for our sample of protostars (including those in the "irregular" category) were the best-fit values in Furlan et al. (2016). Inclinations were drawn assuming the outflow axes were randomly oriented. The maximum disk radius, the mass of the disk, and the radial exponent in the disk density law were left as free parameters to be randomly drawn from those in Table 1. The brightness in the inner 0.2" region for each model was determined by scaling the image flux to correspond to 1.60 μm magnitudes randomly drawn from the distribution of F160W magnitudes in the tabulation of Kounkel et al. (2016). Finally, we sampled the cavity half-opening angles from several different distributions discussed below.

We plot the fraction of models resulting in point sources as a function of infall rate (i.e., envelope density), inclination, and cavity half-opening angle in Figure 16. Here, we subtract the cavity half-opening angle of a source from the inclination, in order to measure the angle of the line of sight with respect to that of the cavity. Where this value is below zero, the line of sight is directly into the cavity and does not pass through the envelope. We find in Figure 16 a strong preference for point-source morphologies in models observed at such inclinations and opening angles. When the inclination minus half-opening angle is positive and near zero, the line of sight toward the central protostar intersects the lower-density, outer regions of the envelope. In this case, the incidence of a point-source morphology increases with decreasing infall rate. Finally, if the infall rate is low, point-source morphologies can be detected at every inclination and cavity half-opening angle combination, although the incidence increases at lower inclinations. As expected, point-source morphologies arise when either the protostar is observed through its outflow cavity or when the envelope is thin. This is consistent with the point-source morphology being dominated by protostars with flat-spectrum SEDs (Figure 6); the flat SEDs are expected for protostars observed at low inclinations or with low envelope densities (Calvet et al. 1994; Furlan et al. 2016).

Zoom In Zoom Out Reset image size

Each iteration ¹⁴ of the Monte Carlo simulation returns the number of point sources without detected nebulosity. We compare this to the number of point sources in our data. Before this comparison, we removed from our sample those sources identified by their SEDs as possible extragalactic contaminants or of an uncertain nature, as well as those without complete SEDs (Furlan et al. 2016), except for sources where HST imaging has revealed a unipolar or bipolar morphology, confirming their protostellar nature. Seventeen sources observed with WFC3 and classified as either nondetections or point sources were removed based on these criteria. Finally, we chose only the point sources observed with WFC3, in order to account for differences in sensitivity. This reduced our sample down to 230 protostars, with 70 point sources.

In Figure 17, \we show normalized histograms of the number of point sources observed for various models of the cavity half-opening angle distributions. In red, we show the simulation results when the cavity sizes are randomly drawn from the values in Table 4. The observed number of point sources is marked with a vertical line. Realizations of 230 protostars with this simulation attain 70 point-source detections or less at rate of 1.02%. We note that our exclusion criteria, described above, reject 11 objects with point-source morphologies from our sample. These objects could not be determined morphologically to be extragalactic contaminants, however, and were removed due to their SEDs. We note, however, that protostars may have extragalactic-like SEDs. HOPS 48, 67, and 301 were classified by Furlan et al. (2016) as extragalactic contaminants based on potential emission features in their Spitzer IRS spectra. In the case of these three sources, however, the features appear to originate in contamination from reflection nebulae or HII regions, and we observe cavities clearly associated with all three with HST WFC3. Thus, we consider our observed number of 70 point sources to be a lower limit.

Zoom In Zoom Out Reset image size

Table 4.  Index of All HOPS Sources in Our Sample

HOPS	R.A.	Decl.	SED	Tbol a	Instrument b	Point	F160W	Volume	Half-opening	Power-
ID	ICRS		Class a			Source c	Morphology	Cleared	Angle	law Fit
				(K)				(%)	(°)	n
29	05:34:49.0	−05:41:42.	I	148.2	WFC3	y	Unipolar	18.4	26.92	2.32
32	05:34:35.4	−05:39:59.	0	58.9	NIC	n	Unipolar	14.4	28.10	1.32
50	05:34:40.9	−05:31:44.	0	51.4	WFC3	y	Unipolar	15.4	27.97	1.48
64	05:35:26.9	−05:09:54.	I	29.7	WFC3	y	Unipolar	0.49	4.07	3.79
66	05:35:26.8	−05:09:24.	flat	264.9	WFC3	y	Unipolar	5.01	17.50	1.12
73	05:35:27.7	−05:07:03.	0	43.0	WFC3	n	Unipolar	0.529	5.28	1.41
81	05:35:27.9	−05:04:58.	0	40.1	WFC3	n	Unipolar	1.66	6.86	6.70
84	05:35:26.5	−05:03:55.	I	90.8	WFC3	y	Unipolar	28.3	35.89	1.76
93	05:35:15.0	−05:00:08.	I	107.3	NIC WFC3	n	Bipolar	5.10	17.3	1.18
135	05:38:45.3	−07:10:55.	I	130.3	NIC	n	Bipolar	22.1	38.41	1.02
136 d	05:38:46.5	−07:05:37.	I	161.7	NIC	n	Bipolar	1.77	8.71	2.53
150	05:38:07.5	−07:08:29.	flat	245.2	WFC3	y	Bipolar	11.6	26.21	1.18
164	05:37:00.4	−06:37:10.	0	50.0	WFC3	n	Unipolar	1.21	7.51	1.77
171	05:36:17.1	−06:38:01.	0	61.8	WFC3	n	Bipolar	8.19	21.34	1.29
185	05:36:36.9	−06:14:57.	I	96.9	WFC3	y	Unipolar	9.11	18.19	2.86
190	05:35:28.4	−06:27:01.	I	385.3	WFC3	y	Bipolar	13.3	26.22	1.45
219	05:41:29.2	−08:43:04.	I	90.0	WFC3	n	Bipolar	37.6	49.73	1.07
233	05:41:52.3	−08:01:21.	I	106.2	WFC3	y	Bipolar	0.71	4.88	3.87
244	05:41:01.9	−08:06:01.	I	127.3	WFC3	n	Unipolar	10.5	23.43	1.43
250	05:40:48.8	−08:06:57.	0	69.4	WFC3	y	Unipolar	29.8	38.39	1.54
273	05:40:20.8	−07:56:24.	I	243.3	WFC3	y	Unipolar	22.8	34.42	1.43
275	05:40:36.3	−07:49:06.	I	146.4	WFC3	y	Unipolar	22.5	35.72	1.26
280 d	05:40:14.9	−07:48:48.	I	121.2	WFC3	y	Bipolar	16.0	28.43	1.46
300	05:41:24.2	−02:16:06.	I	93.7	WFC3	y	Unipolar	1.04	6.78	1.97
322	05:46:46.4	+00:00:16.	I	71.3	WFC3	n	Unipolar	13.5	27.41	1.29
333 d	05:47:22.8	+00:20:58.	flat	240.9	WFC3	y	Bipolar	14.9	29.06	1.50
340	05:47:01.2	+00:26:21.	0	40.6	NIC	n	Unipolar	2.39	10.43	1.83
364	05:47:36.5	+00:20:06.	I	96.7	WFC3	n	Unipolar	5.21	16.90	1.32
367	05:54:36.2	+01:53:54.	I	249.4	WFC3	y	Bipolar	0.91	6.35	1.94

Notes. This abridged version only shows the 30 sources where additional measurements were obtained.

^a Furlan et al. (2016). ^b Objects observed with the instrument NICMOS are noted with the abbreviation "NIC." ^c Denotes (with y/n meaning yes or no) whether the object is observed with an apparent central point source. ^d Averaged parameters are reported for these bipolar sources where both cavities were measured.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

We also perform comparisons to fiducial models assuming a uniform distributions of cavities from 0 to 25, 30, 35, 40, and 45 degrees. The distributions extending beyond 35° include enough large cavities to overpredict the number of point sources. These results indicate that our observations are not significantly biased against the detection of large cavity openings.

Finally, we examined the consequences of outflow cavities that grow with time. We first adopt the relationship between cavity half-opening angle and T_bol found by Arce & Sargent (2006). We used this relationship and the observed distribution of bolometric temperatures of our protostars to derive the half-opening angle distribution we entitle the "Arce Model." We used a linear fit between the infall rate and T_bol to pick a model in our grid on each iteration of the Monte Carlo. This model overpredicts the number of point sources, as it does not take into account the highly evolved protostars with low cavity half-opening angles found in our sample (e.g., Fischer et al. 2014).

In summary, we find that the histogram produced from the observed distribution overlaps with the observed number of sources, although with a 1% probability of predicting the predicted number of protostars or less. If some of our excluded contaminant sources are in fact protostellar in nature, the observed distribution may provide a better match. It is important to emphasize that the result here is that our observed cavity angle distribution is largely consistent with our observed number of point sources. Uniform distributions of half-opening angles extending to 45° overpredict the number of point sources, and we do not find evidence that our observations fail to detect larger cavities. We also find that the uniform distributions with angles <35° reproduce the observed point sources better than our observed distribution does. This suggests that we may be missing small cavities that can be hard to detect, due to higher extinction from their envelopes.

The goal of this study is to assess the impact of jets and winds on protostellar envelopes. This is an essential step toward both understanding how feedback lowers the efficiency of star formation and determining the importance of feedback in halting mass infall and setting the final masses of protostars. Feedback can lower efficiency and mass infall in three ways: by ejecting mass that would have otherwise been accreted, by clearing the envelope, and by entraining gas in the envelope and surrounding cloud into an outflow (e.g., Watson et al. 2016; Zhang et al. 2016). This paper aims to quantify the role of cavity clearing.

One way outflows may halt infall is via the progressive clearing of the envelope as the protostar evolves (Arce et al. 2013). This may be driven, for example, by successive bursts of a wide angle wind (Zhang et al. 2019). The signature of this clearing would be a correlation between cavity size and the evolution of the protostellar SEDs. We find no significant trend between cavity half-opening angle with either T_bol or model inferred $\dot{M}$, both indicators of envelope evolution (Section 4.3). Instead, we find that there is a range of cavity half-opening angles extending from 5°–50° present across the observed range of T_bol and the range of $\dot{M}$ values inferred from model fits (Furlan et al. 2016). This implies that the evolution from dense to thin envelopes is not driven by the progressive growth of the outflow cavities.

To extend this result, we compare our cavity sizes with volume fractions calculated from millimeter and lower-resolution IR studies in Figure 18. We use the tabulated outflow cavity angles and assumed conical cavity shapes to calculate the volume fractions. Our scattered-light measurements extend these by providing a relatively large sample at a common distance observed with a uniform spatial resolution, which eliminates possible biases due to distance, and by detecting a significant number of protostars with relatively high T_bol ( > 100 K) and smaller cavities (<20% of the envelope cleared). The range of volume fractions (and hence, cavity half-opening angles) tabulated in the literature are consistent with those measured from our data, and there is no evidence for large systematic differences between the data sets, despite the different types of observations and methods used to measure the cavity sizes.

Zoom In Zoom Out Reset image size

Arce & Sargent (2006) use millimeter line emission in the blue and red lobes identified in CO maps to measure the cavity angle, assuming a conical outflow geometry. Although we do not share sources (so a direct comparison between the different methods cannot be made), we find that both the size scales probed and the range of observed volume fractions are similar, indicating that there are not large, systematic differences between the two techniques. Arce & Sargent (2006) suggest a correlation between an age diagnostic based on T_bol and the cavity size. This correlation, however, is driven significantly by the Class II objects in their sample, shown in Figure 18. For instance, the Spearman Rank Correlation Coefficient of T_bol and volume fraction decreases in their sample from 0.7 to 0.6 (p = 0.0011 to p = 0.015) without the Class II objects. By definition, these objects only contain small remnants of their protostellar envelopes. They lack well-defined defined cavities, and the lack of envelopes may not necessarily be the result of clearing by progressively expanding outflows cavities.

Seale & Looney (2008) also measure the opening angles of envelope cavities in scattered light detected by Spitzer IRAC. Although this technique has a lower angular resolution than our study and encompasses a sample of objects spanning a much broader range of distances, it has the advantage of being able to detect outflow cavities from Class 0 objects that are apparent in the Spitzer 3.6 μm band. They also find correlation coefficients that indicate no or weak statistically significant correlation between cavity size and age indicators, with two exceptions. The first age indicator that is correlated with cavity size is the IRAC color (Figures 11(b)–(d) of that work), an indicator that may also depend on cavity size, since larger cavities allow more radiation to escape at the wavelengths probed by IRAC. The second correlation is with the age parameter from the Robitaille et al. (2007) model grid. This age is used to set the sampling of cavity angles assuming cavity growth and thus could have induced a correlation. Furthermore, the correlation with the age parameter is relatively weak, with the Pearson product moment at a significance of α = 4% and no evidence for a correlation using Kendall's Tau rank correlation coefficient. Sources from Seale & Looney (2008) with bolometric temperatures in the literature are plotted as triangles in Figure 18.

Other works have found evidence for cavity growth during the Class 0 phase, as suggested by Arce & Sargent (2006). Velusamy et al. (2014) measure the full opening angle near the base of the cavity using the HiRes reduction of Spitzer IRAC images. They find a broken power-law growth showing a clear increase in the sizes of cavities with increasing T_bol from protostars with T_bol < 100 K, but do not reproduce the growth for more evolved objects. ¹⁵ Hsieh et al. (2017) also present a survey of low-luminosity protostars using IRAC images, in addition to CFHT WIRCam Ks-band observations. These authors use the same radiative transfer modeling code described in Section 3, but they use a direct least-squares fit of their model grid to their images to determine the cavity parameters. They find evidence for a similar growth during the Class 0 phase. Although we do not find a similar correlation in our data, the Class 0 phase is dominated by nondetection in our 1.60 μm imaging and the smallest cavities will be harder to detect (Figure 6). Thus, we do not rule out the growth of cavities during the Class 0 phase.

Furlan et al. (2016), using the SEDs and modeling described in Section 4.3, find that the envelopes decrease in density by a factor of 50 as protostars transition from the Class 0 to the flat-spectrum phase. ¹⁶ By the end of the Class I phase, it is thought that most of the stellar mass has been accreted. We should therefore expect the processes that reduce the mass and density of the protostellar envelope to continue through the Class I phase after starting in the Class 0 phase.

The lack of a correlation between the fraction of the volume cleared and the evolutionary indicators, in a sample preferentially probing Class I objects, implies that the evolution of the envelope during the Class I phase is not driven by growth of the outflow cavities. Although envelope clearing contributes up to a 40% reduction, the fact that the envelope density drops by more than an order of magnitude cannot be explained by this clearing alone. Of particular importance are the number of protostars with <15% of the envelope cleared throughout the entire range of T_bol covered. One of the best examples is the protostar HOPS 136, which has a volume cleared of 1.8%. Fischer et al. (2014) found that this protostar was in the late stages of stellar formation. The envelope mass of 0.06 M_⊙ was much smaller than the estimated stellar mass of 0.4–0.5 M_⊙, showing that most of the stellar mass has been accreted. A relatively low-density envelope is inferred from both the SED and the detection of scattered light in the envelope in the HST images (which implies a low optical depth at 1.60 μm). The presence of such protostars with low-density, low-mass envelopes and narrow outflow cavities at the late stages of stellar formation are clear examples of where the clearing of the envelopes by outflows cannot explain the observed low envelope densities.

Our results also put limits on the ability of feedback from outflows to explain the low star formation efficiency. Comparisons of the core mass function and initial mass function suggest that 60%–70% of the core mass will not accrete onto stars (Alves et al. 2007; Könyves et al. 2015), and previous authors have invoked outflows as partially responsible for this effect (e.g., Alves et al. 2007). Assuming the growth of the cavities is monotonic in time, the volume fraction cleared provides a lower limit on the mass fraction cleared by the outflows. From our HST data, the Class I protostars have cleared at most 40% of their volume (Figure 18). Recalling that the mass fraction cleared from a cavity may be as much as 9% higher than the volume cleared (Figure 4), the maximum fraction of mass cleared is 50%. Most of the protostars have cleared a much smaller mass fraction, even those toward the end of their protostellar phase; the median volume fraction cleared for the HST sample is only 10%. These results suggest that the feedback via clearing is not sufficient to explain the small star formation efficiency inferred for dense cores, and other mechanisms should be investigated.

There are other possible ways outflows may reduce star formation efficiency. The gas launched by the star-disk system in a jet or wind can escape the protostar and its envelope. Using estimates of the mass loss rates of 84 protostars, Watson et al. (2016) found that the median fraction of gas launched is 0.09 of the gas accreted (although with a wide dispersion); this may decrease the star formation efficiency by up to an additional 10%. We find a median star formation efficiency of ∼70% given a 10% median volume fraction cleared, a 9% increase for mass fraction cleared, and an additional 10% for mass directly launched and ejected by the central protostar. Only for the largest cavities, which clear up to ∼40% of their envelopes, can the efficiency be as low as ∼40%.

Second, the size of the cavity seen in scattered light may not measure the entire volume of the gas entrained in the outflow. In support of this, Seale & Looney (2008) noticed a possible discrepancy between their scattered-light outflow cavity sizes and the extent of the outflowing gas traced by millimeter line data. This outflowing gas may be slower-moving, denser gas entrained into the outflow that is located outside of the cavities.

In the case of the HH46/47 outflow, Zhang et al. (2016) used ALMA data in the ¹²CO, ¹³CO, and C¹⁸O lines to measure the mass in the outflow, including the slower, denser, entrained gas. They find that the gas mass in the outflow with velocities exceeding the escape velocity is ∼3 times the current stellar mass. If this instantaneous efficiency persists throughout the protostellar collapse, then the entrainment of gas in the outflow may account for the observed inefficiency. Simulations of collapsing cores with turbulence have also been able to achieve star formation efficiencies of 40% (Offner & Arce 2014). Radiative transfer models based on these simulations are needed in order to predict the evolution of cavities and compare them to the cavities measured in this work.

Finally, if outflows are not sufficient to reduce star formation efficiencies to the observed levels or to slow/halt accretion, then other mechanisms must be identified. For example, the collapse of a finite Bonner–Ebert core leads to an exponential tapering in the infall rate (Vorobyov 2010); however, this does not explain the inferred low star formation efficiency of cores. Furthermore, protostellar cores embedded in molecular clouds can draw gas from their surroundings and may not be limited by the mass in the surrounding core (Myers 2009). Oscillating molecular filaments, as suggested by Stutz & Gould (2016) and Stutz (2018) may eject protostars. Alternatively, dynamical interactions in small nonhierarchical systems or clusters may also eject protostars (Reipurth et al. 2010; Bate 2012). Identifying this mechanism should be considered a key problem in star formation, because it plays an important role in determining both the masses of stars and the efficiency of star formation.

We present WFC3 1.60 μm and NICMOS 1.60 μm and 2.05 μm images of 304 protostars and pre-main-sequence stars in the Orion molecular clouds. All of these objects were studied as part of the Herschel Orion Protostar Survey and are well-characterized by their 1.6 μm to 870 μm SEDs (Furlan et al. 2016). In this work, we use the images to resolve light from the central protostar scattered by dust in the envelopes surrounding the protostars, allowing us to probe structures with approximately 80 au spatial resolution. The specific results are as follows:

• 1.  
  We divide the sample into five distinct morphological classes. These morphological classes are nondetections (63), point sources without nebulosity (100), protostars with unipolar cavities (59), protostars with bipolar cavities (31), and irregular protostars (51). Thirteen of these protostars have jets appearing to originate from the protostars, and an additional three have tentative detections of jets. The relative incidence of each morphology depends on SED class: nondetections are dominated by Class 0 objects, protostars with cavities are dominated by Class I objects, and the point sources are primarily composed of flat-spectrum and Class I protostars. The irregular morphological class contains a relatively even mixture of Class 0, Class I, and flat-spectrum protostars. We find that nondetections have the highest bolometric luminosities, while point sources have the lowest.
• 2.  
  For the protostars with observed cavities, we developed an edge detection routine to find the structure of the cavity walls. From this, we fit a power law to the cavity shape and find the best-fit shape (e.g., conical, parabolic, etc.) for 30 protostars in our sample with unipolar or bipolar morphologies. We calibrated this technique against our large model grid to reliably measure the opening of cavities. We find a distribution of cavity half-opening angles ranging from 41–497, while the power-law exponent varies from 1.1 to 6.7 with a median of 1.5. We note that these cavity angles are not correlated with the SED derived angles of Furlan et al. (2016), demonstrating that fitting radiative transfer models to SEDs does not provide reliable constraints on cavity sizes (Appendix F).
• 3.  
  Using the well-characterized SEDs of Furlan et al. (2016), we look for correlations between the observed cavity half-opening angle and evolutionary diagnostics such as SED class and bolometric temperature. Our data show no evidence for a dependence of outflow half-opening angle and volume fraction cleared with any of the evolutionary indicators. Furthermore, several evolved protostars with relatively small cavity sizes are identified. We conclude that there is no systematic growth of the cavity half-opening angle during the Class I phase.
• 4.  
  We find that the incidence of point sources is consistent with both the observed cavity angle distribution and the distribution of envelope densities from Furlan et al. (2016). This implies that the point sources are protostars observed through a line of sight passing through the outflow cavity (hence seeing the protostar directly) or protostars with lower envelope density (as are typical of flat-spectrum protostars). Furthermore, we show that the number of point sources is inconsistent with a significant population of large cavities missed by our survey. Instead, our sensitivity to detecting cavities may decrease toward the smallest opening angles. As a whole, this is evidence that the cavity size distribution we obtain is reasonably complete and representative of the true distribution.
• 5.  
  Our findings indicate that outflow clearing is not the primary mechanism for the dissipation of the envelope during the Class I phase. They further suggest that clearing alone cannot explain the ∼30%–40% star formation efficiencies inferred from core mass functions. Current measurements of the amount of mass directly launched by protostars in winds or jets suggest that this additional factor is not sufficient. Measurements of the molecular gas with millimeter interferometry are needed to determine whether slower, higher-density flows entrained by the outflows are responsible for the halting of infall/accretion and the ∼30%–40% star formation efficiencies. If they are not, mechanisms other than feedback may be required.

This work was supported by NASA Origins of Solar Systems grant 13-OSS13-0094.

Based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program #11548.

Support for programs #11548, #14181, and #14695 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.

R.P. and S.T.M received funding from the NASA Astrophysical Data Analysis program (ADAP) 80NSSC18K1564.

A.S. gratefully acknowledges funding support through Fondecyt Regular (project code 1180350) and from the Chilean Centro de Excelencia en Astrofísica y Tecnologías Afines (CATA) BASAL grant AFB-170002.

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Software: Astropy (The Astropy Collaboration et al. 2013), DrizzlePac (STSCI Development Team 2012), HO-CHUNK (Whitney & Hartmann 1992, 1993), Matplotlib (Hunter 2007), NumPy (Walt et al. 2011), SciPy (Oliphant 2007), TinyTim (Krist et al. 2011).

The images in this appendix show the NICMOS and WFC3 images of protostars that display bipolar or unipolar morphologies, are point sources without associated nebulosity, or are classified as irregular (Figure 19). The complete figure set (262 images) is available in the online journal.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (262 panels)

Tarball of figure set images

Uncertainties for functions of the fitted values described in Section 4.2 were computed without assuming that the fitted parameters are fully independent, particularly n and A. Since θ is found as a function of these two parameters, the uncertainty δ θ is given by

$$\begin{eqnarray}\begin{array}{rcl}\delta \theta & \leqslant & \left|\displaystyle \frac{\partial \theta }{\partial A}\right|\delta A+\left|\displaystyle \frac{\partial \theta }{\partial n}\right|\delta n\\ & = & \displaystyle \frac{{R}_{\max }{\left({{AR}}_{\max }\right)}^{1/n}\left(A\left|\mathrm{ln}\left(\displaystyle \frac{{R}_{\max }}{A}\right)\right|\delta n+n\,\delta A\right)}{{{An}}^{2}\left({R}_{\max }^{2/n}+{A}^{2/n}{R}_{\max }^{2}\right)}.\end{array}\end{eqnarray} \tag{ B1 }$$

We note that this uncertainty is more strongly dependent on uncertainties in n than on those in A.

We calculate the adjusted uncertainty in the power-law coefficient A in Equation (1) by the equation:

$$\begin{eqnarray}&&\delta A=\tfrac{1}{{C}^{n}}\delta {A}_{\mathrm{fit}}+\tfrac{A\,n}{{C}^{n+1}}\delta C+\tfrac{A}{{C}^{n}}\mathrm{ln}(C)\delta n,\end{eqnarray} \tag{ B2 }$$

where C represents the correction factor shown in Figure 9 used to account for the effects of inclination on where cavity edges are detected, δ C is the uncertainty for a given correction value, δ A_fit is the nonadjusted uncertainty in parameter A resulting from least squares fitting Equation (1) to the location of detected cavity edges, and δ n is the uncertainty in n given the same fit.

By approximating the observed cavities as conical outflow cavities of half-opening angles θ, (as computed from A and n from the relation in Equation (4)), the estimated uncertainty in the fraction of the envelope volume subtended by our measured outflow cavities is given as

$$\begin{eqnarray}&&\delta {f}_{\mathrm{vol}}=| \tfrac{\partial {f}_{\mathrm{vol}}}{\partial \theta }| \delta \theta =\sin (\theta )\delta \theta ,\end{eqnarray} \tag{ B3 }$$

where f_vol is given by

$$\begin{eqnarray}&&{f}_{\mathrm{vol}}=1-\cos (\theta ).\end{eqnarray} \tag{ B4 }$$

Values for these uncertainties are discussed in Section 4.2.

Approximately 200 pre-main-sequence stars with disks, or Class II objects, were serendipitously in our WFC3 observations; these are tabulated in Kounkel et al. (2016). We have found that two of these objects are associated with bright, compact nebulosity similar to that found around protostars (Figure 20). MGM 2742 (V2475 Ori) is a binary that is associated with a nebula with an irregular morphology. MGM 925 (V2674 Ori) appears to be seen nearly edge-on and has a clear bipolar morphology.

Zoom In Zoom Out Reset image size

We are able to revise our identification of three objects by their morphology, showing them to be contamination. HOPS 339, shown in Figure 21, is determined to be a disk galaxy. Furlan et al. (2016) describes its SED as mostly flat with a strong 10 μm absorption feature and notes that by SED alone it would not be flagged as a possible extragalactic contaminant. This illustrates the importance of high-resolution near-infrared observations in disentangling galactic contaminants and protostellar objects. One protostellar candidate identified in Stutz et al. (2013), STS2013 038002 (Figure 21), was targeted by the HST program GO 14695. The source appears to be an extended streak; we suggest that it is a background galaxy observed through substantial extinction. The object STS2013 92011 (Stutz et al. 2013) was also targeted by GO 14965. It is shown by WFC3 imaging to be an outflow knot (Figure 22).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In our modeling, we used the dust opacity models adopted by Furlan et al. (2016). To investigate the role of the assumed dust law on the observed morphology, we compared images generated with two dust opacity models from Ormel et al. (2011). The opacity model used in this paper, "icsgra3," is described in Section 3 and adopts a grain coagulation time of 0.3 Myr. We compare this to the opacity model "icsgra2," which in contrasts adopts a time of 0.1 Myr, with consequently more grains of smaller sizes. When comparing model protostars from our grid generated with otherwise identical parameters, we observe that those using "icsgra2" show strongly limb-brightened cavity edge profiles and bright point sources (Figure 23). In contrast, the larger grains present in "icsgra3" result in more forward scattering where the intensity peaks toward the center of the cavity instead of the edges. The "icsgra3" are more consistent with HST observations, which typically show the cavities filled with emission (Figure 5), although there are some examples that show enhanced edges (Figure 13). (See Figure 19 for a complete set of images of all detected sources.) This suggests that the larger grains in "icsgra3" are more representative of our sample (Appendix A). Although beyond the scope of this investigation, future studies of the observed cavity morphologies may provide new constraints on grain properties and their variations.

Zoom In Zoom Out Reset image size

Most of the protostars in this paper have been characterized in detail by Furlan et al. (2016) using modeling in concert with the SEDs of these objects. The models used by these authors differ from those discussed in Section 3, and are fit to the SED from 1.6 μm to 870 μm. The primarily differences that affect the 1.60 μm emission are that our grid uses a finer sampling of high-inclination protostars, a sparser grid of envelope densities, and a different cavity opening angle exponent (2 versus 1.5).

Figure 24 shows histograms of the number of protostars versus inclinations determined by SED fitting. Four different histograms are displayed, one each for bipolar, unipolar, point sources, and nondetections. (Irregular protostars are not shown). Since the bins are chosen to have equal intervals in the cosine of the inclination, a random distribution should result in an equal number of sources in each bin; however, the SED-determined inclination for the overall sample of HOPS protostars peak at 60°–70°, suggesting that the SED-derived inclinations may have systematic biases (Furlan et al. 2016). The distribution of inclinations for point sources is similar in this respect to the overall sample (see Figure 29 of Furlan et al. 2016), except at the highest inclinations. Obscuration by the disk likely accounts for this deficiency, both by decreasing the ability to detect a point source and by decreasing the contrast with the surrounding nebulosity. The unipolar and nondetections also peak around 60°–70°, although they show a deficiency of low-inclination objects. This is expected, because the outflow cavities cannot be detected at low inclinations; the lower obscuration at these angles makes nondetections unlikely.

Zoom In Zoom Out Reset image size

Finally, the bipolar protostars show a broad range of inclinations, even though their observed morphologies require a nearly edge-on perspective. Furlan et al. (2016) showed that there are large systematic uncertainties in the SED-derived inclinations. This figure further demonstrates the limitations of using SED-derived inclinations, particularly for edge-on protostars. In a detailed study of the HST morphology of the HOPS 136 protostar, Fischer et al. (2014) could only find agreement between the SED models and the edge-on morphology by adding a low-density component of dust in the outflow cavity in order to increase the scattering at shorter wavelengths. This suggests that our models are incomplete and therefore underpredict the brightness of protostars in the near-IR; consequently, model fits erroneously favor inclined models where the near-IR emission is less absorbed by the disk.

Our measurements of the cavity half-opening angles can also be used to test the angles derived from SED models. Furlan et al. (2016) fitted models with discrete cavity half-opening angles of 5°, 15°, 25°, 35°, and 45° for cavities with a r^1.5 power-law shape. ¹⁷ Figure 46 of that work shows consistency of the best-fit cavity angle compared to the mode for various criteria of close models; these show there are a wide range of cavity angles that can be fit to a given SED. This is further demonstrated in Figure 25, where we compare half-opening angles inferred from the best SED fits from Furlan et al. (2016) to the cavity half-opening angles derived in this work. In spite of the different cavity shapes assumed, a monotonically increasing correlation should be seen if both are accurate measurements of the cavity size. Since we do not see the expected correlation, and since most protostars have decent-quality model fits for a range of cavity angles in the SED model grid, we find the cavity opening angles are not constrained well by SED modeling. A similar result was found by Seale & Looney (2008).

Zoom In Zoom Out Reset image size