OUTFLOW-DRIVEN TURBULENCE IN MOLECULAR CLOUDS

The majority of stars form in self-gravitating cores within giant molecular clouds (GMCs). Two important issues facing modern theories of star formation are the nature of turbulence in these clouds and the relative inefficiency of star formation. For excellent reviews, see Elmegreen & Scalo (2004) and McKee & Ostriker (2007). These two issues are likely to be related. The star formation rate per free fall (SFR_ff) can be defined as the fraction of mass that is converted into stars within one free fall time at the mean density ${\rm SFR_{ff}}=\frac{\dot{M}_{\star }t_{\rm ff}}{M}$ (Krumholz & McKee 2005). In the absence of some form of support, most of the mass in GMCs would collapse into stars within a free fall time t_ff = [3π(32Gρ)]^−1/2 giving an SFR_ff ∼ 1. Observations however find surprisingly low values of SFR_ffs (typically ranging from 0.01 to 0.1). Theoretical accounts for the low values of SFR_ff will rely on some form of support within the cloud to keep it from collapsing to form stars. Supersonic turbulence is one means by which this can be achieved. Turbulence in molecular clouds is inferred from Larson's laws which are empirical relationships between line widths and size observed in many star-forming regions (Larson 1981). Turbulence has also been inferred from direct measurements of power spectra in molecular clouds across a wide range of scales (Heyer & Brunt 2004). While turbulence can provide an isotropic pressure to support the cloud against self-gravity, both hydrodynamic and MHD turbulence decay quickly (Stone et al. 1998; Mac Low 1999). Thus if clouds are stable, long-lived structures supported against self-gravitational collapse by turbulence, those turbulent motions must be continually driven either internally via gravitational collapse and stellar feedback or externally via turbulence in the interstellar medium (ISM).

Feedback in the form of stellar outflows is one means of driving turbulence in molecular clouds (Norman & Silk 1980). If outflows are well coupled to the cloud then observational studies make it clear that combined outflow energy budgets are sufficient to account for cloud turbulent energy (Bally et al. 1996; Bally 1999; Knee & Sandell 2000; Quillen et al. 2005; Warin et al. 1996). Analytical work by Matzner (2000, 2001, 2007) has explored the role of collimated outflow feedback on clouds. Matzner (2007), in particular, developed a theory for outflow-driven turbulence in which line widths were predicted as functions of a global outflow momentum injection rate. Krumholz et al. (2006) have also considered the nature of feedback via outflows, concluding that these systems provide an important source of internal driving in dense star-forming cores. Early simulation studies also indicated that outflows could drive turbulent motions. More recent work by Li & Nakamura (2006) and Nakamura & Li (2007) have mapped out the complex interplay between star formation and turbulence concluding that outflows could re-energize turbulence. However, studies of single jets by Banerjee et al. (2007) came to the opposite conclusion in the sense that single jets will not leave enough supersonic material in their wakes to act as a relevant source of internal forcing.

In this paper, we address the issue of outflow-driven turbulence by focusing solely on the dynamics of cloud material set in motion by multiple transient collimated jets. We use a higher resolution of the outflows than previous studies to tease out the role of their dynamics in driving turbulence. Following the dimensional analysis of Matzner (2007) we consider a cloud of mean density ρ₀, with outflows occurring at a rate per volume $\mathcal {S}$ with momentum $\mathcal {P}$. This then defines characteristic outflow scales of mass, length, and time:

$$\begin{equation} \mathcal {M}=\frac{\rho _{_0}^{4/7}\mathcal {P}^{3/7}}{\mathcal {S}^{3/7}},\quad \mathcal {L}=\frac{\mathcal {P}^{1/7}}{\rho _{_0}^{1/7}\mathcal {S}^{1/7}},\quad \mathcal {T}=\frac{\rho _{_0}^{3/7}}{\mathcal {P}^{3/7}\mathcal {S}^{4/7}}. \end{equation} \tag{ 1 }$$

Combining these gives other characteristic quantities. Of particular interest is the characteristic velocity

$$\begin{equation} \mathcal {V}=\frac{\mathcal {L}}{\mathcal {T}}=\frac{\mathcal {P}^{4/7}\mathcal {S}^{3/7}}{\rho _{_0}^{4/7}}=\frac{\mathcal {P}}{\rho _{_0} \mathcal {L}^3}. \end{equation} \tag{ 2 }$$

Assuming typical values for ρ₀, $\mathcal {P}$, and $\mathcal {S}$ yields supersonic characteristic velocities indicating that outflows contain enough momentum to drive significant supersonic turbulence. Our simulations seek to explore a detailed realization of this idea. We note that these relations are for spherical outflows and provide appropriate order of magnitude estimates for collimated outflows.

This work is the continuation of a series of studies on outflows and turbulence. In Quillen et al. (2005), observations of dense gas in NGC 1333 showed numerous fossil outflow cavities with total energy and momenta consistent with re-energizing turbulence (Cunningham et al. 2006). Most recently Cunningham et al. (2009) have carried out simulations of a single outflow cavity interacting with a turbulent medium. They conclude that outflow-driven cavities are able to re-energize turbulent motions in their immediate environment provided such turbulent motions already exist to disrupt the cavity. Here, we numerically investigate the degree to which interacting outflow cavities are able to provide such a turbulent environment in regions typical of star-forming clusters. Future studies will include treatment of additional physics including magnetic fields (J. J. Carroll 2009, in preparation), self-gravity, and radiation transport. We note that the studies are relevant to broader question of turbulence driven by continuous versus discrete forcing (Joung & Mac Low 2006).

In Section 2, we describe the numerical model and the properties of individual outflows, in Section 3, we discuss the characteristics of outflow-driven turbulence, and in Section 4, we summarize the implications of outflow-driven turbulence within molecular clouds.

In this study, we use a new Riemann solving MHD code called AstroCUB. The code is second order accurate in space and time and uses a nonsplit CTU method as described by Gardiner & Stone (2005). Our current simulations were purely hydrodynamic and were performed on a periodic cube of length 1.5 pc for 1 Myr at a resolution of 256³ using a polytropic equation of state (γ = 1.0001) to approximate an isothermal gas at 10 K. Highly collimated bipolar outflows of momentum $\mathcal {P}=10^{39.52}$ g cm s⁻¹ were impulsively launched $(t_{\rm launch} =0.002\mathcal {T})$ at a constant rate $\mathcal {S}=10^{-67.20}$ cm⁻³ s⁻¹ into an initially uniform environment of density ρ₀ = 10^−19.60 g cm⁻³. The position and orientation of each outflow was chosen randomly. The values for $\mathcal {P}, \mathcal {S}$, and ρ₀ were chosen to approximate star-forming regions like NGC1333. Using the scaling relations above, these define outflow scales of length $\mathcal {L}=0.360\; \rm pc$, mass $\mathcal {M}=17.3\;M_{\odot }$, time $\mathcal {T}=0.366\;\rm Myr$, and velocity $\mathcal {V}=0.961\;{\rm km\;s}^{-1}$. The simulation domain was chosen to be four times the outflow length scale so the corresponding outflow wave number $\mathcal {K}=\frac{2\pi }{\mathcal {L}} = 4\times k_{\rm min}$, where k_min corresponds to the scale of the box. Outflows were continually launched for three times the outflow timescale for a total of 4³ × 3 = 192 outflows launched during the simulation. The turbulence was then allowed to decay for another 3.7 Myr which corresponds to 2.4 dynamical times $(t_d=\frac{L_{\rm box}}{\mathcal {V}})$. In addition to the outflow-driven turbulence case, we also carried forward an isotropically forced turbulence simulation as a control case. In this simulation, turbulent motions were driven by a constant solenoidal forcing function with a spectrum and magnitude chosen to mimic the driving inferred from the jet simulations. In particular, the momentum density injection rate was set to $\mathcal {SP} = \frac{\rho _{_0}\mathcal {V}}{\mathcal {T}}$ and the forcing spectrum had a power law F_k ∝ k^2.5 beginning at the box scale and peaking at the outflow interaction scale $\mathcal {L}$ as seen in Figure 3.

Each outflow was launched with a 5° half opening angle from a region 10 cells across corresponding to an outflow radius of 5800 AU. The velocity, density, and duration of each jet were modeled on a 0.5 M_☉ star ejecting one-sixth of its mass at a velocity of 240 km s⁻¹ giving a v_c = 40 km s⁻¹ and a P_model = M_⋆v_c = 10^39.6 cm s⁻¹. The mass-loss rate was chosen to be 10⁻⁴ M_☉ yr⁻¹ giving a jet density of 0.21ρ₀ and a duration of 0.83 kyr. The actual momentum injected by each outflow was then measured and found to be about 10^39.52 cm s⁻¹ or 83% of P_model corresponding to v_c = 33.3 km s⁻¹. In addition, the outflows injected about twice the mass, however, the total mass injected by all 192 outflows amounts to less than 3% of the ambient material and does not have a significant effect on the overall dynamics. The values for $\mathcal {P}, \mathcal {S}$, and ρ₀ quoted above are the actual values from the simulation.

3.1. Creation and Maintenance of Turbulence

Figure 1 shows density cross cuts from both the outflow-driven and isotropically driven control simulations. Examination of the outflow-driven simulation shows that initially outflow cavities are able to expand into the quiescent environment without being disrupted and easily grow to pc size structures. Since the outflows are bipolar, the total vector momentum injected by each outflow is zero, however, the total scalar momentum defined as ∫|P|dV grows steadily until the cavity interactions begin to dissipate momentum in the expanding shells. By $1.8\mathcal {T}$ the cavities have entirely filled the domain and by $3\mathcal {T}$ the system has reached a statistically quasi-steady state. Comparison of the late time outflow-driven simulation image and the isotropically forced simulation image show that both have reached states of highly disordered flows with structure present on a variety of scales. As we show below the visual impression of turbulence in both cases is supported by statistical measures of the flows as well as their time evolution. The comparison by visual inspection however is noteworthy because of the "holes" which appear in the outflow-driven density distribution. These are created by fossil outflow cavities and exist until the cavities are subsumed by the turbulent motions (Cunningham et al. 2009) or interact with another cavity. Such shells have been observed in turbulent flows around star-forming regions such as NGC1333 (Quillen et al. 2005) and point to an important morphological, rather than statistical, signature of outflow-driven turbulence.

Zoom In Zoom Out Reset image size

In Figure 2, we show a plot (upper left panel) of the average scalar momentum density versus time for the jet simulation. Initially the scalar momentum density injection rate is $\mathcal {P}\mathcal {S}=\frac{\rho _{_0}\mathcal {V}}{\mathcal {T}}$ since the jets are not interacting. Around $0.5\mathcal {T}$ the expanding cavities begin to interact and scalar momentum begins to be dissipated. By $2\mathcal {T}$ the momentum dissipation rate balances the momentum injection rate and the system reaches a steady state. The average mass weighted velocity in our simulation is very close to $\mathcal {V}$ corresponding to Mach 4.8. At $3\mathcal {T}$ the jets are turned off, and the scalar momentum density dissipation rate is shown to be approximately $\frac{\rho _{_0}\mathcal {V}}{\mathcal {T}}$. Figure 2 also shows the development of the average scalar momentum density for the isotropically forced simulation. It is interesting to note that while both simulations have the same initial momentum injection rate by design, the final average specific momentum in both cases is very close to $\mathcal {V}$ in spite of the varied means of forcing.

Zoom In Zoom Out Reset image size

Figure 2 also shows the time development of the two-dimensional probability density functions for the gas density and Mach number. The jets have stirred up all of the ambient material by ∼ $2\mathcal {T}$ and have reached a statistically steady state by ∼ $3\mathcal {T}$. The steady-state density and velocity probability density functions (pdfs) are fairly log normal apart from the spikes associated with low-density, high-velocity material from successive outflow cavities. Once the jets are turned off at $3\mathcal {T}$ the turbulence decays fairly rapidly. Note that by $4\mathcal {T}$ the pdfs show no evidence of the expanding outflow cavities and are characteristic of a purely turbulent isotropic media. The development of log-normal density distributions in isothermal conditions is one measure of a turbulent flow (Vazquez-Semadeni 1994). The four panels in Figure 2 demonstrate one of the key conclusions of our study: transient outflow cavities can set the bulk of initially quiescent material into random but statistically steady supersonic motions.

3.2. Outflow-Driven Turbulent Spectra

Supersonic turbulence in general involves a cascade of kinetic energy (or scalar momentum) from some forcing scale to smaller scales where energy is eventually dissipated in the form of heat and radiated away. In between the forcing scale and the dissipation scale lies the inertial range where the turbulence is scale free and the energy spectra follows a power law $\mathcal {E}\left(k\right)\propto k^{-\beta }$. Recent simulations of isothermal supersonic turbulence by Kritsuk et al. (2007) have found that β ∼ 2 in the inertial range, and that the lack of eddies at the numerical dissipation scale (Δx) causes the spectrum to steepen for scales ⩽8Δx. This lack of eddies inhibits the cascade of energy from slightly larger scales creating a bottleneck that causes the spectra to flatten above 8Δx up to approximately 32Δx. Above 32Δx the cascade is effectively inertial. Modeling the cascade of energy in realistic astrophysical conditions is therefore a difficult task as the numerical dissipation scale is often several orders of magnitude above the physical dissipation scale in high Reynolds number flows characteristic of astrophysical turbulence.

Here, we define the one-dimensional power spectra $\mathcal {E}$ of the velocity u

$$\begin{eqnarray} &&\mathcal {E}\left(k\right)\equiv \frac{1}{V}\int {|\tilde{\mathbf {u}}\left(\mathbf {k}\right)|^2\delta \left(\left|\mathbf {k}\right|-k\right)d\mathbf {k}}S, \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&\mbox{ where }\tilde{\mathbf {u}}\left(\mathbf {k}\right)=\int _V{\mathbf {u}\left(\mathbf {x}\right)e^{-2\pi i\mathbf {k}\cdot \mathbf {x}}d\mathbf {x}} \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &&\mbox{ so that } \int {\mathcal {E}\left(k\right)dk}=\frac{1}{V}\int _V{\mathbf {u\left(x\right)}^2d\mathbf {x}}=\langle\mathbf {u}^2\rangle.\end{eqnarray} \tag{ 5 }$$

In Figure 3, we plot the one-dimensional velocity power spectra of the outflow-driven turbulence as well as the isotropically forced turbulence. The spectrum of the isotropically forced turbulence follows a power law down to scales ∼8Δx (k/k_min = 32) before the numerical dissipation overwhelms the bottleneck it produces. We note of course that resolution is a critical issue in interpreting turbulent power spectra and we will address these issues in a subsequent paper (J. J. Carroll et al. 2009, in preparation). For the present purposes the results allow us to compare the discrete source, outflow-driven turbulence to the continuous, isotropically forced turbulence at the same resolution and Mach number even though our inertial range is only marginally resolved.

Zoom In Zoom Out Reset image size

Consideration of the right-hand panel of Figure 3 demonstrates that the spectrum of the outflow-driven turbulence is quite different from the isotropically forced turbulence. At large scales, the outflow-driven spectrum rises slowly from the box scale to around $\mathcal {K}$ in a broad maximum before turning over and steepening to β = 2.74 at ∼ $3\mathcal {K}$. Nakamura & Li (2007) also found a steeper power law (β = 2.5) in their simulations of outflow-driven turbulence although the presence of self-gravity and the lower resolution made it difficult to isolate the effect of the anisotropic outflow forcing. The isotropically forced spectrum shows a fairly flat distribution over the forcing range and then falls with a shallower power law (β = 2.45) more appropriate to a Burgers model for compressible turbulence in which the random distribution of shocks gives a spectrum with β = 2 without regard to the cascade of energy as in Kolmogorov turbulence. Matzner (2007) suggested that a cascade of momentum in supersonic turbulence would also produce a β = 2. The β = 2.45 is still somewhat steeper than is to be expected for isotropically forced supersonic turbulence. This may be due to the nature of the solenoidal forcing function as well as numerical dissipation. In order to minimize the effect of numerical dissipation we have therefore fit both the spectrum of the outflow-driven turbulence and the spectrum of the isotropically driven turbulence over the same range in wavenumber. Since the fit is closer to the knee in the outflow-driven turbulence spectra, we expect higher resolution simulations to only enhance the difference in β.

The reason for the difference between the outflow-driven turbulent spectra and the spectra expected to be produced by shock-filled turbulence as well as a momentum cascade model, lies in the nature of the driving. Since the shocks in the outflow-driven turbulence are not randomly distributed (i.e., opposing bow shocks associated with bipolar outflows), the turbulence though shock dominated, should not be expected to follow a Burgers model power law. In addition, because the turbulent turnover timescale at the outflow scale $\frac{\mathcal {L}}{\mathcal {V}}=\mathcal {T}$ is equivalent to the timescale for an outflow to clear out a region of volume $\mathcal {L}^3$, the spectrum at scales below the outflow scale $k > \mathcal {K}$ would not be expected to follow a momentum type cascade model either. The steep spectrum at suboutflow scales (with β ∼ 3 corresponding to a velocity length scaling v(l) ∝ l) is consistent with the Hubble type flow seen in the expanding cavities. The time evolution of the spectra (for many outflows) shown in the right-hand panel of Figure 3 provides some insight. The steep slope of the spectrum at suboutflow scales is present already by $0.70\mathcal {T}$ before the outflows have begun to significantly interact. This is consistent with the interpretation that the steep slope of the spectrum at later times $3\mathcal {T}$ arises from the volume swept-up by the multiple fossil outflow cavities rather than their interaction. In other words, we can interpret the steep spectra as due to expanding outflows that sweep up small-scale (high-k) eddies which effectively removes them from the flow. This effect would reduce the observed power on these scales and steepen the slope of the energy spectrum.

3.2.1. Large-Scale Spectra and Forcing Functions

In order to clarify the behavior of outflow-driven turbulence, a single outflow simulation was run. In the right-hand panel of Figure 3, we present the development of the velocity power spectra for this single jet. The six snapshots are equally spaced in log time. For each line the velocity power spectra peaks at approximately the outflow length scale which grows in time. Note the slight secondary plateau at the outflow width scale. As the expanding cavity grows, the total kinetic energy decreases, and the energy shifts to larger scales. For an expanding outflow in a quiescent medium there is, of course, nothing to limit the growth. Outflows in a turbulent medium however, will loose coherence when the momentum in the outflow is comparable to the turbulent momentum swept up by the expanding cavity. This statement is equivalent to the fact that once the outflows reach the scale $\mathcal {L}$ they interact and are destroyed, converting their directed motion up to random motion. For spherical outflows, this will happen on scales ∼$\mathcal {L}$. Collimated outflows would in theory be able to grow to larger size structures before loosing coherence, however, the momentum in the expanding sides of the cavities can be comparable with that in the head of the jet. In fact, we find that after 1, 2, and 3$\mathcal {T}$, the transverse momentum accounts for 30%, 40%, and 50% of the total momentum, respectively. So a driving scale ⩾ $\mathcal {L}$ is still a fairly good approximation. One would expect the spectrum to fall off steeply for scales above the driving scale $k < \mathcal {K}$. Both panels in Figure 3, however, show that multiple outflows inject their energy over a broad range of scales, and while each outflow injects most of its energy at $L_{\rm of} \sim \mathcal {L}(t)$, the shorter turnover times for scales closer to $\mathcal {L}$ cause the spectra to flatten making the identification of a single outflow-driving scale difficult. To address this issue we consider combining a power-law forcing function $\mathcal {F} \propto k^\alpha$ with a power-law energy spectrum $\mathcal {E} \propto k^{-\beta }$ and a dissipation rate to smaller scales $\Pi \propto \frac{v_l^2}{t_l} \propto kv_k^3 \propto k^{5/2}\mathcal {E}_k^{3/2} \propto k^{3/2(5/3-\beta)}$. With these relations we can look for the steady state solution for the energy at each scale:

$$\begin{eqnarray} &&\frac{d\mathcal {E}_k}{dt}=\mathcal {F}_k-\frac{d}{dk}\Pi _k=0. \end{eqnarray} \tag{ 4 }$$

Using the above assumptions this equation implies α = 3/2(1 − β). In other words, for α < 3/2 a forcing function with a positive slope will result in an energy spectrum that would have a negative slope in the forcing region. This is likely the case in our simulations for scales close to, but shorter than, our nominal driving scale. At the largest scales we found that the outflow turbulent spectra had a positive slope with β = −0.57. Here, a relatively steep forcing function $\mathcal {F}\propto k^{2.4}$ could peak at ∼ $\mathcal {L}$ but reproduce the observed energy spectrum at $L > \mathcal {L}$. Thus, while our simulations do not begin with a specified forcing function and simply allow the interactions of outflow cavities to self-consistently establish a turbulent flow we can still draw broad conclusions about the spectral nature of the imposed forcing.

There have been many simulations of turbulence in both the atomic and molecular ISM (Stone et al. 1998; Kritsuk et al. 2007; Porter et al. 1992; Mac Low et al. 1998). Some of these focus on decaying turbulence while others explore driven turbulence. In the majority of these studies the turbulence is forced with a continuous, isotropic, and predetermined spectra. In this paper, we have focused on turbulence in subvolumes of molecular clouds and considered the more realistic case of discrete, nonisotropic forcing in physical space rather than Fourier space. Our goal has been to address the ability of transient protostellar outflows (fossil cavities) to drive their surrounding media into supersonic turbulent motions. Our results show conclusively that for reasonable values of outflow momenta $\mathcal {P}$ and outflow rate per unit volume $\mathcal {S}$, outflows generate supersonic turbulence. We have found that the spectrum and physical properties of outflow-driven supersonic turbulence differ from that produced by isotropic random forcing. We note first that we are able to broadly identify an outflow forcing scale $\mathcal {L}$ which is of the order of magnitude comparable with that predicted by Matzner (2007). While the density pdfs are log normal consistent with a isothermal turbulent flow (Vazquez-Semadeni 1994), the spectra of outflow-driven turbulence on scales smaller than $\mathcal {L}$ is steeper than that expected from the turbulent cascade model as in isotropically forced turbulence. We interpret this to be the result of overlapping outflow cavities preventing smaller scale structures from forming thereby facilitating a sort of bottleneck. The small-scale spectrum is therefore dominated by the spectrum of outflow cavities themselves which follow a power law consistent with a Hubble-type flow. The departure from β ∼ 2 velocity power spectra has been seen in other models with discrete forcing in the form of supernova-driven turbulence in the ISM (Joung & Mac Low 2006). In those models, a broad shallow plateau associated at large scales was found which then turned over to a steep spectrum at smaller scales which was attributed to numerical dissipation.

Our results add weight to the argument that protostellar outflows can play a significant role in altering their star-forming environments (Nakamura & Li 2007). One must be careful not to overinterpret the current simulations which do not include a number of important physical processes such as magnetic fields and self-gravity. Our goal in this initial study was to focus on one key aspect of the physics at work in low mass star-forming environments (multiple transient protostellar outflows) in order explicate its mechanisms and energetics. The role of self-gravity, in particular, is likely to be important in determining the final turbulent state of the gas (Field et al. 2008) and magnetic fields may be able to transmit energy from regions of protostellar outflow to broader regions of the cloud (Kudoh & Basu 2006). In spite of these important caveats, our results clearly establish a baseline for understanding how outflows can strongly effect stellar birth in clusters which are, after all, the regions where stars actually form.

We thank Chris Matzner, Chris McKee, and Mordecai-Mark Mac Low for extremely useful discussions as well as the referee for their insightful criticisms and helpful suggestions. Hector Arce, John Bally, Pat Hartigan, and Tom Ray were also generous with their time. Tim Dennis, Kris Yirak, Brandon Schroyer, and Mike Laski provided invaluable support and help. Support for this work was in part provided by NASA through awards issued by JPL/Caltech through Spitzer program 20269, the National Science Foundation through grants AST-0406823, AST-0507519, and PHY-0552695 as well as the Space Telescope Science Institute through grants HST-AR-10972, HST-AR-11250, and HST-AR-11252. We also thank the University of Rochester Laboratory for Laser Energetics and funds received through the DOE Cooperative Agreement DE-FC03-02NA00057.