\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


\section{The Minkowski space-time model: does it really represent
the ultimate physical geometry?}

Nowadays, the Minkowski geometry necessarily constitutes the
framework  of any realistic physical geometry including
its curved version in General Relativity, multidimensional schemes etc.

However, the true dimension and structure of real space (and time) is still
unclear, and there is still no answer to the old question: ``Why {\it
``visible''} physical space is 3-dimensional?'' And what in fact
can be said on the origin and properties of {\it Time}?

From general physical viewpoint, the Minkowski space $\bf M$, together with
the Lorentz group as its group of symmetry, possesses a number of
remarkable and well-known properties (causal
structure, universal velocity of propagation of interactions etc.),
is well grounded in theoretical (symmetries of Maxwell and Lorentz equations)
and purely experimental respects.

On the other hand, there are some known drawbacks of the STR paradigm
(absence of the global evolutionary parameter, of phase relations,
of time irreversibility etc.). And the most grievous fact is that
$\bf M$ {\it does not originate from first principles}, its
mathematical structure is rather cumbersome and not distinguished
in any aspect. Besides, assumption  of $\bf M$ as physical
geometry {\it does not offer any key to comprehend the dynamics of fields and
particles or the quantum properties of matter!}

Thus, one has to return back to foundations and seek for an
internally exceptional structure of (extended) physical space-time which
by itself could predetermine the properties of matter!

In this note, we briefly present two approaches to the
establishment of the extended space-time structure and,
combining them, demonstrate a purely classical, geometrical
explanation of the quantum interference phenomena.

First approach (section 2) is the ``6D treatment of Special Relativity''
proposed by Igor A. Urusovskii in 1996-2011 (see~\cite{UrusPIRT} and references
within). Key notions of his approach are: 1) universal
light-like velocity of {\it any} particle in the whole 6D space, and
2) universal internal rotation of particles in the 3D space
orthogonal to the physical 3D one.

Second approach (section 3) is my own {\it ``algebrodynamical''} one (see
references, say, in~\cite{PIRT03,PIRT09}) in which one is also brought to a 6D geometry of
special type. Key constituents arising naturally in the framework
of algebrodynamics are (sections 4-6) 1) the ensemble of identical copies
({\it``duplicons''})
of one and the same particle and, 2) {\it ``dimerous''} structure of
electron which turns to be detectable only at the instants of
merging of some two duplicons.

Finally, rather simple and natural combination of the above two approaches
(section 7) immediately follows in a clear classical picture of the canonical
{\it two-slit interference experiment} (section 8). Concluding remarks are
presented in section 9.

\section{On the 6D geometry of Urusovskii}

Approach presented by I.A. Urusovskii in his original papers~\cite{UrusClass}
is based in fact on the following two axioms (conjectures):

1) {\it universal light-like motion} of all of the  pre-elements of matter
(say, electrons) in the extended (3+3) physical space (so that {\it none particle
is ever at rest}!). The principle was earlier discussed by F. Klein, Yu.B. Rumer,
R.O. di Bartini et al and can either substitute or supplement the principle of
{\it universal velocity of interaction propagation} in Special Relativity.

2) {\it universal periodical regular rotation} of particles along a
curcumference of {\it Compton radius} $R=\hbar/Mc$  in the additional 3D-space
(orthogonal to the ordinary physical 3-space). This conjecture closely relates to
the L. de Broglie's hypothesis on {\it internal frequency} $\nu$
(internal clocks) of elementary particles, $\hbar \nu = Mc^2$, to the
conception of {\it Zitterbewegung} of E. Shr\"odinger etc.

Consequently, Urusovskii naturally assumed that {\it proper time} of a particle (being,
say, at rest in ``observable'' 3-space and rotating, therefore, with constant
speed $c$ in orthogonal 3-space) {\it is proportional to the number of turns
(or to phase increment) in its  rotation in the internal space}.

Then for a particle moving uniformly with velocity $v$, velocity of rotation
$u$ is obtained from the principal relation $v^2+u^2=c^2$ and determined
by the familiar Einstein's factor,
\begin{equation}\label{decel}
u=c\sqrt{1-v^2/c^2}.
\end{equation}
Thus, one comes to the canonical effect of {\it time deceleration} for a moving
particle, even {\it apart of any specification of a transformation symmetry
group}.

As to the latter, Urusovskii himself assumed it to be the Lorentz group
since ``the number of turns = the proper time interval'' should be invariant,
the same for any observer.

Moreover, according to Urusovskii, {\it all principal relations of STR
including the dynamical ones, follow in fact from purely Newtonian physics in
the entire 6D space, under projection onto the ``observable'' 3D space}!
For more details of Urusovskii's approach see his original papers~\cite{UrusClass};
on further development of the theory see, e.g.,~\cite{UrusPIRT} and references
therein.

Let us remark now that the first postulate, on the universal light-like
velocity, is obviously equivalent to the following relation for the
increments of particle's  coordinates $\{{\bf X,Y}\}$ in the physical ${\bf X}$ and
internal ${\bf Y}$ spaces:
\begin{equation}\label{diff}
dx_1^2+dx_2^2+dx_3^2+dy_1^2+dy_2^2+dy_3^2= c^2 dt^2.
\end{equation}
We see thus that {\it one deals in fact with a Euclidean 6D space}, whereas
the {\it global universal time should be considered as the metric on this space},
not as a coordinate! In this respect, the Urusovskii's viewpoint on his theory
as on a ``6D treatment of STR'' is, at least, not evident or uniquely possible.
Note also that recently a number of attempts were undertaken {\it to
reformulate STR on the basis of the 4D Euclidean space} (see, e.g.,~\cite
{Montanus,Pestov}).

As to the second postulate on regular internal rotation (in a 2D plane),
it can be implemented for a particle at rest ($dx_1=dx_2=dx_3=0$)
if one requires, say, in addition to (\ref{diff}):
\begin{equation}\label{rotat}
dy_3=0,~~~dy_1^2+dy_2^2=ds^2=R^2 d\varphi^2 = c^2 dt^2 ,
\end{equation}
where the proper time increment $ds$ turns to be the metric on {\it internal}
space, whereas $d\varphi$ is the corresponding phase increment.

From further results obtained by I. Urusovskii in the framework of his 6D approach
one can distinguish, say, 1) the {\it hydrogen spectrum} with account of its
fine structure, 2) the {\it quark model} of nucleons, 3) a new treatment of
gravitation (as a projection of {\it cosmological force} retaining particles
in a Compton-order vicinity of ${\bf R}^3$), 4) a novel treatment of the
{\it Universe expansion} (resolving a number of paradoxes in standand
cosmology) etc.

Despite of these results (which one can give credence to or not), it is
especially important that Urusovskii's approach allows {\it to link the
space-time geometry with phase relations} (what could probably lead to a
{\it geometrization of quantum theory}) and to directly deduce universal
{\it kinematics} of particles (and, to some extent, their {\it dynamics} as
well) from pure {\it geometry}!

However, in many aspects the Urusovskii's approach is certainly insufficient.
Indeed, any general substantiation of the starting postulates is absent,
the problem of {\it mass spectrum} is not resolved, structure and kinematics
of {\it photons} (which, in the considered framework, should be treated as
non-rotating in ${\bf Y}$) remain quite unclear etc. Evidently, a more general and
mathematically profound approach is needed. In the following section we
briefly present the algebrodynamical scheme that can aspire to such a role.

\section{Main principles and results of algebrodynamics}

In the {\it algebrodynamics}~\cite{AD,Qanalys} one assumes the existence of a
({\it exeptional} in its internal properties) {\it space-time algebra} (STA).
The STA predetermines both the geometry of physical space-time (via the {\it
automorphism group} acting as a Klein's symmetry group of metric), and the
particle-field dynamics (via the {\it analyticity conditions} for
functions-fields over the STA acting as field equations).
Particles themselves may be confronted with {\it singularities}
of functions-fields.

In the role of STA it had been considered the {\it algebra of complexified
quaternions} (biquaternion algebra $\mathbb B$), equivalent to the matrix
$Mat(2,\mathbb C)$-algebra.
The $\mathbb B$-analyticity  conditions -- generalized {\it Cauchy-Riemann
equations} (GCRE) -- {\it turn to be nonlinear as a consequence of noncommutativity}.
Thus, one deals in fact with fields with (self-)interaction!

Some consequences of GCRE are: emergence of gauge and spinor (twistor)
structures, identical satisfaction of Maxwell and Yang-Mills free equations
(for secondary derivative fields), {\it self-quantization} of electric charge
~\cite{Sing} etc.

Note also that the geometry induced by the structure of $\mathbb B$-algebra is
4D in complex numbers = 8D in reals. Whether one (rather artificially!)
reduces the coordinate space to real ${\bf M}$, the whole theory becomes
{\it Lorentz invariant}!

A novel concept of physical Time had been also elaborated in which time
manifests itself as a evolution parameter related to the
condition of {\it local preservation of the primordial $\mathbb B$-field}
(= principal spinor or twistor field in other equivalent representations).

Local evolution turns to be a {\it light-like transfer} of the primary
$\mathbb B$-field. On this way, one comes to the concept of the {\it Flow of
PreLight = Flow of Time}~\cite{Pavlov,Levich}.
Geometrically, particles manifest themselves as {\it caustics of the PreLight
Flow} so that {\it whole matter is of a light-like nature}!

\section{Complex-quaternionic geometry of space-time}

It is noteworthy, however, that none algebra is known which could give rise
precisely to ${\bf M}$. Nonetheless, the considered 4D (in $\mathbb C$) = 8D
(in $\mathbb R$) ${\mathbb B}$-induced geometry has the $SO(3,\mathbb C)$
symmetry group which is 6-parametric and isomorphic to the Lorentz group!

On the other hand, ${\bf M}$ does not contain any {\it phase-like} component,
whereas the ${\mathbb B}$-induced geometry naturally decomposes into a
Minkowski-type {\it macro-geometry} (defined by the modulus parts of complex
coordinates) and {\it micro-geometry} of a ``fiber'', related to their phase
parts~\cite{Mink05} and responsible, probably, for the wave properties of
matter.

Actual ${\mathbb B}$-geometry of (extended) space-time turns to be 6D and
resembles the geometry proposed by I. Urusovskii. It is dynamical in origin
being closely related to universal kinematics of particles in the primordial
complex ${\mathbb B}$-space. Such kinematics can't be realized in real Minkowski
space and results in a rather unexpected picture of Physical World.
One can thus assume that {\it we really live in the complex space-time}!

\section{Complex null cone and the concept of duplicons}.

Complex space-time geometry allows for realization of the {\it one-electron Universe}
conjecture of Wheeler-Feynmann.  On the {\it real} space-time background this is in
fact impossible. Indeed, in real ${\bf M}$, according to the {\it retardation
equation (= equation of the light cone)}, one has only one unique object in the
past that determines the field at a given point  (for $v<c$; for tachions --
more than one, see, e.g.,~\cite{Bolot}).

On the contrary, in complex space $Z=\{z_0,z_1,z_2,z_3\}$ equation of the
{\it complex null cone}
\begin{equation}\label{complexcone}
(z_1-{\hat z}_1(\tau))^2+(z_2-{\hat z}_2(\tau))^2+(z_3-{\hat z}_3(\tau))^2=
(z_0-\tau)^2, ~~~\tau \in {\mathbb C}
\end{equation}
can have a great number of roots (dependently on the form of
${\mathbb C}$-trajectory ${\hat Z}(\tau)$),
in accord, say, with the ``principal theorem of algebra'' (for
polynomial functions ${\hat Z}(\tau)$).

\begin{figure}
\centering
\includegraphics[width=10cm,height=5cm]{Caustics.png}
\caption{\small Duplicons and caustic formation.}
\end{figure}

Therefore, an observer ``living in a ${\mathbb C}$-space'' instantaneously
receives the field signals from a lot of copies locating at the points on the
same ``complex worldline'' (Fig.1). Thus, one actually deals with an ensemble of
identical particle-like formations -- ``duplicons''~\cite{PIRT05,YadPhys}.
Points of observation and location of a duplicon are connected by a $\mathbb C$-null
straight line -- ``ray'' (quite similar to the case of real ${\bf M}$ where
one deals with local light cones).
These lines densely filling ${\mathbb C}$-space form a
fundamental {\it null congruence} for which the ${\mathbb C}$-worldline itself
play the role of its {\it  focal line}~\cite{PIRT05}.

\section{Caustics as lightlike signals and ``dimerous electron''}

At some {\it discrete} instants (for a corresponding position of an observer
in the ${\mathbb C}$-space) a {\it merging} of a pair of duplicons does take
place (Fig.1). This relates to {\it multiple roots} of the ${\mathbb C}$-null cone
equation (\ref{complexcone}). Again, such an event is in principle impossible
in real ${\bf M}$.

At these moments strong {\it amplification} of electromagnetic and other fields
emerges, along a  null straight line connecting points of observation and
of duplicons' merging. Geometrically, these lines represent the {\it caustic
locus} of the fundamental null congruence~\cite{PIRT05}.
Physically, one observes a  {\it lightlike pulse -- a signal} coming from
the points of duplicons' mergings.

Thus, in the framework of a purely classical electrodynamics but contrary to the
case of real ${\bf M}$ background, in complex world there occure
{\it quantum-like discrete radiation processes}!

One obtains therefore a {\it self-consistent relativistic dynamics of an
ensemble of identical and causally connected ``particles''}.
It is evident, however, that (contrary to the Wheeler-Feynmann paradigm)
individual duplicon {\it in principle} can't be considered as a primary
matter pre-element: only a pair of these can be detected, and only at discrete
instants of their merging. We are forced thus to surmise that ``electron''
not only {\it consists} of two ``halves'' - duplicons, but does not in fact
{\it exist} as a unique formation, except at discrete instants of
merging - radiation - detection acts.

Conjecture on {\it electrons as dimerons}~\cite{PIRT09,YadPhys} directly
correlates with modern notions on ``fractal charges'' and ``electron bubbles''~\cite{Bykov},
``light-shining-through-wall'' effects~\cite{tunnel} etc. It allows for
{\it alternative explanation of the wave  properties of matter} (section 8).

\section{Links between complex and Urusovskii's geometries}

In full analogy with the case of real ${\bf M}$, the primary field and the
caustic structure both {\it reproduce themselves along the null complex straight
lines} --  elements of the $\mathbb C$-null cone~\cite{PIRT05}.
For their infinitesimals one has
\begin{equation}\label{inf}
dz_0^2+dz_1^2+dz_2^2+dz_3^2=0.
\end{equation}
Decomposing $dz_\mu = dx_\mu + i dy_\mu$ and separating real and imaginary
parts in (\ref{inf}), one obtains a system of {\it two quadrics} in
${\bf R}^8$,
\begin{equation}\label{quadric}
dx_0^2+dx_1^2+ dx_2^2+ dx_3^2-dy_0^2-dy_1^2-dy_2^2-dy_3^2 =0,
\end{equation}
\begin{equation}\label{quadric2}
dx_0 dy_0 +dx_1 dy_1 +dx_2 dy_2 +dx_3 dy_3 =0 .
\end{equation}
First constraint (\ref{quadric}) defines in fact {\it two 4D Euclidean
(sub)spaces}; these are moreover mutually orthogonal as a result of the second
equation (\ref{quadric2}). We can then define the {\it time increment} $cdt$ as
{\it common metric} on the both subspaces, so that relation (\ref{quadric})
takes the form:
\begin{equation}\label{timemetric}
dx_0^2+dx_1^2+ dx_2^2+ dx_3^2 = dy_0^2+ dy_1^2+ dy_2^2+ dy_3^2 := c^2 dt^2.
\end{equation}
It is natural to assume that fundamental kinematics of matter pre-elements
follows just from the condition of preservation of the primary field and
caustic structure and
satisfies therefore condition (\ref{inf}) or, equivalently, (\ref{quadric2})
and (\ref{timemetric}). According to these, {\it universal ``eternal'' motion
of matter pre-elements occures in both mutually orthogonal spaces $\bf X$ and
$\bf Y$, with one and the same, constant in modulus fundamental
``velocity of light''} $c$.

Thus, complex geometry leads by itself to satisfaction of the first postulate
of the Urusovskii's scheme. Specifically, for a particle which is at rest in
the  ``physical'' 3D subspace of $\bf X$, one obtains making use of
(\ref{quadric2}) and (\ref{timemetric}):
\begin{equation}\label{3Drest}
dx_1=dx_2=dx_3=0,~~\rightarrow~~dy_0 =0,~~~
dx_0^2= c^2 dt^2 = dy_1^2+ dy_2^2+ dy_3^2,
\end{equation}
so that it necessarily moves with the speed of light in the 3D subspace
of $\bf Y$.

The geometry allows moreover for exact {\it rotating configurations}
considered by Urusovskii. To obtain these, one can require, say, in addition
to (\ref{3Drest}):
\begin{equation}\label{rotatconf}
dy_3 =0,~~dy_1^2 +dy_2^2 = R^2 d\varphi^2~~(= c^2 dt^2 = dx_0^2),
\end{equation}
We see therefore that the increment of zeroth coordinate $dx_0$ (generally,
of $dz_0 = dx_0+i dy_0$), invariant under transformations from the symmetry group
$SO(3,\mathbb C)$ of ${\mathbb B}$-algebra, appears in the role of a
{\it proper time interval}~\cite{Mink05}, in full
correspondence with the Urusovskii's picture. It is clear, however, that
much more general configurations are admissible, namely (for a particle
resting in 3D subspace ${\bf X}$), an arbitrary light-like motion in the
orthogonal 3D subspace ${\bf Y}$.

\section{Phase of internal revolution and quantum interference}

For Urusovskii's configurations applied to a duplicon its {\it phase change}
$d\varphi$ along a trajectory is proportianal to the increment of proper time $ds$,
\begin{equation}\label{phase}
d\varphi = ds/R,
\end{equation}
where $R$ is the radius of internal rotation which, according to Urusovskii, is
equal to the Compton length of the duplicon. However, for correspondence with
quantum theory (see below), we shall assume for this radius
\begin{equation}\label{radius}
R=\hbar/2Mc,
\end{equation}
where $M$ is taken to be the mass of the particle (electron) associated with
duplicons.

\begin{figure}
\centering
\includegraphics[width=10cm,height=7cm]{QIntrfN.jpg}
\caption{\small Successive mergings of duplicons in a two-slit experiment.}
\end{figure}

In~\cite{PIRT09,YadPhys}  it had been already demonstrated that, together with
the ``dimerous electron'' conjecture, fundamental relation (\ref{phase})
along with (\ref{radius}) leads to a purely classical explanation of the
canonical two-slit experiment. Specifically, one can consider this as two
successive mergings of a pair of duplicons at space-time points where they can
be detected by an observer (Fig.2). Since these mergings occure in fact in the
extended complex space-time, the coordinates of two duplicons should coincide
whereas their phases could differ to $2\pi N,~~N\in {\mathbb Z}$. Thus, the phase
shift $\Delta \varphi$ for a pair of duplicons between their successive
mergings should be
\begin{equation}\label{phaselag}
\Delta \varphi = \Delta \int d\varphi = 2\pi N,
\end{equation}
or, making use of (\ref{phase}) and (\ref{radius}),
\begin{equation}\label{interf}
\frac{2Mc}{\hbar}\Delta \int ds = 2\pi N.
\end{equation}
This is {\it general relativistic invariant condition for two successive mergings
of duplicons} which substitutes the canonical condition for
interference maxima. Indeed, in the first order of non-relativistic
approximation $v/c <<1$, making use of representation
\begin{equation}\label{approx}
ds=cdt\sqrt{1-v^2/c^2}\approx cdt - v dl/2c,
\end{equation}
with $dl=v dt$ being the path length increment,
one obtains instead of the principal condition (\ref{interf}):
\begin{equation}\label{nonrelativ}
\Delta \int \frac{dl}{\lambda(l)} = N,
\end{equation}
where $\lambda(l)$ denotes the de Broglie's wavelength (generally,
variable along a trajectory),
\begin{equation}\label{wavelength}
\lambda(l):= h/Mv .
\end{equation}

Obviously, equation (\ref{nonrelativ}) corresponds to the familiar
condition for interference maxima (for constant velocity it takes the form
$\Delta L = N\lambda$). However, general relativistic condition in
the algebrodynamics is represented just by equation (\ref{interf}) which may be
rewritten in the form
\begin{equation}\label{interf2}
\Delta S = \frac{N}{2} \Lambda,
\end{equation}
where $S$ is the length of a 4D worldline in real ${\bf M}$ for one of
duplicons (between two successive merging events), and $\Lambda=h/Mc$.

\begin{figure}
\centering
\includegraphics[width=10cm,height=5cm]{QIntrfCompr.jpg}
\caption{\small QM vs AD distributions of probabilities in a two-slit
experiment.}
\end{figure}

Thus, in complex algebrodynamics, in the idealized situation, one deals with
a $\delta$-shaped probability distribution for interference pattern since
duplicons can merge and be thus detected only at some discrete space points
specified by the condition (\ref{interf2}). However~\cite{PIRT09}, due to statistical
errors in a real experiment (during preparation, detection processes etc.)
the probability distribution smears and takes, near a single maximum,
a Gaussian-like form (Fig.3, solid line)
\begin{equation}\label{gauss}
\sim \exp(-\Delta x^2/\Lambda^2),
\end{equation}
while the canonical wave-like distribution is of a form (Fig.3, dashed line)
\begin{equation}\label{canon}
\sim \cos^2(\Delta x/\Lambda).
\end{equation}
In vicinity of a maximum these coincide up to the 3-d order derivative.
Nonetheless, the distinction (see Fig.3) could be revealed
in a suitable quantum interference experiment.

\section{Concluding remarks}

We have demonstrated that quantum interference phenomena can be explained  in
a purely geometrical manner whether one adopts two independent fundamental
conjectures  dealing with 1) universal kinematics of matter pre-elements
(rotation with velocity of light in the additional space) and, 2) dimerous
nature of electrons and, probably, other elementary particles
(detectable at some particular instants of merging of their constituents --
duplicons).

We have shown that these two conjectures are compatible with each other but,
unfortunately, first one does not follow directly from the fundamental
algebraic dynamics in complex space, so that Urusovskii's universal kinematics
remains up to now an intuitive though very attractive hypothesis.

Nonetheless, let us mention that, actually, there exist much more rigid
restrictions on kinematics of pre-particles in the complex-quaternionic space
than those revealed in section 7. These are based on twistor
structure of the primary ${\mathbb B}$-field and make the Urusovskii's
conjecture even more realistic.


\section{Introduction}
\label{I}

Transit networks are fascinating objects to study. While each network and every element in the
network needs to be adapted to a specific situation in the urban agglomeration, it turns out that
nonetheless general  construction principles prevail. The latter fact has been evidenced by
numerous studies using in particular methods of statistical physics, random graph theory, and
other approaches. One of the main fields contributing to this analysis is complex network
science \cite{small_world,complex} (see \cite{Berche12,Holovatch11} for a review). Applying the tools
of complex network theory it was shown, that transit networks
(in this paper we consider urban public transportation networks, PTN) share
common features of other natural and man-made complex networks. In particular,
they are strongly correlated compact structures
(so-called small worlds) \cite{small_world}, that often possess power law scale-free behaviour \cite{vonFerber}.
It has been shown, that they are resilient against random failures while being vulnerable to targeted
attacks \cite{Holovatch11,Berche09}. The study of fractal properties of PTN has a long history
\cite{Thibault87,Frankhauser90,Benguigui91,Benguigui92,Kim03}. The prevailing
part of these former studies concern the density of stations or the total length of track
as function of the distance from the center of the network. In this way PTNs and rail
networks of various cities and regions were analyzed.
Thibault and Marchand \cite{Thibault87}  studied the fractal
properties of different subnetworks within Lyon (regions I,II, III). The fractal dimensions
found were for the suburban rail I: 1.64, II: 1.66, III: 1.88; for the public bus service: 1.00, 1.09,
1.45; the authors also analyse the drainage utility with results: 1.21, 1.30, and 1.79.
Frankhauser \cite{Frankhauser90} obtained $1.58$ for
Stuttgart's rail system.
Benguigui and Daoud \cite{Benguigui91} obtained the fractal dimension for
the Paris railway transportation network as 1.47. In relation to this
result they also discuss DLA and star polymers as known models that
lead to fractal structure. Further results for fractal dimensions were  given by Benguigui \cite{Benguigui92}
for several Rhine towns and Moscow railways: $1.70\pm 0.05$,
Paris metro: $1.80\pm 0.05$.
The fractal dimension of Seul PTN has been estimated as  1.50 for stations
and 1.35 for tracks  \cite{Kim03}.
In a complementary approach, the present authors measured the mean distance
as a function of number of stations traveled \cite{vonFerber}.
For different Berlin subnetworks a fractal dimension in the range
of 1.04 -- 1.22 was found. The Paris metro data leads to a
fractal dimension 1.22 when excluding short distance
contributions.

Some of the results cited above appear to be compatible with the interpretation
of the routes as two dimensional self-avoiding walks (SAWs) with a fractal
dimension of $4/3$ \cite{Nienhuis82}. However, the results appear to vary considerably.
Here, in addition to the fractal dimension of the transit routes we intend to analyse
the distributions of inter-station distances of consecutive stations. We will show
that these distributions appear to have heavy (power-law) tails compatible with
a L\'evy-flight model for the transit routes. The latter approach further presents
a possible explanation for the deviations from the SAW behaviour observed.

A L\'evy-flight model is characterised by the distribution function
of the step lengths $\ell$ of the otherwise random successive steps.
For a specific L\'evy-flight the complementary (the 'tail') of the
cumulative distribution function $P(\ell)$ of the step lengths
$\ell$ follows a power law decay
\begin{equation} \label{1}
P(\ell)\sim \ell^{-\mu}.
\end{equation}
where the exponent $\mu$ characterises the distribution. Note that
that exponents greater than 2 belong to the Gaussian Central Limit
Theorem (CLT) domain of attraction due to the fact that for an
exponent $\mu> 2$ the probability density function (PDF) has finite
variance. Exponents $\mu<2$, however, belong to the L\'evy CLT
domain of attraction. The critical exponent $\mu=2$ is more delicate
and requires the analysis of the growth of the associated truncated
variance function.

In the following, we first describe the empirical data as they were drawn from the
Berlin PTN. Then, we analyze the distribution of the inter-station distances along the routes
and propose a description of PTN routes  and journeys in terms of self-avoiding L\'evy flights.
This is supported by the fact that the scaling of distance traveled as function of the number of
stations deviates from SAWs but is in line with L\'evy flight-like behaviour.

\section{Data analysis}
\label{II}

PTNs are often discussed without
reference to their geographical embedding. The fact that this
subject is left aside by most studies of
PTNs with respect to their complex network behaviour, is due mainly
to the lack of easily accessible data on the locations of stations
and routes. For the present work we have analyzed such data for
stations of the Berlin PTN. Our database consists of 2992 stations that belong to 211 routes.
Furthermore, we have taken into account different means of public
transit, that include subway (U-bahn), tram, bus, and surface high
speed transport (S-bahn). The data was collected in the frames of
a larger project \cite{Berche12,vonFerber}  including a comparative analysis of PTN vulnerability
of fourteen large cities  \cite{Holovatch11,Berche09}.
Here, we complement
available information about the Berlin PTN topology by the information
about the geographical location of each single station. To this end,
Cartesian coordinates for each station were extracted from an
on-line map provided by the Berlin public transport operator. We note that using this
procedure slight deviations between distances calculated on the spherical surface and
those derived from the projection to  a flat surface with Cartesian coordinates occur.
These will be neglected below. With the data on the topological
and geographical properties of the PTN at hand we are in the position to analyze how
these properties are related.

A short glance at the map of a public transport network immediately
reveals that transport lines do in general not follow the shortest
"straight ahead" route between their two end points. Thus the
question arises if there is any other underlying structure or
principle characterizing the observed behaviour of geographically
embedded transport routes. Let us analyse transport routes with
respect to their fractal properties in terms of random walks,
self-avoiding walks and L\'evy flights. In former studies
\cite{vonFerber} we have analyzed the dependence of the mean square
distance between stations $\langle R^2 \rangle$ as function of the
number $N$ of stations traveled. The following power law scaling was
found:
\begin{equation} \label{2}
\langle R^2 \rangle \sim N^{2\nu}
\end{equation}
with the exponent $\nu$ ranging from $\nu = 0.82$ for
bus routes to $\nu = 0.9$ and $0.96$ for subway and tram
routes. The S-bahn data is distorted due to a ring structure
within this sub-network. The obtained values of the exponent lie
in the region $3/4 \leq \nu \leq 1$.
These results are to be compared with the behaviour of self-avoiding walks (SAW)
in two dimensions characterised by an exponent $\nu=3/4 = 0.75$ which corresponds
to the exact result for the fractal dimension of $D_{SAW}=4/3$ \cite{Nienhuis82}.
This observation leads to the hypothesis that the PTN routes may follow  SAW statistics.

\begin{figure}
\centering
\includegraphics[width=57mm]{levy_epjst2012_mixed}
\includegraphics[width=57mm]{levy_epjst2012_powerfit}

\vspace*{1ex} {\bf a. \hspace{9em} b.}
 \caption{Cumulative
distribution for Berlin PTN: separate distributions for each means of transportation
({\bf a}) and the sum of these distributions ({\bf b}). The fit of the log-log plot in ({\bf b})
corresponds to an exponent value of $\mu =2 \pm 0.05$.}
 \label{fig1}
\end{figure}

To further analyze the situation,
let us explore in addition the statistics of the inter-station distances.
In Fig. \ref{fig1} we plot the cumulative distribution $P(\ell)$ of the
distances $\ell$ between consecutive stations along the routes showing
individual curves for each mode of transport as well as an integral curve
including all modes (Fig. \ref{fig1}{\bf a}). As follows
from the linear regions of the corresponding log-log plots, the region
of power law behaviour does not extend much more than one order of magnitude
for each single mode of transport. In the combined curve, this power law
regime is expanded.  The slopes and thus the corresponding exponents are
very similar for the bus and S-bahn routes. As one can deduce from  the plot,
the contribution of tram and U-bahn to the combined curve is rather weak.
In Fig. \ref{fig1}{\bf b} we show the fit of the combined curve to the power
law
\begin{equation} \label{3}
P(\ell)\sim \ell^{-\mu}
\end{equation}
with an appropriate exponent value
{\bf $\mu=2 \pm 0.05$}. Our finding is that at least in part
the distribution $P(\ell)$ displays scale-free properties.
In the following, we want to connect these scale-free properties observed for
large scales with the behaviour of the distribution for smaller distances $\ell$.
To this end we project the distance vectors $\underline{\ell}$ to the $x$ and $y$
axes and analyze the probability density functions (PDF) $p_x(x)$, $p_y(y)$.
Under the hypothesis that the tails of the distributions are governed by
a power law we may expect that the distributions $p_x(x)$,
$p_y(y)$ may be described by a L\'evy $\alpha$-stable distribution.
To fix notation we here specify the characteristic function of the
$\alpha$-stable L\'evy distribution
\begin{equation} \label{3a}
\hat{f}(t;\alpha,\beta,\tau,c) = \exp \Big [ -c^\alpha |t|^\alpha \Big (1-\lambda \beta {\rm sign} (t) \tan \frac{\pi\alpha}{2} \Big ) +i\tau t ) \Big ],
\hspace{1em} (\alpha\neq 1),
\end{equation}
with shape, skewness, location, and scale parameters $\alpha$, $\beta$, $\tau$, and $c$.
The asymptotic behaviour of the L\'evy $\alpha$-stable distribution
function $F(x;\alpha,\beta,\tau,c)$ for $\alpha<2$ is then given by \cite{Meyers10}:
\begin{equation} \label{4}
P(x) \simeq c^\alpha d (1+\beta) x^{-\alpha},\hspace{1em} \mbox{with} \hspace{1em} d=\sin (\pi\alpha/2) \Gamma (\alpha)/\pi.
\end{equation}
This relates the exponents, $\alpha=\mu$. Let us emphasize that the above fit concerns
the entire distribution function, and therefore it takes into account the behaviour at
smaller values of $\ell$. This is complementary to fitting
$P (\ell)$, where only the asymptotic (complementary) behaviour is taken into account.

In Fig. \ref{fig2} we fit the empirically determined PDFs $p_x(x)$,
$p_y(y)$ by an $\alpha$-stable distribution \cite{Veillette12}. The
best fits are obtained for values $\alpha=1.71$ for the
$x$-projection and $\alpha=1.67$ for the $y$-projection. These
values are to be compared with the value of $\mu$ derived above.
Recall that for a perfect $\alpha$-stable distribution the relation
$\alpha=\mu$ should hold. The latter holds in our case only
approximately. While one may not expect  high precision from the
given data, the obtained values of $\alpha$ and $\mu$ support the
conjecture that the inter-station distances are governed by a
non-trivial $\alpha$-stable distribution.

\begin{figure}
\centering
\includegraphics[width=55mm]{levy-x-linPDF}
\includegraphics[width=55mm]{levy-x-logPDF}

\vspace*{1ex} {\bf a. \hspace{9em} b.} \vspace*{1ex}

\includegraphics[width=55mm]{levy-y-linPDF}
\includegraphics[width=55mm]{levy-y-logPDF}

\vspace*{1ex} {\bf c. \hspace{9em} d.}

\caption{Linear plots of the probability density functions (PDF) for
all means of transport projected to the $x$- and $y$- axes (left),
lin-log plots of the same (right) fitted to the $\alpha$-stable
distribution with $\alpha=1.71$ for the $x$-plots and $\alpha=1.67$
for the $y$-plots.}
 \label{fig2}
\end{figure}

In further support of the fit to a L\'evy-stable distribution as described above, we
show in Fig. \ref{fig3} the so-called percentile plots for the
cumulative probability density functions (CDF) defined
as:
\begin{equation} \label{5}
P(x)=\int_{-\infty}^x p(x') dx'.
\end{equation}

\begin{figure}
\centering
\includegraphics[width=55mm]{levy-x-percentile25}
\includegraphics[width=55mm]{levy-x-percentile75}

\vspace*{1ex} {\bf a. \hspace{9em} b.} \vspace*{1ex}

\includegraphics[width=55mm]{levy-y-percentile25}
\includegraphics[width=55mm]{levy-y-percentile75}

\vspace*{1ex} {\bf c. \hspace{9em} d.} \vspace*{1ex}

\caption{Percentile plots of the cumulative $\alpha$-stable probability density distribution
function (CDF) in comparison to the data. Distributions are shown in linear scale for the
step-length projected to the $x$-axis (Figs 3a,b) with exponent $\alpha=1.71$ and to the
$y$-axis (Figs 3c,d) with  exponent $\alpha=1.67$. Figs. 3(a,c) display the behaviour of the
25\%-100\% percentile while Figs. 3(b,d) display the 0\%-75\% percentile.
}
 \label{fig3}
\end{figure}

In the light of these results, let us now re-examine the data on the end-to-end mean square inter-station
distance behaviour, cf. Eq. (\ref{2}). As mentioned above,
the power law scaling with the exponent $\nu$ has originally been interpreted in terms
of SAW behaviour in 2d. Note as well, that the
most reliable fit was achieved for the bus subnetwork which also provides the largest data set
and resulted in an exponent $\nu=0.82$.
Let us re-examine this value assuming as a working hypothesis a L\'evy flight
like behaviour for the inter-station distances.  Analytical and numerical results available
for self-avoiding L\'evy flights show that their scaling exponents differ from the usual
SAW exponents \cite{Grassberger85}. In particular, for the  node-avoiding L\'evy flight the
exponent $\nu$ is known in the form of  an $\varepsilon=(2\mu-d)$-expansion
\cite{Fisher72,Halley85}:
\begin{equation}\label{6}
\nu=\frac{1}{\mu} \Big ( 1 + \frac{1}{4\mu} \varepsilon  +
(19-\frac{5}{4}\mu^2)\frac{\varepsilon^2}{64\mu^2}+\cdots \Big ).
\end{equation}
Flory-type arguments on the other hand lead to \cite{Grassberger85}:
\begin{equation}\label{7}
\nu = \left \{ \begin{array}{ll}
 3/(\mu + d), & \quad d<d_c=2\mu \, , \\
 1/\mu, & \quad d>d_c .
\end{array}
\right .
\end{equation}
Substituting  the values $\mu=1.71$, $\mu=1.67$ into Eq. (\ref{6}) we get $\nu=0.80$, $\nu=0.81$,
whereas following Eq. (\ref{7}) we arrive at  $\nu=0.81$, $\nu=0.82$, both
compatible with the observed behaviour of the bus-subnetwork.

\section{Conclusions and outlook}
\label{III}
The present study extends previous observations of fractality within
public transport routes
\cite{vonFerber,Thibault87,Frankhauser90,Benguigui91,Benguigui92,Kim03}
by showing that also the distribution of
inter-station distances along routes displays non-trivial power law behaviour.
This indicates that the routes may in part also
be described as L\'evy-flights. The latter property may result
from the fact that the routes are planned to adapt to fluctuating
demand densities throughout the served area but also possibly related to
optimisation properties \cite{Humphries10}.

One may object, that the distribution derived above mixes different means of
transport with different mean inter-station distances. However, rather observing
journeys of individual passengers one may acknowledge that they will in
general use two or more means of transport to complete their journeys.
The given combined distribution of inter-station distances allows them
to perform this journey such that it resembles a L\'evy flight. This may allow
for travel optimization.

Self-avoiding L\'evy flights, apart from observing the constraint of non
self-intersec\-tion, evolve randomly. The fact that PT routes at least
within  the present sample appear to display the same scaling
symmetry is quite unexpected. In particular, this behavior seems to
be at odds with the requirement of minimizing passengers traveling
time between origin and destination. The latter argument, however,
ignores the time passengers spend walking to the initial and from
the final stations. Including these, one understands the need for
the routes to cover larger areas by meandering through
neighbourhoods. Given the requirements for a PTN to cover a
metropolitan area with a limited number of routes while
simultaneously offering fast transport across the city one may
speculate that routes scaling like self-avoiding L\'evy flights
may present an optimal solution \cite{Holovatch11,vonFerber_un}.

\section*{Acknowledgments}
This work was supported in part by the  7th FP, IRSES project
N269139 ``Dynamics and Cooperative phenomena in complex physical and
biological environments".