{'fraction_non_alphanumeric': 0.03677881695615207, 'fraction_numerical': 0.03981537714077493, 'mean_word_length': 4.471291866028708, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Systematic Search for Molecular Outflows Toward Candidate Low-Luminosity Protostars and Very Low Luminosity Objects

We present a systematic single-dish search for molecular outflows toward a sample of 9 candidate low-luminosity protostars and 30 candidate Very Low Luminosity Objects (VeLLOs; L_int<0.1 L_sun). The sources are identified using data from the Spitzer Space Telescope catalogued by Dunham et al. toward nearby (D<400 pc) star forming regions. Each object was observed in 12CO and 13CO J = 2-1 simultaneously using the sideband separating ALMA Band-6 prototype receiver on the Heinrich Hertz Telescope at 30 arcsecond resolution. Using 5-point grid maps we identify five new potential outflow candidates and make on-the-fly maps of the regions surrounding sources in the dense cores B59, L1148, L1228, and L1165. Of these new outflow candidates, only the map of B59 shows a candidate blue outflow lobe associated with a source in our survey. We also present larger and more sensitive maps of the previously detected L673-7 and the L1251-A IRS4 outflows and analyze their properties in comparison to other outflows from VeLLOs. The accretion luminosities derived from the outflow properties of the VeLLOs with detected CO outflows are higher than the observed internal luminosity of the protostars, indicating that these sources likely had higher accretion rates in the past. The known L1251-A IRS3 outflow is detected but not remapped. We do not detect clear, unconfused signatures of red and blue molecular wings toward the other 31 sources in the survey indicating that large-scale, distinct outflows are rare toward this sample of candidate protostars. Several potential outflows are confused with kinematic structure in the surrounding core and cloud. Interferometric imaging is needed to disentangle large-scale molecular cloud kinematics from these potentially weak protostellar outflows.


Introduction
The incredible sensitivity of the Spitzer Space Telescope has resulted in the study of a new class of very low luminosity objects (VeLLOs). Formally, a VeLLO is defined as a protostar with an internal luminosity of < 0.1 L ⊙ . The initial Spitzer testing and validation observations for the Cores to Disk Legacy team (c2d) detected a new low-luminosity protostar in the previously classified starless core L1014 (Young et al. 2004). Follow-up observations revealed that L1014 was indeed a VeLLO with a weak molecular outflow (Bourke et al. 2005) and centimeter continuum emission consistent with a thermal jet (Shirley et al. 2007). Since the discovery of L1014, the c2d project has completed a census of the embedded protostellar population in the nearby (D < 400 pc) low-mass star-forming molecular clouds (see Evans et al. 2009 for a summary). The c2d survey is sensitive to bolometric luminosities of  Table 1 for the adopted distances to regions observed). The subset of embedded low luminosity protostars (L bol ≤ 1 L ⊙ ) including VeLLOs were catalogued by Dunham et al. (2008) and provide the most complete collection that may be systematically studied.
The true nature of VeLLOs are a subject of much current research. Some VeLLOs may be protostars in the incipient stages of accretion, the observationally elusive first hydrostatic core phase of their evolution (Larson 1969; see Machida et al. 2008 for a recent theoretical simulation). Alternatively, some VeLLOs may be more evolved protostars which happen to be observed in a low accretion state. Observations of the molecular outflow, a byproduct of magnetized accretion onto a protostar, provide a potential diagnostic to constrain these two possibilities. Molecular outflows may provide a fossil record of the accretion history of the source.
The source IRAM04191+1522 was one of the first discoveries of a low-luminosity source (L int = 0.08 L ⊙ ; Dunham et al. 2006) driving a large-scale molecular outflow (André et al. 1999, Belloche et al. 2002. Since the identification of the outflow around the VeLLO L1014, five more VeLLOs have had outflows detected: L1521F , MMS126 (Stanke et al. 2006), Per 073 (Hatchell & Dunham 2009), L673-7 (Dunham et al. 2010a), and L1148 (Kaufmann et al. 2011. Four out of the seven VeLLOs drive large-scale outflows that were detected with single-dish telescopes (IRAM04191, L673-7, Per 073, and MMS126), two of the sources (L1014 and L1148) have outflows that are detected only with interferometers (not detected with single-dish), and one source (L1521F) has a near-infrared scattered light nebula indicative of outflow activity but no direct detection of the outflow. In addition, there are now a few objects with interferometrically detected outflows, but no detected Spitzer infrared driving source: L1448 IRS2E1 (Chen et al. 2010), Per-Bolo 58 (Dunham et al. 2011), L1451-mm (Pineda et al. 2011), and CB17MMS (Chen et al. 2012). These objects are candidates for first hydrostatic cores although follow-up observations are still needed to confirm their status.
From maps of the emission in the molecular outflow line wings, the dynamical age, mass, momentum, and total outflow force may be estimated (see Cabrit & Bertout 1992). Resulting estimates of the average accretion luminosity are typically larger than the observed source luminosity implying that the accretion rate must have been larger in the past (e.g. Dunham et al. 2010a). This is another version of the famous "luminosity problem" for low-mass protostars, elucidated by Kenyon et al. (1990), where protostars are not luminous enough to account for the expected accretion luminosity based on estimates of the protostellar lifetime. A possible solution to this problem is that the accretion rate is variable with time. Modeling of source counts of protostars in different evolutionary stages from the c2d survey indicates that low luminosity protostars may accrete most of their mass over a period of a few 10 5 years through episodic events (see Dunham et al. 2010b). In order to test this hypothesis, comprehensive surveys for molecular outflows around low luminosity protostas and VeLLOs are needed.
The published searches for molecular outflows associated with VeLLOs have been heterogenous to date toward a very small fraction of the VeLLO population in nearby clouds and cores. In this paper, we present the first systematic survey for large-scale molecular outflows toward the sample of 39 isolated, candidate low-luminosity protostars (including 30 VeLLOs) in the J = 2 → 1 lines of 12 CO and 13 CO (see Table 1) . Our main goal is to characterize the frequency of kinematically distinct large-scale outflows toward VeLLOs and to determine the properties (mass, age, momentum, force, etc.) of any detected outflows. In §2 we describe the observing techniques and our sample selection. In §3 we analyze the properties of 3 detected outflows in the sample as well as describe the properties of the numerous non-detections. In §4 we discuss the nature of the VeLLO outflows.

Observations
Observations were performed on 9 days between December 2009 and March 2011 (see Table 2) using the sideband separating ALMA Band-6 prototype receiver on the 10-m Heinrich Hertz Telescope. This instrument permits simultaneous dual polarization observations of 12 CO J =2 → 1 and 13 CO J = 2 → 1 in the LSB and USB respectively with a 5 GHz IF. Observations were performed in 4-IF mode (dual polarization and sideband separating) using filter bank backends with 128 channels at 250 kHz resolution (0.34 km/s and 0.32 km/s at 220 GHz and 230 GHz respectively).
Pointing was carried out using the planets Jupiter, Mars, and Saturn with a typical rms pointing uncertainty of < 5 ′′ . Position switching observations with a 5 ′ OFF position in azimuth were used to calculate the main beam efficiency (T mb = T * A /η pol ; Mangum 1993). The calculated efficiencies for each observing session are given in Table 2. Sideband rejections were measured to be > 13 dB and were not corrected for. Variation in η pol likely reflect calibration variations and changes in the coupling between the receiver optics and the telescope between observing seasons.

Observing the Sample
Our target objects were selected from the catalogue of low luminosity protostars in the Cores to Disk Legacy sample compiled in Dunham et al. (2008) (their Tables 3 -8).
Throughout this paper, we use the region name followed by the source number in parenthesis listed in the tables of Dunham et al. (2008) (007)). The Dunham catalog contains 199 sources. All isolated sources (> 30 ′′ from a IRAS source) with suitable declinations observable from the Heinrich Hertz Telescope (δ > −30 • ) from Dunham Table 3 (Cores) and  Dunham Table 7 (Ophiuchus), which had not been observed previously, were included in our final target list of 39 sources. We excluded candidates within a beam width, 30 ′′ , of a known IRAS source, since it is unlikely that we will be able to disentangle the outflow from the low luminosity protostar from the stronger IRAS source. We also exclude some sources and regions where known single-dish searches for molecular outflows have already occurred (e.g. IRAM04191 André et al. 1999, L1521F Bourke et al. 2006, L328 Lee et al. 2009, L1014 Bourke et al. 2007, CB130 Kim et al. 2011, and LFAM 26 Bussmann et al. 2007, Nakamura et al. 2011. We exclude sources in Perseus (Dunham Table 4) and Serpens (Dunham Table  8) because those regions have been recently mapped at ≤ 30 ′′ resolution in transitions of 12 CO (Hatchell et al. 2007, Curtis et al. 2010Burleigh & Bieging in preparation). We did not observe the Lupus and Chamaeleon II regions as they are too far south and those regions are the the subject of a southern hemisphere companion survey (Dunham et al., in preparation). The final observed sample consists of 39 candidate protostars, listed in Table  1. Of this sample, 30 sources are classified as VeLLOs L IR < 0.1 L ⊙ and all sources have L IR < 0.45 L ⊙ .
For each source in Table 1, a 5-point map was obtained with 30 ′′ spacing and a minimum of 3 minutes absolute reference position switching. Clean OFF positions were found by manually checking for emission free regions using data from the Dame et al. (2001) survey included in SkyView CO maps 1 . Before using an OFF position we integrated for three to four minutes using a second reference OFF to insure that it was truly free of CO emission. A linear baseline was subtracted from all maps and all data were reduced using custom routines in the GILDAS (CLASS and GREG) software package 2 .

OTF Mapping
We mapped the sources L673-7 (031), L1251-A (045), B59 (022), B59 (023), L1148 (033), L1228 (036), and L1165 (040). The sources were chosen as either known outflow candidates where the outflow properties were uncertain and could be improved or as new outflow candidates based on the existence of line wings in the observed 5-point spectra ( §3.1). For each of these sources we created an on-the-fly (OTF) map with dimensions 5 ′ by 5 ′ , 4 ′ by 4 ′ , or 3 ′ by 3 ′ . For the L673-7 (031) outflow two regions of size 3 ′ by 3 ′ offset by ±120 ′′ in the directions of the red and blue wings were also mapped. All maps had 10 ′′ row spacing with a scan rate of 10 ′′ per second and were scanned boustrophedonically in the RA and DEC directions multiple times to improve the signal-to-noise in the line wings.
On-the-fly maps were reduced and analyzed using custom routines in the GILDAS software package. A linear baseline was removed from all observations before they were scaled using the calculated efficiency. These scaled V pol and H pol observations were combined and then convolved and re-gridded using a Gaussian-tapered Bessel function (Mangum et al. 2008). The resulting beam sizes where 35.5 ′′ for 12 CO and 36.9 ′′ for 13 CO. Integrated intensity maps were created using the convolved and re-gridded spectra. The red and blue wings were partitioned into multiple velocity channels with the minimum channel size being the spectral resolution (δV chan ). Using these files we generated channel maps for each source with the lowest contour at the two to three sigma level of integrated intensity (σ I = σ T δV chan (v 2 − v 1 ) where σ T was found by averaging the baseline rms at several points in the map). These channel maps allowed us to discern any large scale outflow structure for the outflow candidates.

Detections & Non-detections
The line properties for the center spectrum of each source are tabulated in Table 3. The ratio of integrated intensities for 12 CO to 13 CO ranges from 1.06 to 23.3 though the median is just 3.51. This ratio is much lower than the expected value of approximately 50 for optically thin emission (Wilson & Rood 1994); therefore we conclude that the 12 CO emission in these regions is typically very optically thick. This agrees with the general shape of line profiles seen in Figures 1 -10. Some sources display self-absorbed 12 CO line profiles with 13 CO peaking near the 12 CO absorption dip. The 13 CO line profiles are typically well fit by Gaussian line profiles with an average ∆v = 2.12 ± 0.98 km/s. This linewidth is much larger than the sound speed in 10 K molecular gas (0.077 km/s for CO) indicating that the motions are likely characterized by turbulence. Unfortunately, at the velocities of 12 CO line wings, 13 CO is typically not detected; therefore, we are unable to make a direct optical depth correction to the 12 CO line wing emission.
The 5-point maps for each source ( Fig. 1-10) were visually classified according to their spectral features (Table 3). Sources consisting of a single component were flagged with an "s" while sources consisting of multiple components were are denoted by "mc". Just over half of the sources have single component 12 CO spectra (53.3%) while one third of the 12 CO spectra appear to have multiple velocity components (33.3%). Self absorbed sources are distinguished from sources with multiple components by comparing to the 13 CO line: sources with 13 CO which peaks near the 12 CO absorption dip are considered to be self absorbed and are denoted by the flag "sa". The remaining 13.4% of sources show signatures of self-absorption in 12 CO.
Sources displaying evidence of blue or red wings in 12 CO are flagged with "bw" and "rw" respectively. Seven sources (23%) have either a prominent blue wing or a prominent red wing. These are the candidate outflow sources that we explore with OTF observations. The known outflow sources L673-7 (031), L1251 (044), and L1251 (045) are successfully recovered as outflow candidates in our visual flagging scheme. The visual search for outflow signatures is challenging because several sources have confused spectra with multiple velocity components. Nearby IRAS sources, despite being more than one beam (30 ′′ ) away, also contribute to the confusion in some spectra. The sources L673 (027), (028), and (030) are examples of complicated, multi-peaked spectra which appear to show evidence for an outflow, but all three sources are < 50 ′′ from IRAS19180+1116 whose outflow likely dominates the spectra shapes in this region. From the spectra in Figures 1 -10, we visually selected five additional candidates, B59 (022), B59 (023), L1148 (033), L1228 (036), and L1165 (040), to follow up with OTF mapping that show clear evidence of a blue or red wing.
Each candidate protostar in the Dunham catalog is assigned a group number (1 through 6) which indicates the likelihood that the object is an embedded protostar (see Table 1 for the group numbers). Sources in group 1 are verified protostars with the certainty decreasing with increasing group number. The nature of sources in group 6 is very uncertain due to a lack of supporting observations (i.e. no evidence that the source is embedded within a core, etc.). All three of the known outflow sources that we mapped belong to group 1 (L673-7 (031), L1251 (044), and L1251 (045)). There are only two group 1 or group 2 objects in our sample which were not selected as good outflow targets based solely on the appearance of their 12 CO spectra: L673 (027) and L1221 (041). Unfortunately both of these sources are nearby (42 ′′ ) strong IRAS outflow sources. As explained above, the kinematics of the L673 (027) region is severely complicate by the IRAS19180+1116 outflow source. L1221 (041) appears to have extended 4.6µm emission in the NW-SE direction that is indicative of an outflow, but unfortunately, the large outflow source IRAS22266+6845 to the northwest appears to overlap in orientation. The effect of the IRAS source outflow can be seen in the northernmost 5-point spectra for L1221 (041) although no outflow signature is apparent at the position of the VeLLO. This region has been observed with the BIMA interferometer (Lee et al. 2005) and a very tenuous detection of the VeLLO outflow is indicated. Deeper interferometric imaging is needed of L1221 (041) to confirm this result.
Of the remaining 5 candidate outflow sources that we mapped, all of them belong to group 6 except B59 (023) which is in group 3. Interestingly, B59 (023) is the only of those candidates that shows evidence of an outflow being associated directly with the source ( §3.4). Several of the Ophiuchus sources, (191) through (199) excluding (197), were not detected in 12 CO; all of these non-detections are in group 5. This result indicates that those protostellar candidates are not embedded objects and are not associated with molecular gas.
In summary, all of the confirmed outflow candidates reported in the literature and in this survey belong to groups 1, 2, or 3. In the following sub-sections we analyze in detail the properties of the detected and non-detected outflows observed in the OTF maps.

The Fully Mapped L673-7 (031) Outflow
Dunham et al. (2010a) discovered a large outflow surrounding L673-7 (031) (L673-7-IRS), but were not sure that the full extent of the outflow was mapped. Here, we present a new map which extends the region covered by Dunham and fully map the outflow associated with this VeLLO with higher signal-to-noise. We follow the methodology described in Dunham et al. (2010a) to derive the properties of the outflow. We explicitly show the functional dependence on the inclination in the calculations. Also, when appropriate, we always chose the lower limit in a calculation such that the final outflow properties represent lower limits to the actual values.
The outflow properties are first calculated per velocity channel per pixel, then summed over those pixels covered by the outflow and all outflow velocity channels in the intervals [−0.7, 5.3], and [9.3, 16.3] km/s (see Figure 11). Note that this velocity interval is more restrictive than the velocity interval used in Dunham et al. (2010a) and is determined by the velocity interval over which we statistically (> 2σ I ) detect the outflow signature in the map. For each calculated quantity (mass, momentum, etc.), we propagate the statistical errors on a channel-by-channel basis. The dominant source of statistical uncertainty is in the integrated intensity in the map. We do not propagate the distance uncertainty (a systematic error which itself is poorly constrained) into these calculations.
A lower limit of the dynamical time, t d , was found by dividing the extent of the outflow by the center velocity of the outflow wing relative to the v LSR of the protostar. We measure the extend of the outflow to be 290 ′′ compared to the 150 ′′ used by Dunham et al. (2010a). We calculate t d = 4.44 ± 0.23 × 10 4 cos i sin i yr assuming an angular uncertainty of half the beam size (15 ′′ ). Using this value, we find the time averaged mass accretion rate is Ṁ acc = 2.78 ± 0.14 × 10 −7 sin i cos 2 i M ⊙ yr −1 resulting in a total protostellar mass accreted of M acc =Ṁ acc t d = 1.23 ± 0.09 × 10 −2 1 cos i M ⊙ . The resulting lower limit on the time averaged accretion luminosity (L acc = GM accṀacc /R with R = 3R ⊙ ) is L acc = 0.036 ± 0.003 sin i cos 3 i L ⊙ . The ratio of the outflow-derived accretion luminosity to the internal luminosity is 2.1 ± 0.7 sin i cos 3 i where we have assumed a 30% uncertainty on L IR . While our calculated numbers are smaller than the values found in Dunham et al. (2010a), we are now confident that with this new map, the full extent of the L673-7 outflow has been mapped. This independent analysis still supports the hypothesis that L673-7 has undergone episodic accretion; however, the discrepancy between the derived average accretion luminosity and the observed internal luminosity has decreased significantly with our new analysis. It must be kept in mind, though, that we have selected lower limits to all derived outflow quantities except for the inclination. The unknown inclination of the outflow is the largest source of systematic uncertainty. Figure 12 shows that the functional dependence of the outflow dynamical time, force, and time-averaged accretion rate with inclination is monotonic for each quantity. For an inclination of i o = 34.3 • , the correction factor sin i cos 3 i = 1. The probability that the real inclination of the outflow is greater than an inclination i o is given by P (i > i o ) = cos(i o ) for randomly oriented outflows. The inclination for which sin i cos 3 i ≈ 1 corresponds to an 83% confidence interval that the real outflow inclination is larger than i o = 34.3 • . Thus, typical inclination angles will be larger than this value and the derived quantities can easily be factors of a few times larger to over an order of magnitude larger; therefore, any discrepancy between the time-averaged accretion luminosity and the internal luminosity will also increase by the same factor (see §4.1).
3.3. The L1251 (045) Outflow L1251 (045) was first mapped by Lee et al. (2010) (L IR = 0.077); however, only the properties of the larger, more distinct outflow associated with L1251 (044) were calculated (L IR = 0.156). L1251 (044) and (045) are labeled as IRS4 and IRS3 respectively in the Lee et al. (2010) paper. In order to determine the properties of this adjacent outflow we re-mapped the L1251 (045) outflow using higher spacial resolution (30 ′′ versus 48 ′′ in Lee et al.) and with higher signal-to-noise. We use the same method as for L673-7 to derive the outflow properties ( §3.2). The velocity intervals for the blue and red wings span from [−10, −5] km/s and [−1, 0.5] km/s. The smaller velocity difference for the red wing is used to avoid integrating over a second velocity component which peaks +3.8 km/s from the v LSR of the primary component (see Figure 13).
The calculated outflow mass, energy and momentum are slightly larger than for L673-7 despite this outflow being less distinct. Assuming a distance of 300 ± 50 pc (see Dunham et al. 2008) we calculate an outflow dynamical time of t d = 2.69 ± 0.23 × 10 4 cos i sin i yr. This outflow has been powered for just more than half as long (at the same inclination) as the L673-7 outflow; therefore, the derived outflow luminosity and force are higher for L1251-A (045) than for L673-7 (031). We calculate Ṁ acc = 6.3 ± 1.5 × 10 −7 sin i cos 2 i M ⊙ yr −1 . The lower limit on the protostellar mass accreted is 1.7 ± 0.6 × 10 −2 1 cos i M ⊙ , making the lower limit for the accretion luminosity L acc = 0.11 ± 0.04 sin i cos 3 i L ⊙ . The average (minimum) accretion luminosity is a factor of 1.5±0.7 sin i cos 3 i times higher than the observed internal luminosity. For inclinations greater than 34.3 • , the discrepancy will become larger ( Figure 12).

A New Potential Outflow Around B59 (023)
We have identified a new potential outflow around the source B59 (023). This source displays line wings from a blue lobe to the southwest direction from the source (Figure 14). This feature has also been confirmed in CO J = 3 → 2 observations with the ASTE 12m telescope (Dunham et al., in preparation). Unfortunately, the red wing is confused by a powerful outflow associated with 2MASS source J17112317-2724315 to the northeast (see Brooke et al. 2007). It is not clear whether the blue shifted emission actually emanates from B59 (023), or whether we have mapped its full extent.
Since we are only able to discern one wing of this purported outflow, we report the outflow properties found for the blue wing multiplied by two. B59 (023) is assumed to be at the same distance (D = 125 pc) adopted by the c2d team for the Pipe Nebula and Ophiuchus complex. The derived outflow luminosity and force are significantly higher than for L673-7 (031) and L1251 (045), in part due to the apparent small dynamical time. B59 (023) is one of the nine sources in our survey that is technically not a VeLLO (L IR = 0.343 L ⊙ ). The (minimum) average accretion luminosity is 1.5 ± 0.5 sin i cos 3 i times higher than the observed bolometric luminosity; however, we caution that the derived parameters may be overestimated if the outflow indeed has a longer dynamical time. Given that the region is too kinematically confusing for trustworthy analysis using only a single dish, we believe this source to be a prime candidate for follow-up interferometric observations.

OTF Map Non-detections
All of the remaining OTF maps were kinematically confused and we were unable to discern any large scale, distinct bipolar outflow structure originating from our sources. The outflow lobes are not obvious from analysis of the channel maps or individual spectra in the HHT maps (see Figures 15-18). Thus, the remaining four sources cannot be classified as outflow candidates based on their 12 CO channel maps. B59 (022) is west of (023) and its spectra are also severely confused by the outflow from IRAS17081-2721 (visible across the northern edge of the channel maps). There does appear to be a red lobe to southwest of B59 (022), but the corresponding blue lobe has strong red wing emission overlapping in the same position (see the top two rows of the channel map in Figure 15). The channel maps for L1148 (033) show no distinct structure in the CO wing emission ( Figure 16). Outflow sources can been seen in the maps of L1228 (036) and L1165 (040) but the emission is clearly offset and not associated with the candidate protostars in either case (Figures 17 and 18). We calculate the mass sensitivity of our survey using the same technique as for the detected outflows but over an empty portion of the map. Using the map L1148 (033) and assuming a conservative outflow extent which covers half the map, the 3σ limit for a mass detection is 1.38 × 10 −3 M ⊙ which is consistent with the uncertainties for the detected outflows. We display the maps of these non-detections for completeness.

Analysis & Caveats
Of the 39 sources observed in this survey, three have previously known large-scale outflows. We initially identified five additional protostars as outflow candidates; but of those five new candidates, only B59 (023) showed evidence of a newly detected candidate outflow although confusion with the strong nearby IRAS17081-2721 outflow makes estimation of its outflow properties difficult. The VeLLO single-dish outflow detection rate is only 2 out of 30 sources or 7% in this survey. Recently, Hatchell & Dunham (2009) and Curtis et al. (2010) separately mapped and analyzed the outflow properties of protostars in the Perseus molecular cloud. Outflow signatures from only one VeLLO was uniquely identified in their maps (Per 073 in Hatchell & Dunham; L IR = 0.018 L ⊙ ) although many sources with luminosities approaching 0.1 L ⊙ were detected in both surveys (e.g., Per 071 in Dunham Table 4 which is source 65 in Curtis et al. 2010; L IR = 0.105 L ⊙ ). If we include the other previous VeLLO single-dish detections from the literature (IRAM014191, Per 073, and MMS126), then the single-dish detection rate increases to 13% for 38 VeLLOs observed in cores. These detection rates are low and suggest that large scale outflows from VeLLOs in the Dunham catalog are rarely a dominant kinematic feature in the surrounding molecular cloud.
We compare the properties of the detected VeLLO outflows in this survey with those of other known VeLLOs and low-luminosity objects observed from the Dunham catalogue. The outflow properties are taken from the individual paper references for each source (IRAM01491 André et al. 1999, MMS 126 Stanke et al. 2006, L1014 Shirley et al. 2007 Figure 19 were reported directly in each paper; therefore, in those cases the quantities were derived using the same equations as in this paper. If a range of values were given, the average value is plotted. Errobars were rarely reported in the original references and are not shown in the figure. The plotted values in Figure 19 do not make optical depth corrections or inclination corrections. We used the literature values with the exception of MMS 126 from Stanke et al. (2006) whom made an average optical depth correction (factor of 3.5) and assumed an inclination of 57.3 • (we have reversed those corrections for comparison).
There is strong evidence that the outflow properties change as protostars evolve from the Class 0 to the Class I phase; namely, Class 0 sources drive faster, more massive molecular outflows than Class I sources (Arce & Sargent 2006). The outflow mass, momentum, and force decrease from the Class 0 to the Class I phase (see Curtis et al. 2010 for a recent summary). The survey of Bontemps et al. (1996) studied low-mass protostars in the same 12 CO J = 2 → 1 transition that we used in this survey toward a sample that spanned bolometric luminosities from 0.2 to 41 L bol among sources from a heterogenous population of cores and clouds. Bontemps et al. found a correlation between the outflow force and bolometric luminosity with Class 0 sources having larger outflow force than Class I sources.
Recently, Curtis et al. (2010) performed a similar analysis on the protostellar population within a single cloud (Perseus) with sources extending to lower luminosities and found similar results for Class 0 and I sources. The slopes and offsets of the Curtis correlations differ from Bontemps correlations, but the trend of the Class 0 outflow force correlation lying above the Class I correlation holds (see Figure 19).
We have plotted the outflow force of the protostars observed in this paper along with VeLLOs and low-luminosity protostars from the Dunham catalogue that were observed in the literature (Figure 19). The Class 0 (solid line) and Class I (dashed line) correlations observed by Curtis et al. (2010) are over-plotted. The three single-dish detections in this paper are generally more consistent with the Class 0 correlation than the Class I correlation. In general, however, when VeLLO outflows from the literature are included, single-dish detections span both the Class 0 and Class I correlations. The VeLLOs that are only detected interferometrically (blue squares are L1014-IRS and L1148-IRS) are more consistent with the Class I correlation. This is a fair comparison since outflows from both VeLLOs were searched for and not detected with single-dish telescopes (Crapsi et al. 2005;Kauffmann et al. 2011, this work -L1148 (032) ). This result is consistent with the conclusions drawn by Shirley et al. (2007) that L1014-IRS is a more evolved protostar which was detected with Spitzer while in a low accretion state.
The time average accretion rate, time-averaged accretion luminosity, total outflow mass, and dynamical time of the outflow are also plotted in Figure 19. All of the detected VeLLO outflows are less than 10 5 years old. There is a strong positive correlation between the observed outflow mass and the observed dynamical time of the outflow. More luminous sources also tend to have higher mass outflow rates (and therefore higher accretion rates) although the correlation is less well defined. This is not surprising since the accretion luminosity is linearly proportional to dM/dt.
The time-averaged accretion luminosity onto the protostar is larger than the observed internal luminosity for every single-dish detected protostar except for one source (MMS 126). The median ratio of the outflow-derived, time-averaged accretion luminosity and the observed internal luminosity of the protostar at the current epoch (current meaning when observed by Spitzer ) for all VeLLO outflows detected with single-dish telescopes is 1.8. The average value of the ratio is 17 but is significantly skewed by IRAM04191 which has a ratio of 89 (André et al. 1999). This result is direct evidence for time-variable accretion. As we shall discuss in the next section, there are several caveats that must be addressed, but the L acc /L int values plotted are likely lower limits and should only increase when corrections for optical depth, excitation, and inclination ( Figure 12) are made.
These results may be directly compared to theoretical models of mass accretion. One such model is an episodic accretion for VeLLOs where accretion bursts are short-lived events precipitated by instabilities in the protostellar disks. The recent paper of Dunham & Vorobyov (2012) investigate this scenario using hydrodynamic models of disks with gravitational instabilities (Vorobyov & Basu 2010) coupled with radiative transfer models to analyze the time evolution of the observed protostellar core. In reality, the instability does not have to be purely gravitational in origin but also may be magnetic in origin (e.g. MRI Armitage et al 2001, Zhu et al. 2010, Simon et al. 2011Tassis & Moushovias 2005). The general nature of the variation of the mass accretion rate with time is an exponentially decreasing function punctuated with short-lived bursts that vary by two orders of magnitude (see Figure 6 of Dunham & Vorobyov 2012). The average mass accretion rate in their example model is 3.7 × 10 −6 M ⊙ /yr over the first 10 5 years after the formation of the protostar. This is a factor of 4.2 times higher than the typical accretion rate in their example model at 10 5 years. These numbers are very similar to the observed ratios of the accretion luminosity derived from the outflow properties to the observed protostellar luminosity in this paper. This is only a single example model and models with a range of differing initial conditions would produce a wider range of average accretion rate ratios. With the caveats in mind which indicate that the observed L acc /L int and therefore <Ṁ > /Ṁ now is likely underestimated, this result is consistent with a picture of episodic accretion; but other smoothly varying accretion models, such as a model with a purely exponential drop in the accretion rate with time, cannot be ruled out by these observations alone. Advantages of the episodic accretion model are that it is physically motivated and that it is also successful in explaining the distribution of protostars in a bolometric luminosity -bolometric temperature diagram (i.e. it physically explains the original "luminosity problem"; see Figure 8 of Dunham & Vorobyov 2012). The results from this outflow survey provide additional supporting evidence but not definitive evidence of the episodic accretion model.

Caveats
Despite our sensitive maps in 12 CO and 13 CO J = 2 → 1, it may still be very difficult to distinguish the terminal bow shock of the outflow in a complex kinematic environment. Indeed, many of our protostellar candidates are situated in very confusing regions where the observable signature of a weak outflow can easily be lost on large scales. In this sense, the dynamical time is a lower limit to the actual age of the outflows. Parker (1991) and Fukui (1993) have statistically studied outflow ages and compared to the dynamical time and found an up to order of magnitude discrepancy (see Curtis et al. 2010 for a discussion). The impact of this underestimate of dynamical time on the outflow force is more complicated because more mass would also be swept up by the larger outflow. Therefore, the outflow force is unlikely to decrease in direct proportion to the increase in dynamical time. The effect of a underestimate of the outflow terminal bow shock on the estimates of the average accretion luminosity is also more complicated. If the outflow extends beyond the region estimated from our maps, then the larger outflow age means that more mass has been accreted. However, that change coupled with an increase in the estimate of the dynamical age means that the overall change in the average accretion luminosity may not change drastically from what we have estimated in this paper.
This caveat must be weighed against the uncertainty in the outflow optical depth, excitation temperature, and inclination which will only increase the derived outflow quantities. For instance, if the outflow inclination is increased from 34.3 • to 70 • , then the outflow force would increase by a factor of 10. Typical corrections for higher excitation temperatures or optical depths are additional factors of a few. Throughout the calculation of outflow parameters in §3, the lower limits were consistently chosen and L acc was calculated for i = 34.3 • ( sini cos 3 i ≈ 1) which corresponds to an 83% confidence interval that the real outflow inclination is larger. Given these factors, it seems very unlikely that the outflow forces plotted for VeLLOs in Figure 19 will decrease.
We find at least factors of 1.5 and 2 difference between the average accretion rate estimated from outflow properties and the observed protostellar luminosity. These ratios are significant at the > 2σ level for all three sources. The typical optical depth corrections (factors of a few) and larger inclinations (also factors of a few) will only tend to increase this difference. The accretion luminosity itself is estimated from very rough values of the outflow jet speed, entrainment efficiency, mass loss rate in the wind versus the mass accretion rate, and assumed protostellar radius. A simultaneous change in several of the the assumed parameters listed above (e.g. lowering the typical jet speed from 150 km/s to 75 km/s, increasing the ratio ofṀ ẇ Macc = 0.2, and increasing the protostellar radius from 3 to 9 R ⊙ ) would be required to reconcile the two luminosity estimates for typical optical depth, excitation temperature, and inclination effects.

Conclusions
We have systematically surveyed 30 VeLLOs and 9 low-luminosity objects in 12 CO and 13 CO J = 2 → 1 with the HHT from the Dunham et al. (2008) catalog of candidate lowluminosity protostars in nearby (D < 400 pc) clouds and cores. From five-point maps, we identified five sources as new candidates outflows and made extended on-the-fly 12 CO J=2 → 1 maps of these regions. Only B59 (023) showed evidence of a candiate outflow, despite the red wing being obscured by emission from IRAS 17081-2721. The remaining sources do not show convincing evidence for a large scale outflow in 12 CO maps. When combined with the grid maps for all 39 sources, these results indicate that the distinct kinematic signatures of large-scale outflows are rare (13% single-dish detection rate for VeLLOs) toward sources in the Dunham et al. (2008) catalog. Every candidate protostar from the Dunham catalog with a single-dish detected outflow belongs to the highest confidence group classifications (groups 1, 2, and 3). We also re-mapped two previously detected VeLLO outflows (L673-7 (031) and L1251-A (045)) to find the full extent of the outflows and better determine their outflow properties. Our results indicate that the time averaged accretion luminosity implied by the outflow properties is at least a factor of two larger than the observed internal luminosity. For larger inclinations of the outflow axis, higher gas excitation temperatures, and optical depth corrections this discrepancy will increase. This result is consistent with a picture of episodic accretion for these objects.
In light of the results from our single dish survey, interferometric observations are needed to complete the search for VeLLO outflows. Interferometric observations have the advantage of spatially filtering out much of the large scale structure that is potentially causing confusion in our single-dish survey. The combined results from an interferometric survey and this single dish survey will provide a census of outflow properties and statistics for VeLLOs that will illuminate the history of accretion towards these sources. Tassis, K., & Mouschovias, T. C. 2005, ApJ, 618, 783 Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896Young, C H., et al. 2004, ApJS, 154, 396 Zhu, Z., Hartmann, L., Gammie, C. F., et al. 2010, ApJ, 713, 1134 This preprint was prepared with the AAS L A T E X macros v5.2.                b Group indicating likelihood of being an embedded protostar.
c All source names are preceded by the prefix "SSTc2d ".
d Name of the core region from Evans et al. (2007) in which this source is located.
e Integrated luminosity using all available detections between 1.25 (2MASS J-band) and 70 µm. Entries in italics denote that the calculation only extends to 24 µm due to no flux information available at 70 µm.  Table 3. Spectral properties for the center position of each source.
b Ellipses indicate non-detections.