PROTOSTELLAR OUTFLOW EVOLUTION IN TURBULENT ENVIRONMENTS

Both turbulence and outflows appear to be ubiquitous phenomena associated with star formation. 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). Though the issue of whether the observed spectra are a consequence of gravitational fragmentation or mechanical forcing remains a topic of current debate (Field et al. 2008). Outflows driven by newly formed stars also appear uniformly across star-forming environments. Outflows appear to form very early in the formation of a star (the Class 0 phase) and continue until a sizable fraction of the stellar mass has been assembled. In addition, outflows occur across the entire spectrum of stellar masses from brown dwarfs (Cabrit 2007; Whelan et al. 2005) all the way to high mass O and B stars (Shepherd 2003).

Given the apparent ubiquity of both turbulence and outflows, one can ask about their respective importance in terms of impact on the star formation process. In recent years, the role of turbulence in moderating star formation has become an active area of research (Mac Low & Klessen 2004; Krumholz & McKee 2005). Some workers have argued that turbulent pressure in the form of a space-filling isotropic distribution of eddies provides support for molecular clouds against gravitational collapse (Goldreich & Kwan 1974; Fleck 1980). In this way, turbulence is thought to be the agent which establishes equilibrium and provides a means for molecular clouds to be long lived structures. One difficultly with this view, however, is the recognition that both hydrodynamic and magnetohydrodynamic turbulence decay on a timescale of order the cloud free-fall time (Stone et al. 1998; Mac Low 1999). Thus if clouds are stable, long-lived structures some means must be found to support them against self-gravitational collapse (Blitz et al. 2007).

When discussing issues of jet feedback one can distinguish between three scale lengths. "Microscale" feedback concerns the effect of outflows on their own launch scale (L < 10 AU). "Mesoscale" feedback relates to the impact of jets on environment associated with infalling envelopes (L < 10⁴ AU). Considerable work has focused on these domains, i.e., the way outflows can sculpt envelopes (Delamarter et al. 2000; Cunningham et al. 2005) or their influence on the accretion onto the star (Krumholz et al. 2005). Finally, "macroscale" feedback, which is the concern of this paper, is associated with the effect of outflows on the scale of star-forming clusters or the entire molecular cloud itself. The role of outflows on large scales became of interest after it was discovered that some outflows extended to parsec scales which is of order cloud size or at least the size of the"clump" from which star clusters form within a cloud (Bally & Devine 1994; Bally et al. 1996; Stanke et al. 1999; McGroarty et al. 2004). Young star clusters such as NGC 1333 are the domains of interest in this problem. The impact of the mechanical feedback driven by the formation of multiple stars on their cores and parent cloud remains a subject of intense debate, and this is precisely the location where issues in outflows and turbulence overlap.

Numerous observational studies of star-forming regions have shown that, in general, the energy in active outflows is large enough to support a cluster or cloud against collapse if that energy can be coupled to cloud turbulence (Bally et al. 1996; Bally & Reipurth 2001; Knee & Sandell 2000; Quillen et al. 2005; Warin et al. 1996). These investigations support the idea that outflows can provide a source of internal forcing to drive turbulence. Theoretical treatments of feedback from multiple spherical outflows creating a self-regulating star-forming system was first explored by Norman & Silk (1980). Analytical work by Matzner (Matzner & McKee 2000; Matzner 2001, 2007) have 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.

The first simulation-based study of the problem was presented by Mac Low (2000) whose results indicated that outflows could drive turbulent motions. Decay rates for outflow powered turbulence in this study were ∼4× faster than that of uniformly driven turbulence. These early simulations were of relatively low resolution and further high-resolution work is needed. Recent works by Li & Nakamura (2006) and Nakamura & Li (2007) have mapped out the complex interplay between star formation outflows and turbulence. Their simulations include a self-consistent formulation of driving outflows from newly formed stars, and they concluded that protostellar outflows were a viable means of generating turbulence in star-forming clusters. However, studies of single jets by Banerjee et al. (2007) came to the opposite conclusion. Exploring the volume filling averages of postshock material, they concluded that single jets would not leave enough supersonic material in their wakes to act as a relevant source of internal forcing.

Many studies of protostellar turbulence have focused on active outflows. These may not, however, be the only means of coupling outflows with turbulence. Numerical studies (Cunningham et al. 2006a) have shown that direct collisions of active outflows are not effective at transferring jet/outflow momentum to cloud material (by active, we mean that the momentum flux of the jet or wide-angle wind remains roughly constant throughout the interaction). Along similar lines, Quillen et al. (2005) explored the relation between turbulent cloud motions and observed outflow activity in NGC 1333. Their results showed that the location of active outflows did not effect the degree of turbulence in the cloud. Quillen et al. (2005) were, however, able to identify many subparsec sized, slowly expanding cavities. These they identified as dense wind-swept shells of gas that expand into the cloud after the driving source of the outflow has expired. The net mechanical energy required to open all these cavities represents a significant fraction of that required to drive the cloud turbulence. The authors argued that these cavities were relics of previously active molecular outflows. Because cavities have a longer lifetime than those opened by active and younger outflows, they would be more numerous than higher velocity, younger outflows. These cavities could provide the coupling between outflows and turbulence in molecular clouds. Thus one means of re-energizing turbulence in clouds may be through"fossil outflow" cavities which are disrupted after they have slowed to speeds comparable to the turbulent velocity of the ambient cloud. In Cunningham et al. (2006b), high-resolution Adaptive Mesh Refinement (AMR) simulations of transiently driven outflows were carried out. These simulations, which included a treatment of molecular chemistry and cooling, tracked the development of fossil outflow cavities driven by transient outflows out to near parsec scales. The simulations provided observational predictions in terms of position–velocity diagrams that were consistent with the momentum scaling relations used in Quillen et al. (2005).

While the potential for protostars to act as a means of internal forcing remains attractive, many outstanding questions remain open. The energetics of turbulent resupply must be understood in terms of a coupling coefficient between outflow and cloud. The nature of outflow interactions must also be specified as Cunningham et al. (2006a) have already shown that active outflows are not efficient in setting large volumes of gas into random motions. Thus, the emphasis shifts to fossil cavities and their interactions with each other and/or with pre-existing turbulent medium. We expect that turbulence can be generated by collisions between fossil outflows because over the lifetime of a cluster the cavities will tend to fill the cluster volume Cunningham et al. (2006a). Individual fossil cavities need not experience collisions but will still couple to the ambient medium. The mechanisms through which the momentum and energy of collimated outflow structures is "isotropized" has yet to be articulated. Something of a chicken and an egg problem exists in this regard. In Quillen et al. (2005) expanding outflow cavities where observed down to speeds a few times that of the turbulence v_outflow ∼ 2v_turb ∼ 1 km s⁻¹. They argued that the cavities would be subsumed by turbulent eddies once they decelerated to speeds of order v_turb. The expectation that the directed momentum and energy of a single decelerating cavity can be given up to turbulence must be explored. This is the question we seek to address in this paper.

In the current work, we perform a series of a simulations meant to take the work of Cunningham et al. (2006b) to the next level of complexity. This paper is meant as an critical intermediary step to full-scale simulations of multiple interacting transient outflows. Here we consider the role of individual transient jets which are ejected into a decaying turbulent medium. We note that the evolution of a radiative jet in a turbulent medium has yet to be explored, and our simulations are relevant to YSO jet morphology and the re-energization of turbulence via jets and the disruption of the outflow cavities. Finally, our simulations are relevant to issues addressed in the work of Banerjee et al. (2007) and their conclusions that jets cannot drive supersonic turbulence. We note that the question of the steady state rate of decay of such flows remains outside of the scope of the present work as this question could only be addressed by models that include a steady injection of multiple outflows which are evolved for longer timescales than the numerical models considered in this paper. This issue will be addressed in a forthcoming paper (Carroll et al. 2009).

We have carried forward a series of numerical simulations to investigate the interaction of young stellar jets with a turbulent environment. Each of the simulations is carried forward on a periodic domain spanning L_x × L_y × L_z = 4000 ×  2000 ×  2000 AU discretized to a 504 × 252 × 252 computational grid. An environment of density ρ_a = 50 cm⁻³ with temperature T_a = 200 K is set into a state of decaying turbulence from an unspecefied source. Possible candidates for supplying background turbulence within star-forming environments include an initial turbulent state present during the formation of the cloud from compressive motions in a turbulent interstellar medium (Elmegreen 2000), differential rotation of galactic disks, gravitational collapse, the ionizing radiation, winds and supernovae from massive starts, and protostellar outflows (Mac Low 1999). The spectrum of initial motions can be written as

$$\begin{eqnarray} v_x({\bf x}) =& \sum _{i,j,k} A_{x i,j,k} \sin ({\bf k}_{i,j,k} \cdot {\bf x} + \phi _{x i,j,k}) \nonumber \\ v_y({\bf x}) =& \sum _{i,j,k} A_{y i,j,k} \sin ({\bf k}_{i,j,k} \cdot {\bf y} + \phi _{y i,j,k}) \\ v_z({\bf x}) =& \sum _{i,j,k} A_{z i,j,k} \sin ({\bf k}_{i,j,k} \cdot {\bf z} + \phi _{z i,j,k}),\nonumber \end{eqnarray} \tag{ 1 }$$

with wavenumbers

$$\begin{equation} {\bf k}_{i,j,k} = \frac{2\pi i}{2000\;\hbox{AU}} {\bf \hat{x}} + \frac{2\pi j}{2000\;\hbox{AU}} {\bf \hat{y}} + \frac{2\pi k}{2000\;\hbox{AU}} {\bf \hat{z}}. \end{equation} \tag{ 2 }$$

The summation is taken over all i, j, and k between 0 and 31, excluding the i = j = k = 0 term. The longest wave mode seeded in each spatial direction is therefore equal to the length of the short edges of the computational grid, and the shortest wave mode in each spatial direction is resolved across approximately 16 computational zones. The phase angles ϕ_*,i,j,k are chosen at random between 0 and 2π. The direction of each wave mode is determined by choosing two of the amplitude components, A_xi,j,k and A_yi,j,k at random between −1/2 and 1/2. The third amplitude A_zi,j,k is computed to satisfy the constraint that the initial velocity field be solenoidal. Each velocity amplitude A_i,j,k is normalized so that the initial velocity power spectrum follows a power-law dependence on wavenumber E_Ko(|k|) ∼ |k|^−β with net kinetic energy prescribed by the desired root-mean-squared (rms) turbulent speed (v_turb) as $V \int E_{Ko}({\bf k}) d{\bf k} = \frac{1}{2} \rho _a v_{\rm turb}$, where V is the volume of the computational domain, E_K(k) is the power spectrum of velocity perturbations,

$$\begin{equation} E_{K}({\bf k}) = \left| \mathcal {F}\left(v_x({\bf x}) \right) \right|^2 + | \mathcal {F}(v_y({\bf x}))|^2 + \left| \mathcal {F}\left(v_z({\bf x}) \right) \right|^2, \end{equation} \tag{ 3 }$$

and $\mathcal {F}$ is the Fourier transform operator, $\mathcal {F}_{{\bf k}} = \int f_{{\bf x}} e^{i {\bf k} \dot{\bf x}} d{\bf x}$. It is also convenient to define the one-dimensional power spectrum as

$$\begin{equation} E_{K}(k) = \int E_{K}({\bf k}) \delta (|{\bf k}|-k) d{\bf k}. \end{equation} \tag{ 4 }$$

The simulations which begin with an initially turbulent flow field are all seeded with an initial rms turbulent velocity v_turb = 10 km s⁻¹, yielding an initial rms turbulent Mach number of 6.03, and velocity power spectrum index of β = −2, characteristic of highly supersonic Burger's-type turbulence (Porter et al. 1992; Mac Low 2003).

Our simulations employ periodic boundary conditions at each boundary interface. This choice of boundary conditions prevents boundary artifacts from polluting the spectrum of motions in the grid. Jet-driven cavities are therefore launched from a grid-embedded region rather than from grid boundaries. This choice of domain setup allows the simultaneous use of periodic grid edges and jet launching. The jet-launch region extends from 125 AU < x < 250 AU, and (y − 1000 AU)² + (z − 1000 AU)² < r²_j, where r_j = 250 AU. The launch region maintains a constant density ρ_j = 250 cm⁻³ and temperature T_j = 200 K profile and enforces a time-dependent jet velocity as

$$\begin{equation} v_x(t) = \left\lbrace \begin{array}{@{}l@{\qquad}l@{}} v_o \exp \left(-t/t_{\rm off} \right) & t<2 t_{\rm off} \\ 0 & {\rm otherwise}, \\ \end{array} \right. \end{equation} \tag{ 5 }$$

with initial velocity v_o = 150 km s⁻¹. The jet velocity is smoothed following a quadratic form so that the outer radius of the launch region has a velocity that is 90% of the velocity at the center of the jet.

The choice of physical parameters outlined above encompasses a compromise between the desire for computational efficacy and the most observationally relevant molecular cloud conditions. Giant molecular cloud complexes exhibit characteristic densities of 100–10⁴ cm⁻³, with line widths of 1–10 km s⁻¹, and characteristic temperatures in the range of 10–30 K (Mac Low & Klessen 2004). Our parameters lie toward the lower end of the density range, the higher end of the observed turbulent speeds, and are somewhat hotter. While the AstroBEAR code used to perform the calculations includes a module for nonequilibrium H₂ chemistry and the collisional ionization dynamics for H and He (Cunningham et al. 2005), we have chosen to reduce the computational expense of the calculation by modeling the system as a gas with molecular weight of μ = 1m_H, with a simplified cooling law that depends only on temperature. The initially T_a = 200 K, μ = 1m_H gas has an equivalent thermal pressure to that of a molecular gas with coronal abundance of He and T_a = 86 K, μ = 2.33m_H, still somewhat higher than the central regions of typical molecular clouds though perhaps comparable to cloud edges that are more exposed to background X-ray and cosmic-ray heating sources. The somewhat large choice of initial turbulent speed compensates for the the somewhat high initial temperature by retaining the correct turbulent Mach number. All else fixed, lower densities and higher ambient temperatures reduce the rate of cooling behind shock fronts which, in turn, reduces the number of iterations required to evolve the stiff cooling source term with adequate numerical precision. Our choice of parameters is nonetheless characteristic of strong cooling and does lead to the formation of a thin Mach disk at the head of the simulated outflow (relative to the grid resolution) and a collimated jet bow shock structure. In addition to requiring more iterations to evolve the cooling terms, simulations with more realistic physical parameters require greater resolution to capture comparable resolution of the cooling behind strong shocks with the number of grid zones scaling roughly as [δnδT⁻¹]³ and, accounting for the smaller numerically stable time step, the computational effort scales as [δnδT⁻¹]⁴.

A total of five simulations have been performed using the AstroBEAR astrophysical fluid code to carry forward the integration of the Euler equations using the MUSCL–Hancock shock capturing scheme and Marquina flux function as described by Cunningham et al. (2009) with an energy sink term which models the effects of optically thin radiative energy loss as described by Cunningham et al. (2006a). These simulations did not use the AMR or magnetic field capacities of AstroBEAR. The suite of simulations includes three models in which jets with varying decay parameters are driven into a decaying turbulent environment (run 1–2), one model which follows the decay of the turbulent medium in the absence of disruption by a jet (run 0), and one model in which a nondecaying jet is driven into a quiescent environment (run J). In Table 1, we summarize the initial conditions and results of each simulation. The fifth column from the left tabulates the cooling efficiency as measured by the net energy loss at the end of the simulations, the sixth column tabulates the total energy injected into the grid at the jet-launch region over the duration of the simulations relative to the initial turbulent energy seeded in the grid, and the rightmost column tabulates the least-squared fit to the velocity spectrum power-law index at the end of the simulation $\beta |_{t=t_{f}}$ for 10^1.5 AU < k < 10^2.3 AU. In our simulations, the power-law spectral dependence expected for inviscid turbulent flow begins to break for larger k due to numerical dissipation and for smaller k due to the inability to represent modes with characteristic length greater than the shorter edges of the simulation domain.

Table 1. Simulation Summary

	toff	Turbulence	Final Time tf	Cooling Efficiency	$\frac{E_{K,{\rm jet}}}{E_{Ko}}$	$\beta |_{t=t_f}$
Run 0	No jet	✓	922 yr	36.1%	⋅⋅⋅	−2.09
Run 1	8 yr	✓	922 yr	49.0%	0.285	−2.24
Run 2	100 yr	✓	221 yr	9.45%	3.56	−2.38
Run 3	Nondecaying jet	✓	161 yr	22.7%	71.0	−2.43
Run J	Nondecaying jet	⋅⋅⋅	161 yr	22.1%	⋅⋅⋅	−2.46

Download table as:  ASCIITypeset image

The simulations presented in this work are of considerably shorter timescales than that of actual protostellar outflows. This is primarily due to computational constraints imposed by the need to simultaneously resolve the overall outflow, the outflow injection radius R_j, and turbulent eddies with scale length <R_j. The simulations of extinct outflow sources presented here are interpreted as models of the fossil relics of extinct protostellar outflow activity in molecular clouds. The simulated fossil cavities are intended to be roughly comparable to actual protostellar outflow cavities at comparable time relative to the outflow shut-off time, t/t_off.

The right panel of Figure 1 shows a crosscut about the mid-plane of the simulation domain at the end (t = t_f) of each simulation. The initial and jet-launch conditions of each simulation are summarized in the left three columns of Table 1. The right three columns of Table 1 list diagnostic parameters at the end of each simulation. The "cooling efficiency," which is taken as the fraction of the total energy budget lost from the simulation domain via radiative cooling, is listed for each run. This quantity measures the radiative loss of energy from the system that may have otherwise been available to sustain turbulent motion. Note that this cooling efficiency increases for simulations that run for longer physical time. This is because these simulations have more time for radiative energy loss to occur. It also increases for simulations with longer jet decay times as these jets feed energy into their bow shocks for longer times and have stronger cooling compared with jets that turn off earlier.

Zoom In Zoom Out Reset image size

The propagation speed of the Mach disk structure at the head of the simulated outflows can be estimated analytically for the case of the continuously driven jet models (run 3 and run J). Blondin et al. (1990) consider the momentum flux across the inward and outward facing surfaces of the Mach disk structure at the head of jet-driven flows to derive an analytic estimate for jet bow shock propagation speeds,

$$\begin{equation} v_{\rm bs}=v_j \left[1+\left(\frac{\rho _j}{\rho _a}\right)^{-1/2}\right]^{-1}. \end{equation} \tag{ 6 }$$

This expression evaluates to v_bs = 104 km s⁻¹ for the continuously driven jet model parameters considered in this paper. The shocks delineating the Mach disk at the head of the continuously driven jets (run 3 and run J) propagate ∼3500 AU over t_f = 161 yr, corresponding to propagation speeds of 103 km s⁻¹, in excellent agreement with the analytic estimate. Comparison of the two simulations shows that the turbulence does not effect the average propagation speed of driven jets.

Figure 2 shows the time evolution of the total mechanical energy in the grid for each simulation. Because the mechanical energy of the continuously driven jet quickly overwhelms that of the turbulent environment, the time evolution of the mechanical energy for run 3 and run J is indistinguishable in Figure 2. For decaying jets, run 1 and run 2, the net mechanical energy increases until the driving source begins to decay at t ∼ t_off. Interestingly, for t>t_off, the mechanical energy decays in a roughly similar form as that of the control simulation, run 0, which follows the decay of the turbulent environment alone.

Zoom In Zoom Out Reset image size

3.1. Jet Morphology

The morphological effects of the turbulent medium on the jet-driven bow shock in run 3 are readily apparent despite the comparatively slow turbulent speed in the ambient medium (v_turb is only 10% of the jet propagation speed). This is particularly apparent in the change in behavior of the vortex ring in run J located near x ∼ 2750 AU in comparison to run 3. Specifically, the formation of a clearly defined vortex ring is suppressed for the case of the identical jet conditions that are driven into a turbulent medium in run 3. Furthermore, the entire bow shock from the base of the jet to its head shows significant disruption compared to the control case. This is most apparent, and can be most easily analyzed quantitatively, at the Mach disk at the jet head. In the jet which is driven into a quiescent environment (run J), a thin shock-bounded Mach disk structure forms at the head of the outflow. In the case, where the same jet-launch conditions are imposed into a turbulent environment (run 3), the Mach disk is observed to broaden significantly. Figure 3 presents a comparison of the Mach disk width (W_M) for both runs. The figure was constructed by measuring the extent of the Mach disk at its widest point using compression in the jet head as a tracer (ρ>750 cm⁻³). The figure shows that the Mach disk width is increasing in the turbulent simulation by more than a factor approximately 4 over that in the control case.

Zoom In Zoom Out Reset image size

To quantify this change run J shows a Mach disk width of W_M ∼ 250 AU. Dividing by the age of the jet yields an average expansion rate of the Mach disk of v_exp = W_M/t_f = 7.1 km s⁻¹. For run 3, the jet driven into the turbulent environment, the Mach disk structure shows W_M ∼ 800 AU. The average expansion rate is v_exp = W_M/t_f = 22.7 km s⁻¹. Note that the additional expansion of the Mach disk in the turbulent case exceeds the initial ambient turbulent velocity, v_turb = 10 km s⁻¹. Thus, it appears that more than just the entrainment of ambient turbulent motion is at play in this behavior.

The change in the behavior of the Mach disk highlights the fact that even structures which are denser and propagate more quickly than the ambient turbulence are significantly disrupted in its presence. This result suggests an interpretation of jet evolution in a turbulent medium which can account for the eventual disruption of fossil cavities in decaying jets. It is well known that initially laminar shock flows can be disrupted via the action of various instabilities associated with strong cooling, i.e., radiative (Sutherland et al. 2003) and thin shell modes (Vishniac & Ryu 1989; Vishniac 1994). Given this susceptibility to fluid instability, the motion of the ambient turbulence with eddies occurring at a variety of scales should provide a space-filling environment of multimode perturbation seeds. Because the turbulent field is composed of a continuous spectrum of eddies, it will seed the outflow shell with perturbations that correspond to the fastest growing modes of any instability to which the outflow-swept shell may be subject. These seeds act to initiate instability growth in the bow shock delineating the cavity and facilitate the coupling of the energy and momentum of the original outflow to the turbulent environment.

Consider for example a shell of shocked gas of width h with internal sound speed c. The growth rate for modes of the nonlinear thin shell instability for an initial displacement L can be approximated as (Hueckstaedt et al. 2006):

$$\begin{equation} \Gamma \sim C_d^{-1/2} c k\sqrt{k L}, \end{equation} \tag{ 7 }$$

where C_d is of order 1. The growth of the instability requires the displacement L be of order the shell width, L ∼ h. Since the shells are thin, h ≪ r_j, and the timescale for the bow shock to fragment will be short compared to a characteristic hydrodynamic timescale of the flow, which in this case we could take to be the crossing time of the jet across the cluster scale t_cl ∼ r_cl/v_j. Thus, considering only the nonlinear thin shell instability with L ∼ h and k ∼ 1/h, we see that

$$\begin{equation} \Gamma t_{\rm cl} \sim \frac{c}{v_j} \frac{r_{\rm cl}}{h} \sim M^{-1} \frac{r_{\rm cl}}{h}. \end{equation} \tag{ 8 }$$

Given typical Mach numbers in YSO jets of M ∼ 100 and r_cl ∼ 1, we have Γt_cl ≫ 1. Thus, perturbation seeds provided by the turbulence will lead to the fragmentation of the bow shock by the time the cavity has reached its maximum extent.

Putting the issue of the generation of turbulence aside for a moment, we see that issues associated with the propagation of active jets are raised when we consider the disruption of the bow shock. Figure 4 shows a three-dimensional visualization of run 3. Note the fragmentation of the outflow into clumpy filaments. In particular, note the fragmentation of the Mach disk associated with the jet bow shock complex into clumps. Each clump drives a smaller scale bow shock in its wake, resulting in the formation of filamentary density enhancements in the Mach disk. Similar irregular structures, including clumps and filaments in the Mach disk, have been observed in a number of jet systems including HH 1/2 and HH 47 as has been discussed by Bally et al. (2002) and Hartigan et al. (2005), respectively. Yirak et al. (2007) have noted similar bow shock structures when the head of a propagating jet is impeded by stationary density homogeneities embedded within the ambient environment. Thus, it appears that filamentary fragmentation of protostellar jet bow shocks may be characteristic of an inhomogeneous environment in general, regardless of the form of such inhomogeneities.

Zoom In Zoom Out Reset image size

Observations by Bally & Reipurth (2003) and Bally et al. (2006) also reveal shock delineated knotty deflected structures within protostellar jet beams. Such deflections are puzzling when they appear discontinuous, with the jet showing a sharp deflection only some ways down the beam. While Yirak et al. (2007) find that sufficiently dense inhomogeneities can cause similar jet deflections, the motions of the turbulent environment in the present simulations show no comparable effect on the continuously (run 3) or long-driven (run 2) jets. However, the turbulent motions of the ambient environment do sufficiently disrupt the lateral edges of the jet bow shock structure to expose the driving jet beam to the influence of nearby turbulent eddies. It is interesting to note that the turbulent eddy scales, which are represented in the computational domain (L ⩽ 4000 AU), do act to compress the edges of the long decay time jet beam (run 2). Thus, it is possible that supra-jet-scale flow eddies (L ≫ bow shock radius), which would have higher velocity, may drive bending and global disruption of jet structures. On the issue of large-scale disruption, we note the work of Lebedev et al. (2004) in which the deflection of supersonic jets via interaction with a large-scale transverse flow was considered. This study explored the issue both experimentally and numerically, and addressed the issue of a jet interacting with side-wind after propagating some distance from the jet source. These studies raise the possibility of an encounter with a large eddy as the cause of discontinuous jet bending seen in observations. However, the effect of very large-scale eddies is not considered in the present study, as their inclusion would require a much larger computational volume and many more grid points to simultaneously achieve adequate resolution of the jet and the large-scale flow.

3.2. Turbulence and Kinetic Energy Power Spectra

The left column of Figure 1 shows plots of the velocity power spectra for each simulation at several times over the course of the simulations. The power spectra have been constructed by binning the three-dimensional Fourier transform of the computational domain into spherical slices of spectral width Δk = 2π/Δx = 0.792 AU⁻¹.

First we note the degree to which the outflows interact with the background turbulence. Only the jet-driven cavity with the shortest driving decay time in run 1 becomes completely disrupted by, and subsumed into, the background turbulence. However, the more slowly propagating radial edges of the cavity bow shock in the longer jet decay time runs all show some degree of coupling with the ambient turbulence. In the case of the continuously driven jet in the turbulent environment (run 3), approximately one cavity bow shock radius (R_BS) crossing time (t_BS ≈ R_BS/v_turb) has elapsed by the end of the calculation. In that simulation the oldest portion of the outflow limb, which is closest to the outflow source, shows the most disruption qualitatively.

Moving now to the power spectra, we see that the longer lived jet sources of run 2 and run 3 show increasing power in the largest flow eddies (smallest wavenumber). This indicates that the outflow cavities in these simulations provide sufficient power at length scales comparable to the width of the simulation domain to support turbulent motions in the ambient flow against energy decay. Even though the rapidly decaying jet (run 1) injects comparatively less energy, we also find more energy in the largest eddies compared to the case of turbulence without a jet. Much of the mechanical energy injected into the domain by the outflow source appears as large-scale motion (smallest wavenumber) in this simulation as well.

The increasing power at small wavenumber exhibited by the jet models can be interpreted as the net result of the decay of turbulent energy at large wavenumber along with the resupply of turbulent energy by the injected outflow toward smaller wavenumber. The more powerful outflow source in run 3 resupplies energy toward smaller wavenumbers faster than the turbulent cascade of drains energy toward larger wavenumbers. In the case of the outflow cavity models (run 1 and run 2), modes with the smallest wavenumber increase even as the total kinetic energy in the grid decreases (Figure 2). We therefore identify the driving scale of turbulent resupply associated with the scale of cavity propagation as k_drive ∼ 1/L_BS. Furthermore, the simulations that drive outflows exhibit steeper power-law slopes than the decaying turbulence case. At the end of the respective simulations, the decaying turbulence model in run 0 reveals a spectrum with power-law index β = −2.09 compared to the jet-only model with β = −2.46. These cases bracket the slope of the spectra for the outflow with ambient turbulence models, which show a clear steepening with increasing outflow energy (Table 1). This steepening can be interpreted as an outflow-induced suppression of subcavity scale flow eddies. The power spectra therefore indicate that outflow cavities act to support the ambient turbulence in two ways. First the outflows power eddies of comparable extent to the length of the cavity L_BS and second, the smoothly driven outflow sources under consideration in this work drive smooth cavity interiors, which are devoid of eddies of smaller scale than the spatial extent of the cavity, thereby inhibiting the cascade of energy to suboutflow scales.

We now address the issue of cavity disruption. Note that the density distribution of run 3 in the neighborhood of the jet is characteristic of a uniform, laminar flow. After a short time, the energy injected by the jet greatly exceeds that of the ambient turbulence. Therefore, the power spectra at the end of this simulation are dominated by the jet at all scales, and the spectra of run 3 appear indistinguishable from that of the jet-only control in run J. The density distribution crosscuts shown in the right column of Figure 1, however, reveal that the outflow cavity driven by the short, nearly impulsive jet in run 1 has been completely subsumed into turbulent eddies. By the end of the simulation, the flow field for the short jet pulse (run 1) is qualitatively similar to that of the turbulence-only control (run 0). This is the behavior that was expected by Quillen et al. (2005) and Cunningham et al. (2006b) for long extinct outflow structures embedded in turbulence. Furthermore, the net mechanical energy of the disrupted outflow cavity in run 1 decays at a rate comparable to that of the turbulence-only control in run 0. This implies that the disrupted outflow cavity becomes turbulent itself, as it evolves under the influence of the turbulent environment.

The time-dependent evolution of the outflow cavity driven by the short jet pulse is illustrated in Figure 5. The figure is composed of several plots, which show cross cuts about the simulation midplane of the logarithm of the Mach number (left), the logarithm of the gas density in cm⁻³ (center), and the spatial distribution of the fraction of outflow ejected gas as followed by an advected contaminant (right) at t = 92,  230,  461,  691, and 922 yr from top to bottom. The color mapping of the plots of Mach number is constructed to reveal regions of supersonic flow only. The images in the right column showing the advected "jet contaminant" reveal that by the end of the simulation the ejected gas has been widely distributed by the turbulent flow across the computational domain. The end of the simulation at t = 922 yr corresponds to approximately one turbulent crossing time across the shorter edges of the computational domain of L/v_turb ∼ 2000 AU/10 km s⁻¹ = 948 yr. By approximately one bow shock radius crossing time, t ∼ 461 to 691 yr, the initially laminar jet flow has evolved into a turbulent flow indistinguishable from the turbulent motions of the cloud. Furthermore, the flows with the highest Mach number at the end of the simulation correspond to regions containing the highest concentration of jet gas. This indicates that one crossing time after the outflow source expires the transform of bulk, directed fossil cavity motion into a fully turbulent flow is complete.

Zoom In Zoom Out Reset image size

3.3. Protostellar Jets and the Production of Turbulence

The question of if, where, and how protostellar outflows drive turbulence in their environments remains an open issue. In this work, we have attempted to address one facet of the problem. In the recent work of Banerjee et al. (2007), the problem of outflow-driven turbulence was also addressed; however, these authors came to a very different conclusion than we have in this and previous works. In Banerjee et al. (2007), single-outflow cavity structures were used as a means of determining the efficacy of outflows as a source of supersonic turbulence. In that work it was argued that regardless of the momentum imparted to the cloud by an outflow, effective outflow–turbulence coupling will only occur if a significant volume fraction of the outflow cavity retains supersonic flow speed at the time it is disrupted into the cloud. Using numerical simulations, Banerjee et al. (2007) found that most of the volume overrun by either transient or continuously driven isothermal jet models remains subsonic after passing through the bow shock, and thereby conclude that outflows are unlikely to act as a significant source of supersonic turbulence in molecular clouds. The results in our paper point toward the opposite conclusion. In what follows, we discuss the differences between the basic scenario for jet-driven turbulence presented in Banerjee et al. (2007) and the ideas that underpin this work and our previous studies (Cunningham et al. 2006a; Quillen et al. 2005; Cunningham et al. 2006b). In addition, we address specific issues raised by the simulations in Banerjee et al. (2007).

The basic premise of our work is that fossil cavities driven by transient protostellar outflows couple strongly to their surrounding cloud environments and drive turbulence within those clouds. We conjecture that the coupling can occur in two ways. First to produce turbulence in an initially quiescent medium interactions in the form of collisions between fossil cavities are required. Second, once fully developed turbulence exists, individual fossil outflows can still give their energy and momenta back to the turbulent flow through the disruption of the cavity shell. The simulations presented here address the second process. With regard to the initial generation of turbulence via interacting fossil cavities we note a few key points. First, as shown in Cunningham et al. (2006a), in a typical cluster environment one can expect that every parcel of ambient gas will be overrun by an outflow at least once. One can also estimate occurrence of collisions. For clouds containing a volumetric density of N_* protostellar driven bipolar outflows each of radius R, length L, and lifetime t_out, the probability P of cavity overlap increases as

$$\begin{equation} P = \left[ N_* \pi R^2 L \frac{t_{\rm out}}{t_{\rm cloud}} \right]^2, \end{equation} \tag{ 9 }$$

where t_cloud is the age of the cloud, which we assume fixed. Thus, the older the outflow the higher the probability that cavities will overlap. In the work of Cunningham et al. (2006a), direct numerical simulations showed that active outflows, where jet beams collided, were less efficient at stirring ambient material than the sum of two similar outflows that do not collide. This was because the beams could not pass through each other and forced a shock mediated redirection in which significant kinetic energy of the flow was lost to radiation. The conclusion of that study was not, however, that collisions of all kinds were ineffective. Instead the authors conjectured that interactions between slower, longer lived (and larger) fossil cavities could be the source of turbulent driving. Indeed the important difference between active outflows and fossil outflow cavities is that in the later case overlapping outflows can penetrate through each other causing disruption but not redirection. This is due to the fact that the fossil cavities are hollow and contain no active momentum-bearing jet beam. We note that the work of Matzner (2007) provides an analytical description of just such a process. Using dimensional arguments, Matzner (2007) derives a length scale at which impulsively driven (spherical) outflows will interact. Given an ambient cloud of density ρ and volumetric outflow rate S, and with each outflow driven with momenta I, Matzner finds that the turbulence will be driven at a scale length

$$\begin{equation} L_{\rm turb} \sim \frac{I^{1/7}}{\rho ^{1/7} S^{1/7}}. \end{equation} \tag{ 10 }$$

Typical values of ρ, I, and S appropriate for NGC 1333 yield an interaction scale L_turb ∼ 1. Simulation studies presented in Li & Nakamura (2006) and Nakamura & Li (2007) support the basic scenario that multiple outflows can drive turbulence in a protstellar cloud. In addition, Carroll et al. (2009) have directly confirmed Matzner's prediction with a longer effective scale length obtained due to jet collimation of L_eff ≈ 3L_turb. Thus, we conclude that multiple overlapping, interacting fossil cavities can drive turbulence in an initially quiescent medium.

In this paper, we have shown that interactions need not occur for a fossil cavity, already in a turbulent medium, to resupply its environment with turbulent energy. Our identification of the disruption and eventual subsumption of the fossil cavities demonstrates the pathways through which an outflow interacting with only ambient turbulence can act to support a turbulent cascade. This can be seen in the time evolution of the Mach number of the short jet pulse driven outflow cavity (run 1) shown in the left column of Figure 5, revealing that the outflow cavity becomes significantly disrupted and subsumed into the ambient turbulence before becoming subsonic. Thus, in the present simulations, the seeding of flow instability via ambient turbulence provides the mechanism for disruption.

The scenario through which protostellar outflows drive turbulence proposed and studied in this work contrasts with that of Banerjee et al. (2007). The primary reason for the difference is that the two works strive to answer different questions. Whereas Banerjee et al. (2007) seek to determine the efficacy of isolated outflows in a quiescent media, the present work focuses on the efficacy of outflows toward sustaining an initial state of supersonic turbulence. Our scenario may also be contrasted with the assumptions of Banerjee et al. (2007) who began their simulations with a single isothermal jet propagating into a quiescent ambient medium. They seek to test if a single jet can leave behind enough supersonic material to power turbulence on its own. Models arising from an initially turbulent state should be viewed as more realistic as clouds are never quiescent. Molecular clouds will always be born in an environment characterized by a range of motions on different size scale due to formation, i.e., either gravitational collapse or collisions between large-scale flows. In this regard, we recall that both observational and theoretical studies of molecular outflows show mass–velocity relations that scale as M(v) ∝ v^k with =1.3 < k < 2.5 (Lada & Fich 1996; Smith et al. 1997; Lee et al. 2001). Thus, while the bulk of material will be moving at low velocity, much of the gas set in motion is supersonic. This is because the ambient environment is laced with a volume-filling supersonic flow which will form a "floor" or background into which the outflow propagates. Thus, the assumptions underlying the claims that jets cannot drive turbulence should be considered in light of the environments these jets inhabit. Such consideration as shown here and in the work of Nakamura & Li (2007) and Carroll et al. (2009) shows that jets can be quite effective in driving and maintaining turbulent cascades.

We also note other factors that may cause lower estimates of flow speeds and momentum imparted to entrained gas in the models of Banerjee et al. (2007). The outflows presented in this paper are driven at high Mach number, M_j ≈ 90, relative to the sound speed of the ambient gas, which is to be expected in YSO jets. The continuously driven jet models in this paper achieve a characteristic bow shock compression ratio of 150, consistent with the analytic estimate given by equation (5) of Blondin et al. (1990). In contrast, the simulations of Banerjee et al. (2007) consider isothermal hydrodynamic jets driven by flows with Mach numbers, M_j, of 5 and 10 (with a single M_j = 20 case interacting with a cloud). Because the M_j = 5 and 10 jet models are driven with density contrast ρ_j/ρ_a = 1, the inward and outward facing shocks at the jet Mach disk propagate at half the driving speed of the jet with shock Mach number M_s of 2.5 and 5. Consideration of Rankine–Hugoniot jumps across a one-dimensional shock yields the prediction that these models would achieve compression ratios of M²_s, or 6.25 and 25, consistent with the density plots presented in that study. Thus lower compression ratios are achieved than is expected from higher Mach number jet models. Higher Mach number models exhibit more highly compressed flow structures that take longer to decelerate to a subsonic regimes.

Our models also consider the effects of radiative line cooling, which is more characteristic of actual protostellar outflow environments. The simulations in Banerjee et al. (2007) utilize nearly isothermal flow conditions via a choice of the polytropic index γ = 1.0001. While useful for approximating the compressive effects of cooling, it is important to note the limitations of such an approach which can, in general, produce differences from a condition in which cooling is applied. While the isothermal approximation will produce high compressions near shocks, which are characteristic of flows subject to rapid radiative energy loss, it does not capture the effect of mechanical cooling in rarefactions. Such effects are particularly important for transient outflow models that employ a time-dependent jet launching. This effect can be illustrated by considering the one-dimensional impulsive withdrawal of a piston with constant speed v_p from a uniform gas with sound speed c_o. The leading front of the rarefaction has sound speed c, given as

$$\begin{equation} c = \min \left[c_o - \frac{\gamma -1}{2}|v_p|,0\right]. \end{equation} \tag{ 11 }$$

This simplifies to c = c_o for models that use an isothermal equation of state. While isothermal models adequately approximate the effect of rapid radiative energy loss and high compression across strong shocks, such models do not capture the effect of mechanical cooling across rarefactions.

Strong rarefactions which ensue after the driving source of the outflow expires are characteristic of outflow cavity/transient outflow models (Cunningham et al. 2006b). Because post-rarefaction isothermal flows are not subject to the effect of mechanical cooling, the sound speed in these flow regions is overestimated and will likely cause a systematic underestimation of the Mach number in such regions of the flow field. This effect biases the results of such models in a manner that is adverse to the development of supersonic motions. These effects should be studied in future works to ascertain their impact on the volumetric distribution of Mach numbers in single jets driven into a quiescent medium.

The results presented in this paper show that outflow-driven cavities contribute to turbulence in their immediate environment, provided there exists some mechanism to disrupt the (mostly) laminar flow patterns that would result from an undisturbed outflow in a quiescent environment. The particular form of the disruption in this work is that of decaying turbulence in the ambient environment with an initial turbulent velocity that is nominally 10% that of the propagation speed of the outflow bow shock. We have shown that even structures which are denser and propagate more quickly than the ambient turbulence are significantly disrupted, despite their higher momentum density. Outflow structures can be disrupted by the development of instabilities such as radiative (Sutherland et al. 2003) and thin shell (Vishniac & Ryu 1989; Vishniac 1994) modes which have growth times far shorter than the evolutionary timescales of the flow. The wide spectrum of motions present in the ambient turbulence provides a multimode spectrum of perturbation seeds that initiate the growth of these instabilities. Extinct outflow cavities thus become fully subsumed into the background turbulence in the time required for these disruptions to cross the cavity. Observation of outflow-driven cavities similar to the work carried out in Quillen et al. (2005) therefore becomes increasingly confused with the ambient turbulence during this transition. For the case of highly collimated outflow structures, which are prototypical of low-mass star formation, the relevant crossing distance is the radius of the roughly cylindrical outflow bow shock. The disruption time for low-mass protostellar outflows is therefore short relative to the turbulent crossing time of their parent cores. After the disruption time, the outflow cavities act to support turbulence in the local environment on a driving scale comparable to the size of the outflow cavity. We conclude, therefore, that protostellar outflows provide an efficient form of dynamical feedback in low-mass star-forming cores as suggested by Norman & Silk (1980).

We postulate that the action which disrupts an outflow cavity need not be fully developed turbulence in the ambient environment. The interaction of an extinct outflow cavity with another outflow or outflow remnant would produce similar results. The result of Cunningham et al. (2006a) showed that large radiative energy loss from the interaction of multiple active outflow actually inhibits the mechanical support. Here, we emphasize cavity–cavity interactions can also drive turbulence and this mechanism, as shown by Carroll et al. (2009). Future work should also focus on ecological studies of multiple star formation (Li & Nakamura 2006; Nakamura & Li 2007). Simulation of such systems over long timescales with sufficient resolution could determine the rate of steady-state turbulent decay resulting from a steady rate of outflow formation.

We thank Chris Matzner and Chris McKee for extremely useful discussions. Hector Arce, John Bally, and Pat Hartigan 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. Andrew Cunningham received support under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under contact DE-AC52-07NA27344. We also thank the University of Rochester Laboratory for Laser Energetics and funds received through the DOE Cooperative Agreement No. DE-FC03-02NA00057.