Investigations of Protostellar Outflow Launching and Gas Entrainment: Hydrodynamic Simulations and Molecular Emission
We investigate protostellar outflow evolution, gas entrainment, and star formation efficiency using radiation-hydrodynamic simulations of isolated, turbulent low-mass cores. We adopt an X-wind launching model, in which the outflow rate is coupled to the instantaneous protostellar accretion rate and evolution. We vary the outflow collimation angle from and find that even well collimated outflows effectively sweep up and entrain significant core mass. The Stage 0 lifetime ranges from 0.14-0.19 Myr, which is similar to the observed Class 0 lifetime. The star formation efficiency of the cores spans 0.41-0.51. In all cases, the outflows drive strong turbulence in the surrounding material. Although the initial core turbulence is purely solenoidal by construction, the simulations converge to approximate equipartition between solenoidal and compressive motions due to a combination of outflow driving and collapse. When compared to simulation of a cluster of protostars, which is not gravitationally centrally condensed, we find that the outflows drive motions that are mainly solenoidal. The final turbulent velocity dispersion is about twice the initial value of the cores, indicating that an individual outflow is easily able to replenish turbulent motions on sub-parsec scales. We post-process the simulations to produce synthetic molecular line emission maps of CO, CO, and CO and evaluate how well these tracers reproduce the underlying mass and velocity structure.
Subject headings:stars: formation, stars:low-mass, stars:winds, outflows, ISM: jets and outflows, turbulence
During formation, stars launch high-velocity, collimated mass outflows that impact the local gas and global cloud environment. Since the complex physics of outflow launching occurs on sub-AU scales and happens during the protostellar phase while the forming star is enshrouded by dense dust and gas, the process is observationally difficult to investigate. High-resolution observations by ALMA, which can probe AU scales may shed light on this process (Arce et al., 2013), but current observational data fail to probe the appropriate scales. Consequently, the details of the “central engine” that powers outflows is not well understood.
Mass outflows from forming stars bear on a number of fundamental open questions in star formation. By entraining and unbinding infalling material, outflows reduce the efficiency at which gas accretes onto protostars (Matzner & McKee, 2000; Arce & Goodman, 2002; Arce & Sargent, 2005). This has implications for the stellar characteristic mass and stellar initial mass function (Krumholz et al., 2012). Outflow activity may also partially explain the difference between the distribution of prestellar core masses and stellar masses, which differ by a factor of (Enoch et al., 2007; Alves et al., 2007; Rathborne et al., 2009). All molecular clouds are observed to be turbulent, but the origin of this turbulent energy is unknown. Protostellar outflows are one mechanism that may replenish turbulent motions (Bally et al., 1999; Swift & Welch, 2008; Nakamura et al., 2011; Nakamura & Li, 2007; Carroll et al., 2009; Cunningham et al., 2009; Carroll et al., 2010; Wang et al., 2010; Li et al., 2010; Hansen et al., 2012).
Several analytic theories have been proposed to explain outflow characteristics (Shu et al., 1994; Pelletier & Pudritz, 1992; Li & Shu, 1996; Matzner & McKee, 1999). These theories are based on the strong coupling between magnetic fields and accreting gas. As a result of the subsequent magnetic enhancement and field winding during the collapse process, mass is ejected at high-velocities along the field lines.
Numerical simulations provide an alternative avenue to explore outflow launching and characteristics. However, the complexity of the physics and the wide range of scales involved (0.01 AU-1 pc) prohibit a first-principles approach in all aspects of the problem. Consequently numerical simulations focus on one of three regimes. Studies investigating the outflow engine model 100 AU scales with sub-AU resolution (Lovelace et al., 2010; Lii et al., 2012; Price et al., 2012; Tomida et al., 2013). Calculations focusing on the formation of an individual star and role of outflow activity model larger volumes at lower resolution, e.g., 1AU (Duffin & Pudritz, 2009; Seifried et al., 2011; Machida & Hosokawa, 2013). In some cases, models may adopt a simplified model for the outflow launching (Lee et al., 2000, 2001; Rosen & Smith, 2004; Banerjee et al., 2007; Offner et al., 2011). Some of these calculations are able to follow the launching of outflows self-consistently, but generally resolution is not sufficiently high to create a well collimated jet on sub-AU scales (e.g., Duffin & Pudritz 2009; Seifried et al. 2011; Machida & Hosokawa 2013) and calculations may be too computationally expensive to follow for long evolutionary times (e.g., Price et al. 2012; Tomida et al. 2013). Finally, studies investigating the interaction of outflows with the parent cloud and subsequent impact on cluster properties focus on pc scales and adopt a sub-grid model for outflow launching (Nakamura & Li, 2007; Carroll et al., 2009, 2010; Offner et al., 2012; Hansen et al., 2012; Wang et al., 2010; Li et al., 2010; Krumholz et al., 2012). Such models are generally motivated by previous analytic work and coupled to hydrodynamic quantities such as the stellar mass and accretion rate.
In this paper we focus on the intermediate scales (20 AU -0.2pc) on which outflows interact with their local environment. We perform numerical simulations of the collapse of turbulent, dense low-mass cores. We adopt a sub-grid model for the outflow launching, which allows us to vary the outflow opening angle and investigate the impact on protostellar accretion and evolution. In section 2 we describe the simulation methodology. In section 3 we describe the simulation outcomes, including the distribution of gas and protostellar properties. We then post-process the simulations to obtain molecular line emission in section 4. We summarize our results in section 5.
2. Numerical Methods
We perform the simulations using the ORION Adaptive Mesh Refinement (AMR) code (Klein, 1999; Truelove et al., 1998). We include self-gravity and compute the gas temperature by solving the equations of radiative transfer in the Flux-Limited Diffusion (FLD) approximation (Krumholz et al., 2007). This treatment assumes the gas and dust are perfectly thermally coupled, and we adopt a standard solar metallicity dust composition (Semenov et al., 2003). In very hot gas, the dust sublimates and temperatures can no longer be computed with the FLD solver. Instead, we include the treatment for atomic line cooling described by Cunningham et al. (2011), which implicitly solves for gas temperatures exceeding K.
We initialize the calculations with a 10 K sphere of radius pc and uniform density g cm, which corresponds to a total mass of 4 . This cold core is pressure confined by a hot, low-density medium. The velocity field of the cold, dense sphere is initialized with a random field of perturbations having Fourier modes . The Mach number of the gas, , is chosen so that the core satisfies the linewidth-size relation (McKee & Ostriker, 2007). The calculations use a basegrid of 64, where the sphere is resolved to level 2, and we insert five or seven AMR levels such that the minimum cell size is 26 or 6.5 AU. Additional refinement occurs if either of three criteria is met. To avoid artificial gravitational fragmentation, the grid is refined if the gas density exceeds the local Jeans density, , where is the local sound speed and is the Jeans number (Truelove et al., 1997). Any gas with density is always refined to level 2 (AU) whether or not it exceeds the local Jeans density. These two criteria ensure that the outflow cavity walls and entrained material are always followed at higher resolution. Finally, higher resolution is added if there is a strong gradient in the radiation energy density: . This criterion facilitates the convergence of the radiation solver and forces cells in the radiatively heated region near the protostar to be refined even if the gas densities are low. Protostars are always contained within cells on the finest level, which is where the accretion and wind launching takes place. We adopt outflow boundary conditions so that outflowing material can freely escape the domain. The simulation properties are summarized in Table 1.
As the calculation proceeds the initial turbulence decays, allowing the gas to gravitationally collapse. When the Jeans density is exceeded on the maximum level a star particle is inserted (Krumholz et al., 2004). These particles represent individual protostars and follow a sub-grid stellar evolution model that includes accretion luminosity down to the stellar surface, Kelvin-Helmholz contraction, and nuclear burning (Offner et al., 2009). The protostars are also endowed with a model for the launching of protostellar outflows based upon (Matzner & McKee, 1999). This model has previously been used in other outflow studies (Cunningham et al., 2011; Offner et al., 2011; Krumholz et al., 2012; Hansen et al., 2012; Offner et al., 2012). The model is characterized by three dimensionless parameters that specify the outflow ejection efficiency, outflow velocity, and momentum distribution. The mass ejection rate, , gives the fraction of accreting gas that is launched in the wind. This fraction is observationally uncertain (e.g., Plunkett et al. 2013), but the disk wind (Pelletier & Pudritz (1992) and X-wind (Shu et al., 1988) models predict 0.1-0.33. Here, we adopt = 0.2-0.3. Consequently, of the infalling gas accretes onto the star, while is launched in an outflow.
For the outflow launching velocity, we adopt the Keplerian velocity near the stellar surface, , where is the protostellar mass, is the protostellar radius, and is a model parameter relating to the wind acceleration. The X-wind and disk wind models both predict , while surveys of outflows suggest a range of values, , with no apparent dependence on spectral class (Cunningham et al., 2011). Here we adopt such that . In this model, the launch velocity increases fairly strongly with mass, and so the outflow speed increases over the course of the simulation. We set the wind temperature to K, the temperature of an ionized wind.
The outflow direction is set by the direction of the angular momentum vector of the protostar. This is in turn determined by the time-integrated angular momentum of the accreting gas, a quantity that depends upon the turbulent properties of the core and evolves over the calculation. Thus, the wind direction is not fixed but is self-consistently dictated by the hydrodynamic evolution of the accretion flow. However, we do not find large variations in the outflow direction or precession, which has been observed in some outflows (Ybarra et al., 2006; Wu et al., 2009). The outflow axis is generally parallel to the angular momentum vector of the accretion disk.
We define the effective opening angle of the wind, , following Matzner & McKee (1999). They write the angular momentum distribution of the outflow in terms of the polar angle measured from the protostar’s rotational axis :
Based on observed low-mass protostars, Matzner & McKee (1999) infer that and suggest a fiducial value of = 0.01, which produces strong collimation. Here, we adopt 0.01, 0.03, and 0.1. Small values of result in most of the outflow momentum being deposited in the grid within a few cells along the rotational axes.
We inject the wind into the refined cells with radial distances from the protostar. By construction, the injection region is exterior to the accretion region, so that the wind does not impede accretion along the disk midplane.
We evolve each of the calculations for 0.5 Myr, which corresponds to approximately four core free-fall times. Observationally, the embedded phase lasts Myr (Evans et al., 2009). Thus, we follow the core evolution until most of the initial core mass is either accreted or expelled by the outflow.
At the final time, the calculations each have a single protostar. Additional protostars may form during the evolution, but protostars that approach within four cells of the primary and also have a mass are merged. We find that one or more additional protostars do form in our studies, generally in the early stage of collapse before the formation of an organized Keplerian disk around the primary, but their masses remain below and they don’t survive as independent protostars.
To check for convergence, we perform an additional calculation, th0.1fw0.3h, with 7 AMR levels and AU, which otherwise has identical model parameters to th0.1fw0.3 (see Table 1). Model th0.1fw0.3h produces one fewer fragment at early times and the fragment masses are smaller than those in model th0.1. This leads us to conclude that higher resolution models would converge to a single or, at most a close, binary system.
Figure 1 shows the density distribution for models th0.1fw0.3 and th0.1fw0.3h, which are very similar. Figure 2 shows the protostellar masses, wind launched, domain gas mass and mass of high-velocity gas as a function of time. The trends are very similar, although the protostar in th0.1fw0.3 is formed earlier with a slightly higher mass. The outflow launching and protostellar luminosity depend upon the protostellar radius, which in turn depends upon the interior stellar state. A slight difference in mass in the two cases leads to one protostar progressing to deuterium burning before the other, which impacts the protostellar radius and introduces small differences between gas velocities. The domain mass evolution in the two runs is sufficiently similar that the black thin line for run th0.1fw0.3 is hidden beneath the th0.1fw0.3h line.
In principle the amount of entrainment may be enhanced or suppressed depending on the simulation resolution. Run th0.1fw0.3h contains more cells at higher resolution and has a larger resolution gradient between the protostar and the ambient medium. However, Figure 2 shows that the fraction of high velocity material, , is similar in the two runs, which suggests that gas entrainment is not strongly sensitive to resolution to the extent we vary it here.
For completeness we also consider run cl.th0.01, which follows the formation of a cluster of protostars and employs the same outflow launching model. The details of the simulation are described in Offner et al. (2012). We use this simulation to contrast the turbulence generated from an ensemble of outflows with that of an isolated protostar.
|ModelaaThe model name, gas mass, simulation resolution on the maximum level, effective opening angle, fraction of accreted mass ejected in a wind, and the final time of the simulation. The molecular cores have initial radii of 0.064pc, temperature 10 K, and turbulence given by the linewidth-size relation km s.||()||(AU)||(rad)||(Myr)|
|cl.th0.01bbTurbulent clump simulation of a forming star cluster, which has pc, initial temperature 10 K, and km s. See Offner et al. (2012) for additional details.||180||128||0.01||0.3||0.3|
|ModelaaModel name, final stellar mass, final remaining gas on the domain, remaining gas that has K, the star formation efficiency , the final 3D velocity dispersion, the Stage 0 lifetime, and the star formation rate per free fall time.||()||()||()||(km s )||SFR|
3.1. Outflow Morphology
The simulated outflows exhibit a variety of morphologies over time. Figures 3 and 4 show volume renderings of the gas velocity at six times for models with and . The initial turbulence in the core produces significant asymmetry between the upper and lower outflow lobes as well as asymmetry about the launching axis. Since the outflow itself is driven symmetrically, most of these differences arise from the interaction between the outflow and envelope gas. For example, the top, right panel of Figure 4 shows that the lower lobe breaks out of the core earlier, while the upper lobe remains confined. Similar outflow asymmetries were also found in Offner et al. (2011).
At intermediate times the km s gas appears relatively similar in the two cases, despite the different launching angles. This is because the outflow successfully sweeps up and entrains core material in both cases. In fact, much of the high-velocity (red) gas shown is entrained material, which has mixed with the launched material. The figures show that gas at wider angles moves systematically slower. Sometimes slower moving, less collimated material is attributed to a “disk wind” mechanism, wherein the outflowing material is launched at larger radii with slower velocities (Pelletier & Pudritz, 1992). However, here there is only a single launching mechanism that gives rise to a continuum of outflow morphologies and properties.
At late times, when most of the core gas has been accreted or dispersed, little entrainment occurs and the opening angle differences are more apparent. For example, the last panel of Figure 4 shows a very collimated outflow compared to that in Figure 3.
The evolution of the gas velocity in the surrounding material is quite striking. Most of the initial core turbulence decays within 0.2 Myr and the core gas has velocities km s, which is uncolored in Figures 3 and 4. However, the fraction of diffuse gas with velocities km s increases with time. This is generally warm () gas that has been shocked but not unbound by the outflows. This demonstrates the potential for outflows to drive turbulence in their surroundings, and we explore this more quantitatively in Section 3.5.
3.2. Mass Evolution
The momentum distribution of the outflow and its interaction with the surrounding core gas have direct implications for the protostellar accretion rate and mass. As discussed in §2, we define the outflow launching velocity, , as a function of the protostellar properties. Figure 5 shows the launching velocity versus time for the four models. The curve discontinuities correspond to different interior stellar states with different amounts of deuterium burning; three of these states are indicated for model th0.01 (see the Appendix in Offner et al. 2009 for more details). At early times the stellar state is somewhat sensitive to accretion rate variations. Jumps in the wind launching rate are due to changes in the stellar radius, which depends upon the deuterium burning state. However, after the first 0.05 Myr, the accretion details have little impact. The gradual increase in from Myr is due to the slow rise in . At late times, the model converge to similar values because the protostellar masses and radii are comparable. More comprehensive stellar evolution models including variable accretion likewise predict that early accretion creates minimal scatter in protostellar properties after Myr (e.g., Hosokawa et al. 2011).
We find that the outflow evolution and gas entrainment is not linear in even though the simulations begin with identical initial core masses and turbulent velocity distributions. Some of the differences that result are due to differences in the early protostellar accretion, which depend on the extent the wind launching impacts the accretion flow. In the wider angle case, the infalling material interacts more with the outflow, leading to additional fragmentation. The widest angle wind run undergoes the most fragmentation episodes (), while the narrowest angle run experiences only one. If different numbers of small fragments accrete onto the protostar, this will translate into small changes in the primary mass and later evolution.
The early phase of accretion and disk building lasts for 0.06 Myr. During this time streams of gas feeding the inner region deposit material with different angular momenta, producing significant changes in disk structure and total angular momentum. After this time, the gas settles into a stable Keplerian accretion disk of radius 200 AU. We expect this size to be an upper limit since magnetic fields are not included (Commerçon et al., 2011; Seifried et al., 2011; Myers et al., 2013; Krumholz et al., 2013).
Figure 7 shows the column density distribution at different times for run th0.01. The other runs show a similar evolutionary progression. As shown in Figure 7, the outflow breaks out of the core around 0.2 Myr. The outflow broadens and entrains additional gas from 0.2-0.3 Myr. By 0.35 Myr, the initial core has been almost completely disrupted or accreted. Thereafter, the protostar accretes the remaining turbulent, cold gas (green), more slowly.
Figure 6 shows the accretion rates for the four runs as a function of time. The outputs have slightly different but comparable output intervals, with a median spacing of yr. All cases exhibit a decline in the mean accretion rate by more than an order of magnitude over the 0.3 Myr accretion time. Variability is largest at early times, which corresponds to the period of disk building and clumpy accretion. A decline corresponding to the disruption of the envelope is clearly visible in the accretion history of th0.01 (top panel in Figure 6) and in the others to a lesser degree. Some features in the average accretion, which span 0.03-0.05 Myr, correspond to periods when the protostar accretes a clump of material that is falling inward. An example of this is marked on the bottom panel of Figure 6. The interaction between the outflow and the envelope creates a significant amount of clumpiness in the residual cold gas as shown in Figure 7.
At late times, the accretion disk is stable and accretion variability generally remains within a factor of of the average accretion rate. Radiative heating acts to suppress large scale instability, and the accretion rate from the envelope onto the disk is sufficiently low that we do not expect global instability to occur (Kratter et al., 2010; Offner et al., 2010). However, the simulations do not have sufficient resolution to resolve small planet-size clumps if they occur (e.g., Vorobyov 2013; Tsukamoto et al. 2013).
The differences between the runs is illustrated by Figures 8 and 9, which show the evolution of the protostellar masses, gas mass, and outflow mass as a function of time. Neither the evolution of protostellar mass nor the amount of high-velocity gas is monotonic with opening angle. This is partially due to differences in the early fragmentation and evolution and partially because of the initial turbulence within the core impacts the entrainment and evolution of the outflows. Small changes in the orientation of the outflow and the protostellar mass result in different results.
As expected, the run with the highest outflow efficiency has the lowest protostellar mass and the most mass on the domain at the end of the run. The run that has , th0.1, also exhibits less mass at high-velocities because the gas is less collimated and more of the outflow momentum is distributed at wider angles with lower velocities. Run th0.1 also concludes with more gas mass remaining on the domain. This suggests that the wider angles winds in the prescription are not as efficient at unbinding and expelling core gas. Counterintuitively, the narrower opening angle cases, th0.01 and th0.03, conclude with less gas on the domain. They have higher amounts of high-velocity gas even though the total mass launched is slightly lower. This demonstrates that narrowly collimated outflows can entrain and expel dense material. This is a different result than found by Banerjee et al. (2007), who investigate a highly-collimated jet interacting with a magnetized medium. However, these authors adopted smooth, non-turbulent initial conditions; density inhomogeneities is an important component in coupling the outflow momentum and surrounding gas (e.g., Cunningham et al. 2009; Wang et al. 2010).
3.3. Protostellar Evolutionary Stage
One of the main goals of observing individual protostars is to determine the evolutionary stage of the source. The earliest evolutionary period, Stage 0, during which most of the accretion occurs, is defined to be the length of time during which the envelope mass, , is greater than the protostellar mass, (Andre et al., 2000; Enoch et al., 2008). While may be relatively easy to obtain from continuum emission measurements, the protostellar mass is nearly impossible to determine because protostellar masses cannot be directly measured. In rare cases, high-resolution observations of accretion disk properties can provide indirect estimates of the protostellar mass (Tobin et al., 2012), but in most cases the stage must be inferred through indirect means such as spectral energy distribution (SED) modeling (Whitney et al., 2003; Robitaille et al., 2007). A more observationally convenient method of source classification involves using the effective bolometric temperature of the SED to separate colder, more extincted sources, which are allegedly younger, from warmer sources, which appear to have less surrounding gas and are thought to be older. However, projection effects (Whitney et al., 2003; Offner et al., 2012) and non-monotonic accretion rates (Dunham & Vorobyov, 2012) can make younger sources appear older than they are and vice versa. For example, a source viewed along the outflow cavity will have less extinction along the line-of-sight and appear younger than it would if it were viewed through an edge-on accretion disk. Similarly, a source undergoing an accretion burst will appear brighter and may look more evolved. Comparing SEDs to analytic models for their evolutionary stage combined with some direct imaging can help improve age estimates, but SED classes are often assigned for individual sources without imaging and are thus likely poorly correlated with evolutionary stage.
In contrast, when protostars are treated as an ensemble, the individual errors associated with projection and variability may average out, allowing determinations of the stage lifetime (Evans et al., 2009). Using the statistics of protostars in local regions, the Class 0 phase, which represents the earliest and highest accretion phase, is inferred to have a lifetime of 0.1 Myr (Enoch et al., 2009; Maury et al., 2011). Evans et al. (2009) estimated a combined Class 0 and Class I lifetime of 0.5 Myr, which suggests a Class I lifetime of 0.3-0.4 Myr. However, these times are very uncertain because they depend upon an assumed disk lifetime, which is Myr.
Simulations provide one avenue for exploring the underlying Class and Stage lifetime. Offner et al. (2012) demonstrated that at early times, the stage inferred for simulated forming protostars on the basis of SED characterization and modeling is generally correct, although, inferred properties such as the envelope and protostellar mass could be incorrect by more than a factor of 2. Machida & Hosokawa (2013) inferred Stage 0 lifetimes ranging from yr.222Machida & Hosokawa (2013) refer to this as the “Class 0 lifetime”, but because their definition is based on the protostellar mass and not on SED characteristics such as the bolometric temperature or spectral slope, the times they report are more accurately the Stage lifetimes.
As plotted in Figure 9 and indicated in Table 2, we find Stage 0 lifetimes of 0.14-0.23 Myr. These are slightly larger than found by Machida & Hosokawa (2013), but they are comparable to estimates of the observed Class 0 lifetime. The range in values underscores that the length of the physical stage depends upon a range of initial properties, including outflow collimation, the initial turbulence and core properties. In calculations of isolated cores determining the core mass as a function of time is trivial. However, in simulations of clusters (e.g., Hansen et al. 2012), the envelope mass is connected with the cloud environment and may change with time due to additional accretion. Thus, in clustered conditions, simulations, like observations, must wrestle with the challenge of defining a “core”.
3.4. Star Formation Efficiency
Understanding the efficiency at which molecular gas turns into stars is an important theme in star formation, which has repercussions for the evolution of molecular clouds and the origin of the stellar IMF. On cloud scales, only a few percent of the gas turns into stars per dynamical time (Tan et al., 2006; Krumholz & Tan, 2007). This appears to be a consequence of a combination of large scale supersonic turbulence, magnetic support, and stellar feedback. The efficiency of dense gas on the core scale is much higher. Since the decay time for turbulence is short and magnetic diffusion and reconnection reduce the field strength, the efficiency of dense gas is mainly limited by feedback and stellar multiplicity. Comparison between the stellar IMF and the observed distribution of core masses suggests an efficiency of (Enoch et al., 2008; Alves et al., 2007), which naively implies that one-third of the gas in a core is converted into stellar mass. This is consistent with theoretical estimates (Matzner & McKee, 2000); although, if a core produces a binary or multiple star system the actual efficiency per core will be higher (Holman et al., 2013).
Observationally, the derivation of the star formation efficiency is complicated by a number of factors including whether cores are bound, how many (if any) stars they form, and time-dependent effects. Consequently, represents a fairly uncertain ensemble average of the efficiency. Numerical simulations of isolated cores allow an unambiguous estimate of the star formation efficiency. In these calculations we find efficiencies of 0.36-0.43 at 0.5 Myr as defined by the total mass in the protostar relative to the initial envelope mass (see Table 2). If we ignore the remaining cold gas on the domain at 0.5 Myr, some of which may accrete onto the protostars if the calculations were run longer, we find 0.41-0.51. These estimates suggest lower and upper limits on the efficiency at , respectively. These values are comparable to the estimated observed core efficiency, although this is very uncertain. Since each of these calculations formed only a single star, it makes sense that the efficiencies are larger than the average value obtained when comparing the core mass function to the IMF.
The efficiency of star formation is more meaningful when considered together with some characteristic timescale, because, in principal, if there is no cloud dispersal mechanism then nearly all the gas will turn into stars on a sufficiently long timescale. The star formation rate (SFR) per free fall time is defined as the fraction of mass that turns into stars per free fall time (Krumholz & McKee, 2005). Here, the SFR, where is the protostellar mass at 0.5 Myr and is the time over which the accretion occurred. Unimpeded global collapse gives SFR=1. As shown in Table 2, we find that the SFR given that the free fall time of the initial core is Myr.333If the initial starless phase in the calculation is included then the SFR. These calculations, which only consider an isolated core, yield a fairly low SFR. This implies that a combination of outflow feedback and core turbulence can contribute an order of magnitude to the global star formation inefficiency.
3.5.1 Velocity Dispersion
A number of authors have investigated the ability of outflows to contribute to the global energy budget in numerical simulation of clouds (Nakamura & Li, 2007; Nakamura et al., 2011; Carroll et al., 2009; Wang et al., 2010; Hansen et al., 2012). However, little attention has been devoted to the study of outflow-driven turbulence on core scales, which we address here.
Figures 3 and 4 qualitatively demonstrate that the outflows in these simulation successfully drive turbulence on core scales. Figure 10 shows the root-mean-squared gas velocity dispersion, , as a function of time for the three runs. At early times, declines due to turbulent dissipation until the protostar forms. Once the outflow is launched, increases non-monotonically as the high-velocity outflow gas interacts with the surrounding dense gas. The largest changes occur during the phase in which the outflow has not fully broken out of the core. At late times the velocities reach a quasi-steady state of about twice the initial velocity dispersion. We conclude that a single outflow is sufficient to offset the turbulent decay on sub-pc scales and can inject significant energy throughout the local region.
3.6. Solenoidal and Compressive Motions
Any velocity field can be deconstructed into two orthogonal components: a compressive field () and a solenoidal field (). Physically, these components indicate what fraction of the motions are “squeezing” versus “stirring.” The compressive mode acts in a similar sense to gravitation in that it forces gas to higher densities. For star formation, the fraction of compressive motions may have bearing on the gas density and column density distributions (Federrath et al., 2010).
In order to compute the solenoidal and compressive modes, we flatten the AMR data to a fixed grid (level 2). This allows us to easily perform Fourier Transforms of the data. Although this means that we neglect smaller scale motions, we find that the results are not significantly different if we adopt resolution instead, which implies that the larger scale modes dominate the totals.
In these simulations the turbulence is initialized using a purely solenoidal random field. During the starless phase of the simulation, during which the turbulence decays and gravitational contraction begins, the solenoidal motions decline until compressive motions dominate (bottom panels of Figure 10). The launching of the outflow increases the solenoidal component, essentially by reducing gravitational collapse and creating circulation of material to larger scales. After 0.1 Myr, the ratio of the solenoidal to compressive velocity dispersion approaches one and the modes appear to be in equipartition. Some of the compressive motion is due to ongoing gravitational contraction, as can be seen by the changing ratio for . This suggests that protostellar outflows inject turbulent motions that are more solenoidal than compressive.
The bottom panel of Figure 10 shows the ratio of the solenoidal to compressive rms velocity along the three cardinal directions. The details of the distribution depend upon the orientation of the outflow with respect to the view.444Note that the individual ratios do not add to one because The top panel of Figure 10 shows the individual components of the solenoidal and compressive motions. These are higher than the total because , which means that at any given position or can be much greater than .
3.7. Core and Clump Turbulence Comparison
In order to compare with turbulent outflow driving in a clustered environment, we analyze turbulent properties in a larger simulation, cl.th0.01, which is forming a cluster of stars (see Table 1). The clump simulation adopts periodic boundary conditions to model a pc piece of a larger molecular cloud. Consequently, a number of protostars form in fairly close proximity to one another, which is illustrated in the top panel of Figure 11. Like the previously discussed core simulations, the initial turbulence is allowed to decay and the only kinematic input is from protostellar outflows. Unlike the previous calculations the gas is not centrally condensed and gravitational motions are less ordered and more localized. Since gravity contributes to the compressive mode, we can use clth0.01 to examine the outflow driving in a context where gravity is less dominant. We note that the cl.th0.01 simulation has periodic boundary conditions, and high-velocity material cannot leave the domain, which leads to an excess of turbulence in comparison to the case where outflowing material is allowed to escape the domain. As shown in the top panel of Figure 11, the outflow driving is not directly discernible from the gas column density distribution.
The bottom panels of Figure 11 show and the solenoidal and compressive components as a function of time. As in the isolated case, the total dispersion and the dispersion of the solenoidal and compressive components increase once stars begin forming and launching outflows. However, the solenoidal and compressive modes are not as similar, and the outflows appear to drive much more solenoidal motion than compression. This could be because the cloud is not centrally condensed or globally collapsing, which is the case in the isolated core runs, and consequently, gravity contributes less compressive motion. This trend is illustrated by the solenoidal to compressive ratio, which is greater than unity at all times and in all three directions.
4. CO Molecular Emission Maps
In this section, we investigate the outflow morphology and velocity distribution inferred from emission in several observational tracers. In order to generate synthetic emission maps, we post-process the simulations at selected intervals with radmc-3d (Dullemond, 2012). We use the non-LTE Large Velocity Gradient approach (Shetty et al., 2011), which computes the level populations given input 3D density and temperature fields. We flatten the AMR data to a fixed 256 resolution (0.001 pc). To interpolate over velocity jumps between grid cells we use “doppler catching” with 0.025, which interpolates the velocity field such that velocity changes between points do not exceed 0.025 times the local line width. For the warm, high-velocity gas component, the thermal line width can exceed 1.4 km s, so this small value is helpful for hot cells that border denser, cold gas. We adopt a constant “microturbulence” value of 0.05 km s in order to include sub-resolution line broadening caused by unresolved turbulence. We assume a constant molecular abundance except in the hot (K) where the molecules are assumed to have much lower abundances or be completely dissociated. The molecular excitation and collisional data are adopted from the Leiden atomic and molecular database (Schöier et al., 2005). We use the Rayleigh-Jeans approximation () of the Planck function to compute the brightness temperature corresponding to the output flux.555It is standard for observers to employ this approximation to obtain the observed brightness temperature even if the emission is not in the Rayleigh-Jeans limit.
As the second most abundant molecule in star-forming regions after H, CO is the most commonly used tracer for cm gas. CO is also useful as a tracer of the lower density, molecular component of protostellar outflows (Arce et al., 2007). The CO isotopologues, CO and CO, have lower abundances and higher critical densities and so are useful tracers for examining slightly higher density gas, where CO is optically thick. Here, we assume a constant CO abundance of for cells with temperatures below 900 K. Gas above this temperature is either strongly shocked, which may dissociate the CO or low-density, high-pressure confining gas, which was not part of the original molecular cloud core. We adopt CO/CO=62 and CO/H =.
Figure 12 shows the integrated emission for CO(1-0), CO(1-0), and CO(1-0). The CO emission is very bright and clearly traces the initial core gas as well as the outflow. Clumps and streams of gas can be seen being ejected from the inner region. CO only traces the outflow cavities and dense core center; very little of the initial core or ejected material is apparent. The CO picks out only the inner cavity structure, and the cavity appears more rounded than in the other two tracers.
Figure 13 shows the integrated CO(1-0) emission in a set of velocity channels for four different times. At the earliest time, the high-velocity gas is very localized. Once the outflow has broken free of the core, the emission takes on a distinctive “v” shape, which has been observed in a variety of observed outflows (e.g., IRAS3282, Arce & Sargent 2006). The outflow emission characteristics at Myr are also similar to those reported by Quillen et al. (2005) who mapped NGC 1333 in CO and found cavity sizes of pc and velocity widths of a couple .
The amount of gas that appears at higher velocities depends partly on the orientation of the outflow with respect to the line-of-sight. The figure shows a view in which the outflow is titled degrees with respect to the line-of-sight. At smaller angles of orientation, the outflow cavity appears more collimated. The emission with velocities close to 10 km s is mainly concentrated near the protostar. Much higher velocity gas is present than appears in the maps because much of the highest velocity gas has temperatures exceeding the cutoff of K. At the launch point, the wind is assumed to be ionized ( K). This gas quickly begins to cool and much of the high velocity gas is in fact around a few thousand Kelvin. In principle this hot, high-velocity gas may emit in CO if some of it is entrained and heated core gas, since CO does not dissociate until K. However, given that the formation timescale of CO is much longer than the dynamical time of the outflow, we expect the CO abundance of this material to be much lower and the contribution to the total CO (1-0) emission from this hot molecular gas to be small.
While there are a few instances of discrete clumps of material, the outflows do not exhibit the apparent episodicity exhibited by outflows like HH 46/47 (Arce et al., 2013). Although the simulated protostars do experience some accretion variability, it is generally less than a factor of a few over consecutive time steps. Observed episodic clumps may be created by a very high-velocity jet component, which is also not well resolved here or is not reproducible with our simple outflow model.
Toward late times, the higher velocity gas is more diffuse and patchy. The outflow is not directly apparent in the channel maps, instead the emission traces the residual turbulent core gas. Whether this complex velocity structure is resolvable in observations depends on the total cloud velocity dispersion, which may hide small features, and the beam resolution.
The CO isotopologues are somewhat imperfect tracers of the velocity dispersion. Figure 14 shows the average second velocity moment as a function of time for two simulations computed directly from the simulated data and from the CO, CO, and CO emission. For the raw simulation data, the second velocity moment is simply the mass-weighted root mean squared velocity. For the emission data, the average second moment is defined as the average of the map second moments:
where is the brightness temperature, is the first moment at location and is the number of pixels in the map with non-zero dispersion. The same data is shown where the emission has been convolved with a 5” Gaussian beam assuming that the source is 250 pc distant. This is analogous to an outflow in Perseus observed with the Combined Array for Research in Millimeter-wave Astronomy (CARMA).
Generally, the velocity dispersion increases slightly as the outflow expands into the core and then decreases as the envelope is accreted or expelled (see also Figure 10). The velocity dispersions of the isotopologues are generally expected to follow , where the smaller line width tracers probe relatively smaller sizes and higher density gas. Here, we find that the line widths vary as expected. Convolving the emission with a beam preserves the order of the line width but tends to make the dispersions slightly larger. None of the tracers track the mass-weighted velocity dispersion well. The CO dispersion is the most similar, although the actual line-of-sight dispersion may be either larger or smaller than the observed value. The difference between the synthetic line width and the simulation velocity dispersion is mainly due to varying excitation, which is a function of the local temperature, density, and optical depth. Since we assume constant abundances for the CO isotopologues, we expect that in this comparison the intensity of the synthetic emission maps is better correlated with the underlying gas density than would be true for actual observations.
Figure 15 shows a map of the second moment for each of the isotopologues at one simulation snapshot. The dark areas in the plot indicate sightlines with no cold, dense gas. The CO traces the outflow best, and clumps of low-density, high-velocity gas are apparent in the upper outflow cavity, which is distinctly V-shaped, and on the left side. The other two tracers highlight more of the wider angle material, which is part of the outer outflow cavity. This material has second-moments of 1-2 km s. The CO probes a smaller fraction of the gas, but its emission and estimated velocity dispersion still strongly resemble that of CO.
We perform a set of radiation-hydrodynamic simulations of isolated forming protostars with a model for outflow launching in order to investigate outflow evolution, turbulent driving, and star formation efficiency. Our outflow model is based on the X-wind magnetized wind model, in which the outflow rate is coupled to the protostellar accretion rate and evolution. The gas temperature is calculated using flux-limited diffusion including radiative feedback from the forming protostar and atomic cooling at high temperatures. In the simulations, we vary the outflow collimation angle from and find that even well-collimated outflows effectively entrain gas and drive turbulence. Due to the inclusion of radiative heating, after the early infall phase, accretion onto the protostar is generally smooth and only varies by a factor of a few on timescales of yr. Consequently, while the gas accretion and outflow launching is somewhat variable, the outflow does not undergo large episodic bursts.
The final turbulent velocity dispersion is about twice the initial value, which demonstrates that an individual outflow easily replenishes turbulent motions on 0.1 pc scales. Although the initial turbulence in the simulations is purely solenoidal, we find that the resulting gas motions are approximately equally solenoidal and compressive. Some of the gas compression is due to gravitational collapse, so it appears that outflows drive motions that are predominantly solenoidal rather than compressive. We confirm this conclusion by analyzing an additional simulation of a cluster of forming protostars, which is not centrally condensed.
We post-process the simulations with radmc-3d to produce synthetic molecular line emission maps in CO, CO, and CO. The emission morphologies of the simulated outflow cavities appear similar to observed outflows. The line width is anti-correlated with the tracer critical density, and as expected, CO is a better probe of the outflow gas. The effective dispersions of the different tracers are similar but not well-correlated with the actual dispersion computed using the 3D density information. We find that the results are qualitatively similar if the emission is convolved with 5” beam, although the inferred dispersions are larger.
In future work we plan to address outflow entrainment and turbulent driving including the effects of magnetic fields. Some numerical simulations can produce magnetically launched outflows self-consistently. However, either resolution is insufficient to probe a highly collimated jet component or the calculation is evolved for a very short time due to computational time constraints. Most magnetic calculations also assume perfect coupling between the magnetic field and gas (ideal MHD), and no 3D calculation includes all the necessary physics modeling magnetic diffusion and reconnection. More accurate and complete future studies are needed to fully understand outflow launching, entrainment, and turbulence.
- Alves et al. (2007) Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17
- Andre et al. (2000) Andre, P., Ward-Thompson, D., & Barsony, M. 2000, Protostars and Planets IV, 59
- Arce & Goodman (2002) Arce, H. G. & Goodman, A. A. 2002, ApJ, 575, 928
- Arce et al. (2013) Arce, H. G., Mardones, D., Corder, S. A., Garay, G., Noriega-Crespo, A., & Raga, A. C. 2013, ApJ, 774, 39
- Arce & Sargent (2005) Arce, H. G. & Sargent, A. I. 2005, ApJ, 624, 232
- Arce & Sargent (2006) —. 2006, ApJ, 646, 1070
- Arce et al. (2007) Arce, H. G., Shepherd, D., Gueth, F., Lee, C., Bachiller, R., Rosen, A., & Beuther, H. 2007, Protostars and Planets V, 245
- Bally et al. (1999) Bally, J., Reipurth, B., Lada, C. J., & Billawala, Y. 1999, AJ, 117, 410
- Banerjee et al. (2007) Banerjee, R., Klessen, R. S., & Fendt, C. 2007, ApJ, 668, 1028
- Carroll et al. (2010) Carroll, J. J., Frank, A., & Blackman, E. G. 2010, ApJ, 722, 145
- Carroll et al. (2009) Carroll, J. J., Frank, A., Blackman, E. G., Cunningham, A. J., & Quillen, A. C. 2009, ApJ, 695, 1376
- Commerçon et al. (2011) Commerçon, B., Hennebelle, P., & Henning, T. 2011, ApJ, 742, L9
- Cunningham et al. (2009) Cunningham, A. J., Frank, A., Carroll, J., Blackman, E. G., & Quillen, A. C. 2009, ApJ, 692, 816
- Cunningham et al. (2011) Cunningham, A. J., Klein, R. I., Krumholz, M. R., & McKee, C. F. 2011, ApJ, 740, 107
- Duffin & Pudritz (2009) Duffin, D. F. & Pudritz, R. E. 2009, ApJ, 706, L46
- Dullemond (2012) Dullemond, C. P. 2012, RADMC-3D: A multi-purpose radiative transfer tool, astrophysics Source Code Library
- Dunham & Vorobyov (2012) Dunham, M. M. & Vorobyov, E. I. 2012, ApJ, 747, 52
- Enoch et al. (2009) Enoch, M. L., Evans, N. J., Sargent, A. I., & Glenn, J. 2009, ApJ, 692, 973
- Enoch et al. (2008) Enoch, M. L., Evans, II, N. J., Sargent, A. I., Glenn, J., Rosolowsky, E., & Myers, P. 2008, ApJ, 684, 1240
- Enoch et al. (2007) Enoch, M. L., Glenn, J., Evans, II, N. J., Sargent, A. I., Young, K. E., & Huard, T. L. 2007, ApJ, 666, 982
- Evans et al. (2009) Evans, N. J., Dunham, M. M., Jørgensen, J. K., Enoch, M. L., Merín, B., van Dishoeck, E. F., Alcalá, J. M., Myers, P. C., Stapelfeldt, K. R., Huard, T. L., Allen, L. E., Harvey, P. M., van Kempen, T., Blake, G. A., Koerner, D. W., Mundy, L. G., Padgett, D. L., & Sargent, A. I. 2009, ApJS, 181, 321
- Federrath et al. (2010) Federrath, C., Roman-Duval, J., Klessen, R. S., Schmidt, W., & Mac Low, M.-M. 2010, A&A, 512, A81
- Hansen et al. (2012) Hansen, C. E., Klein, R. I., McKee, C. F., & Fisher, R. T. 2012, ApJ, 747, 22
- Holman et al. (2013) Holman, K., Walch, S. K., Goodwin, S. P., & Whitworth, A. P. 2013, MNRAS, 432, 3534
- Hosokawa et al. (2011) Hosokawa, T., Offner, S. S. R., & Krumholz, M. R. 2011, ApJ, 738, 140
- Klein (1999) Klein, R. I. 1999, Journal of Computational and Applied Mathematics, 109, 123
- Kratter et al. (2010) Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 708, 1585
- Krumholz et al. (2013) Krumholz, M. R., Crutcher, R. M., & Hull, C. L. H. 2013, ApJ, 767, L11
- Krumholz et al. (2012) Krumholz, M. R., Klein, R. I., & McKee, C. F. 2012, ApJ, 754, 71
- Krumholz et al. (2007) Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J. 2007, ApJ, 667, 626
- Krumholz & McKee (2005) Krumholz, M. R. & McKee, C. F. 2005, ApJ, 630, 250
- Krumholz et al. (2004) Krumholz, M. R., McKee, C. F., & Klein, R. I. 2004, ApJ, 611, 399
- Krumholz & Tan (2007) Krumholz, M. R. & Tan, J. C. 2007, ApJ, 654, 304
- Lee et al. (2000) Lee, C., Mundy, L. G., Reipurth, B., Ostriker, E. C., & Stone, J. M. 2000, ApJ, 542, 925
- Lee et al. (2001) Lee, C., Stone, J. M., Ostriker, E. C., & Mundy, L. G. 2001, ApJ, 557, 429
- Li & Shu (1996) Li, Z. & Shu, F. H. 1996, ApJ, 468, 261
- Li et al. (2010) Li, Z.-Y., Wang, P., Abel, T., & Nakamura, F. 2010, ApJ, 720, L26
- Lii et al. (2012) Lii, P., Romanova, M., & Lovelace, R. 2012, MNRAS, 420, 2020
- Lovelace et al. (2010) Lovelace, R. V. E., Romanova, M. M., Ustyugova, G. V., & Koldoba, A. V. 2010, MNRAS, 408, 2083
- Machida & Hosokawa (2013) Machida, M. N. & Hosokawa, T. 2013, MNRAS, 431, 1719
- Matzner & McKee (1999) Matzner, C. D. & McKee, C. F. 1999, ApJ, 526, L109
- Matzner & McKee (2000) —. 2000, ApJ, 545, 364
- Maury et al. (2011) Maury, A. J., André, P., Men’shchikov, A., Könyves, V., & Bontemps, S. 2011, A&A, 535, A77
- McKee & Ostriker (2007) McKee, C. F. & Ostriker, E. 2007, ARA&A, 45, 565
- Myers et al. (2013) Myers, A. T., McKee, C. F., Cunningham, A. J., Klein, R. I., & Krumholz, M. R. 2013, ApJ, 766, 97
- Nakamura et al. (2011) Nakamura, F., Kamada, Y., Kamazaki, T., Kawabe, R., Kitamura, Y., Shimajiri, Y., Tsukagoshi, T., Tachihara, K., Akashi, T., Azegami, K., Ikeda, N., Kurono, Y., Li, Z., Miura, T., Nishi, R., & Umemoto, T. 2011, ApJ, 726, 46
- Nakamura & Li (2007) Nakamura, F. & Li, Z. 2007, ApJ, 662, 395
- Offner et al. (2009) Offner, S. S. R., Klein, R. I., McKee, C. F., & Krumholz, M. R. 2009, ApJ, 703, 131
- Offner et al. (2010) Offner, S. S. R., Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 725, 1485
- Offner et al. (2011) Offner, S. S. R., Lee, E. J., Goodman, A. A., & Arce, H. 2011, ApJ, 743, 91
- Offner et al. (2012) Offner, S. S. R., Robitaille, T. P., Hansen, C. E., McKee, C. F., & Klein, R. I. 2012, ApJ, 753, 98
- Pelletier & Pudritz (1992) Pelletier, G. & Pudritz, R. E. 1992, ApJ, 394, 117
- Plunkett et al. (2013) Plunkett, A. L., Arce, H. G., Corder, S. A., Mardones, D., Sargent, A. I., & Schnee, S. L. 2013, ApJ, 774, 22
- Price et al. (2012) Price, D. J., Tricco, T. S., & Bate, M. R. 2012, MNRAS, 423, L45
- Quillen et al. (2005) Quillen, A. C., Thorndike, S. L., Cunningham, A., Frank, A., Gutermuth, R. A., Blackman, E. G., Pipher, J. L., & Ridge, N. 2005, ApJ, 632, 941
- Rathborne et al. (2009) Rathborne, J. M., Lada, C. J., Muench, A. A., Alves, J. F., Kainulainen, J., & Lombardi, M. 2009, ApJ, 699, 742
- Robitaille et al. (2007) Robitaille, T. P., Whitney, B. A., Indebetouw, R., & Wood, K. 2007, ApJS, 169, 328
- Rosen & Smith (2004) Rosen, A. & Smith, M. D. 2004, MNRAS, 347, 1097
- Schöier et al. (2005) Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369
- Seifried et al. (2011) Seifried, D., Banerjee, R., Klessen, R. S., Duffin, D., & Pudritz, R. E. 2011, MNRAS, 417, 1054
- Semenov et al. (2003) Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611
- Shetty et al. (2011) Shetty, R., Glover, S. C., Dullemond, C. P., & Klessen, R. S. 2011, MNRAS, 412, 1686
- Shu et al. (1988) Shu, F. H., Lizano, S., Ruden, S. P., & Najita, J. 1988, ApJ, 328, L19
- Shu et al. (1994) Shu, F. H., Najita, J., Ruden, S. P., & Lizano, S. 1994, ApJ, 429, 797
- Swift & Welch (2008) Swift, J. J. & Welch, W. J. 2008, ApJS, 174, 202
- Tan et al. (2006) Tan, J. C., Krumholz, M. R., & McKee, C. F. 2006, ApJ, 641, L121
- Tobin et al. (2012) Tobin, J. J., Hartmann, L., Chiang, H.-F., Wilner, D. J., Looney, L. W., Loinard, L., Calvet, N., & D’Alessio, P. 2012, Nature, 492, 83
- Tomida et al. (2013) Tomida, K., Tomisaka, K., Matsumoto, T., Hori, Y., Okuzumi, S., Machida, M. N., & Saigo, K. 2013, ApJ, 763, 6
- Truelove et al. (1997) Truelove, J. K., Klein, R. I., McKee, C. F., Holliman, II, J. H., Howell, L. H., & Greenough, J. A. 1997, ApJ, 489, L179+
- Truelove et al. (1998) Truelove, J. K., Klein, R. I., McKee, C. F., Holliman, II, J. H., Howell, L. H., Greenough, J. A., & Woods, D. T. 1998, ApJ, 495, 821
- Tsukamoto et al. (2013) Tsukamoto, Y., Machida, M. N., & Inutsuka, S. 2013, ArXiv e-prints
- Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9
- Vorobyov (2013) Vorobyov, E. I. 2013, A&A, 552, A129
- Wang et al. (2010) Wang, P., Li, Z., Abel, T., & Nakamura, F. 2010, ApJ, 709, 27
- Whitney et al. (2003) Whitney, B. A., Wood, K., Bjorkman, J. E., & Wolff, M. J. 2003, ApJ, 591, 1049
- Wu et al. (2009) Wu, P.-F., Takakuwa, S., & Lim, J. 2009, ApJ, 698, 184
- Ybarra et al. (2006) Ybarra, J. E., Barsony, M., Haisch, Jr., K. E., Jarrett, T. H., Sahai, R., & Weinberger, A. J. 2006, ApJ, 647, L159