\section{Introduction} \label{intro}

The incredible sensitivity of the \textit{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$ \lsun .  The initial \textit{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
$L_{bol} \sim 0.004 \left(\frac{D}{140 pc}\right)^2$ \lsun\ and
can identify essentially all low-mass protostars with accretion
rates down to 
$\frac{dM_{acc}}{dt} = 1.3 \times 10^{-9} \frac{R_{proto}}{1 \mbox{$R_{\odot}$} } \frac{0.1 \msun }{M_{proto}} 
\left( \frac{D}{140 pc}\right)^2$ \msun /yr (see Table 1 for the adopted
distances to regions observed).  
The subset of embedded low luminosity protostars ($L_{bol} \leq 1$ \lsun ) 
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$ \lsun ; Dunham et al. 2006) 
driving a large-scale molecular
outflow (Andr\'e 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 (Bourke et al. 2006), 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
\textit{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 \rightarrow 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 \S2 we describe the observing techniques and our 
sample selection.
In \S3 we analyze the properties of 3 detected outflows in the sample as
well as describe the properties of the numerous non-detections.  In \S4 we  
discuss the nature of the VeLLO outflows.


\section{Observations} \label{obs}

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 \rightarrow 1$ and $^{13}$CO\  J = $2 \rightarrow 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$\mbox{\arcsec} . 
Position switching observations with a $5^{\prime}$ OFF position
in azimuth were used to calculate the main beam efficiency
($T_{mb} = T_A^* / \eta_{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 $\eta_{pol}$ likely reflect calibration variations and changes 
in the coupling between the receiver optics and the telescope between observing
seasons.


\subsection{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)
(i.e. B35A (007)).
The Dunham catalog contains 199 sources.
All isolated sources ($> 30$\mbox{\arcsec}\ from a IRAS source) with suitable declinations
observable from the Heinrich Hertz Telescope ($\delta > -30$\mbox{$^{\circ}$} ) 
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\mbox{\arcsec} , 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\'e 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 $\leq 30$\mbox{\arcsec}\ resolution in transitions of $^{12}$CO\ (Hatchell et al. 2007, Curtis et al. 2010;
Burleigh \& 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$ \lsun\ and all sources
have $L_{IR} < 0.45$ \lsun .


For each source in Table 1, 
a 5-point map was obtained with 30$\mbox{\arcsec}$ 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\footnote{The SkyView Virtual Observatory http://skyview.gsfc.nasa.gov}. 
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\footnote{http://iram.fr/IRAMFR/GILDAS/}.


\subsection{OTF Mapping} \label{otf}

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 (\S3.1). 
For each of these sources we created an on-the-fly (OTF) map 
with dimensions 5$\mbox{\arcmin}$ by 5$\mbox{\arcmin}$, 4$\mbox{\arcmin}$ by 4$\mbox{\arcmin}$, or 3$\mbox{\arcmin}$ by 3$\mbox{\arcmin}$.  
For the L673-7 (031) outflow two regions of size 
3$\mbox{\arcmin}$ by 3$\mbox{\arcmin}$ offset by $\pm$120$\mbox{\arcsec}$ in the directions 
of the red and blue wings were also mapped. 
All maps had 10$\mbox{\arcsec}$ row spacing with a scan rate 
of 10$\mbox{\arcsec}$ 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$_{\mathrm{pol}}$ and H$_{\mathrm{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\mbox{\arcsec}\ for $^{12}$CO\ and  36.9\mbox{\arcsec}\ 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 ($\delta 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 ($\mathrm{\sigma_{I} = \mean{\sigma_{T}}\sqrt{\delta V_{chan}(v_{2}-v_{1})}}$ where $\mathrm{\mean{\sigma_{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.


\section{Results} \label{results}

\subsection{Detections \& Non-detections} \label{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 $\Delta v = 2.12 \pm 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$\mbox{\arcsec} )
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$\mbox{\arcsec}\ 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$\mbox{\arcsec} ) 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 \mu$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 (\S3.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.


\subsection{The Fully Mapped L673-7 (031) Outflow} \label{L673}


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\sigma_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$\mbox{\arcsec}$ compared to the 150$\mbox{\arcsec}$ used by Dunham et al. (2010a).
We calculate $\mathrm{t_{d}=4.44\pm0.23\times\eten{4} \, \frac{\cos i}{\sin i}\ \mathrm{yr}}$ assuming an
angular uncertainty of half the beam size (15\arcsec ).
Using this value, we find the time averaged mass accretion rate is
$\mathrm{\langle\dot{M}_{acc}\rangle=2.78\pm0.14\times\eten{-7} \, \frac{\sin i}{\cos^{2} i}\ \msun\ \mathrm{yr}^{-1}}$ resulting
in a total protostellar mass accreted of $\mathrm{M_{acc}= \dot{M}_{acc}t_{d}} = \mathrm{1.23\pm0.09\times\eten{-2} \, \frac{1}{\cos i}\ \msun}$.
The resulting lower limit on the time averaged accretion luminosity ($L_{acc} = GM_{acc}\dot{M}_{acc}/R$ 
with $R = 3$\mbox{$R_{\odot}$} ) is
$\mathrm{L_{acc}=0.036 \pm 0.003 \, \frac{\sin i}{\cos^{3} i}\ \lsun}$.  
The ratio of the outflow-derived accretion luminosity to the
internal luminosity is $2.1 \pm 0.7 \, \frac{\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{\mbox{$^{\circ}$}}, the correction factor 
$\frac{\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 $\frac{\sin i}{\cos^{3} i} \approx 1$  corresponds to an 
83\%\ confidence interval that the real outflow inclination is larger than $i_o$ = 34.3{\mbox{$^{\circ}$}}. 
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 \S4.1).  



\subsection{The L1251 (045) Outflow} \label{IRS4}


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$\mbox{\arcsec}\ versus 
$48$\mbox{\arcsec}\ in Lee et al.) and with higher signal-to-noise. 
We use the same method as for L673-7 to derive the outflow properties 
(\S3.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 $\pm$ 50 pc (see Dunham et al. 2008)
we calculate an outflow dynamical time of 
$\mathrm{t_{d}=2.69\pm0.23\times\eten{4} \, \frac{\cos i}{\sin i}\ \mathrm{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  $\mathrm{\langle\dot{M}_{acc}\rangle=6.3\pm1.5\times\eten{-7} \, \frac{\sin i}{\cos^{2} i}\ \msun\ \mathrm{yr}^{-1}}$. 
The lower limit on the protostellar mass accreted is 
$\mathrm{1.7\pm0.6\times\eten{-2} \, \frac{1}{\cos i}\ \msun}$, 
making the lower limit for the accretion luminosity
$\mathrm{L_{acc} =  0.11 \pm 0.04 \, \frac{\sin i}{\cos^{3} i}  \, \lsun}$.  
The average (minimum) accretion luminosity is a factor of $1.5 \pm 0.7 \, \frac{\sin i}{\cos^{3} i}$ 
times higher than the observed internal luminosity.  For inclinations greater than
34.3\mbox{$^{\circ}$} , the discrepancy will become larger (Figure 12). 


\subsection{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 \rightarrow 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$ \lsun ).  The (minimum) average accretion luminosity is
$1.5 \pm 0.5 \frac{\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.


\subsection{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\sigma$ limit for a mass detection is $1.38\times \eten{-3} \msun$ which 
is consistent with the uncertainties for the detected outflows. We display the 
maps of these non-detections for completeness.



\section{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$ \lsun ) although
many sources with luminosities approaching $0.1$ \lsun\ 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$ \lsun ).  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\'e et al. 1999, 
MMS 126 Stanke et al. 2006, L1014 Shirley et al. 2007, L1448-IRS2E1 Chen et al. 2010,  
Per 071 073 105 106 Curtis et al. 2010, L1251A-IRS3 Lee et al. 2010,
L1148-IRS Kauffmann et al. 2011, Per-Bolo-58 Dunham et al. 2011, 
L1451-MM Pineda et al. 2011, LFAM 26 Nakamura et al. 2011, CB17 MMS Chen et al. 2012).  
Not all quantities plotted in 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\mbox{$^{\circ}$}\ (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 \rightarrow 1$ transition that we used in this survey
toward a sample that spanned bolometric luminosities from $0.2$
to $41$ \lbol\ 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 \textit{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 \textit{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\'e 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. 2011; Tassis \& 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 \times 10^{-6}$ \msun /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
$<\dot{M}>/\dot{M}_{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.  




\subsection{Caveats}


Despite our sensitive maps in $^{12}$CO\ and $^{13}$CO\ J = $2 \rightarrow 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$\mbox{$^{\circ}$}\ 
to $70$\mbox{$^{\circ}$} ,
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 \S3, the lower limits were consistently
chosen and $L_{acc}$ was calculated for $i = 34.3$\mbox{$^{\circ}$}\ 
($\frac{sin i}{cos^3 i} \approx 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\sigma$ 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 $\mathrm{\frac{\dot{M}_{w}}{\dot{M}_{acc}}=0.2}$,
and increasing the protostellar radius from 3 to 9 $R_{\odot}$) 
would be required to reconcile the two luminosity estimates for typical
optical depth, excitation temperature, and inclination effects. 










\section{Conclusions}

We have systematically surveyed 30 VeLLOs and 9 low-luminosity objects
in $^{12}$CO\ and $^{13}$CO\ J = $2 \rightarrow 1$ with the
HHT from the Dunham et al. (2008) catalog of candidate low-luminosity 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 $\rightarrow$ 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.


\section*{Acknowledgments}

We thank the staff of the Arizona Radio Observatory including the
operators Bob Moulton and John Downey for their assistance with
the observations.  Partial funding for KS was provided by the
Arizona Space Grant Consortium.   YLS is partially supported by
NSF grant AST-1008577.



\newpage

