COLLECTIVE OUTFLOW FROM A SMALL MULTIPLE STELLAR SYSTEM

Thomas Peters; Pamela D. Klaassen; Mordecai-Mark Mac Low; Martin Schrön; Christoph Federrath; Michael D. Smith; Ralf S. Klessen

doi:10.1088/0004-637X/788/1/14

1. INTRODUCTION

Molecular outflows in high-mass star-forming regions appear to differ from those around low-mass stars not just in their strength, but also in their lack of collimation. Beuther & Shepherd (2005) suggest that more massive stars appear to have less collimated outflows, although conclusive evidence for this awaits high-resolution observations of massive star-forming regions by ALMA. The nearest massive star-forming regions with strong outflows, such as DR 21 (Davis & Smith 1996) and Cepheus A (Narayanan & Walker 1996) are observed to have complex jet structures rather than the single, well-collimated jets typically observed from low-mass stars.

Our published study of the interplay of magnetic fields and self-gravity in the accretion flow around stars with masses M_* > 15 M_☉, which we call high-mass stars in subsequent discussion, suggests that accretion disks around such stars may be vulnerable to disruption by gravitational torques from their own accretion flows, preventing them from driving magnetocentrifugal jets at all (Peters et al. 2011). However, we determined that observed outflows have substantially higher momentum than either the large-scale magnetic tower jet driven by the accretion flow (Peters et al. 2011), or the ionization-driven outflow (Peters et al. 2012).

Zinnecker & Yorke (2007) review the theoretical controversy over whether massive stars form by monolithic collapse in isolated cores or accretion from a cluster environment containing many other stars. Multiple low- and intermediate-mass stars form in the accretion flow onto high-mass stars in analytic (Kratter & Matzner 2006) and numerical models (Bonnell et al. 2004; Smith et al. 2009b; Peters et al. 2010a, 2010b; Wang et al. 2010). The exceptions to this are numerical models starting with initial core masses of 100–300 M_☉ rather than 1000 M_☉ (Krumholz et al. 2007b; Commerçon et al. 2011; Myers et al. 2013), particularly with strongly centrally concentrated initial density profiles, which Girichidis et al. (2011) show to suppress fragmentation regardless of other physical processes. Indeed, the lack of isolated young O stars has long been noticed observationally (Gies 1987; Gvaramadze et al. 2012), although whether every OB star formed in a group does remain contested (Oey & Lamb 2012).

Peters et al. (2012) suggested that the combined effect of jets driven from the inner disks of the secondary stars could explain the observed properties of outflows from massive-star forming regions, even absent an outflow from the most luminous and massive star. Jets generally can be approximated to have power proportional to the masses of their stellar sources (Matzner 2002). In a region with a Salpeter (1955) initial mass function, the cumulative jet power will thus be proportional to the cumulative mass, which is indeed dominated by the low-mass stars regardless of whether the highest-mass stars actually have a jet or not. The collective action of the intermediate-mass stars should begin to dominate the outflow from the region even before enough mass has accumulated on the most massive object for it to begin emitting significant amounts of ionizing radiation, or for all of the low-mass stars to have formed: Peters et al. (2010a) show that the lower-mass stars form later than the more massive ones in their simulations.

The combined action of jets was first simulated in a massive collapsing region by Li & Nakamura (2006) and Nakamura & Li (2007). This work was extended by Wang et al. (2010), who coupled the outflow momentum feedback to sink particles. Cunningham et al. (2011) and Krumholz et al. (2012) presented models with combined radiative and outflow feedback. These studies included decaying turbulence in the collapsing gas, but neglected rotation, thus maximizing the support provided by the outflows against collapse by minimizing their alignment. However, simulations of massive star formation including turbulence show the angular momentum vectors of protostars in dense groups are initially closely correlated on length scales of a few tenths of a parsec, and only become more randomly oriented as the objects grow in mass and accrete more distant material (Jappsen & Klessen 2004; Fisher 2004). Even a model including not only turbulence, but also a better treatment of radiative heating, and higher spatial resolution than our model also finds such almost planar structures, extending at least to the size of the star forming region in our model of 1500 AU (see Figure 2 of Krumholz et al. 2007a).

We describe a numerical simulation of this early stage of outflow driving in an initially rotating cloud with turbulence only induced by its own gravitational instability. Since we do not include background sources of turbulence, our core collapses to the center and builds up a rotationally flattened structure there, in which the entire stellar group forms, thus minimizing the support, but maximizing outflow alignment. We expect that more realistic turbulent initial conditions would lead to fragmentation of our cloud on large scales and the formation of several collapsing regions with globally misaligned angular momentum vectors, but at least initially aligned protostellar spin axes in the densest star-forming regions. There is observational evidence for outflow alignment from well-separated protostars with distinct outflows on scales of a few tenths of a parsec in the vicinity of DR 21 (Davis et al. 2007) and in Source G of W49A (Smith et al. 2009a), which has suggested that they formed from the same flattened rotating cloud.

We compare the result to the observed outflows from the nearby young, massive star-forming regions Cepheus A and DR 21. In later work, we will describe the subsequent development of the H ii region within the outflow.

In Section 2, we describe our numerical model and in Section 3, we describe results, including simulated observations. In Section 4, we compare our results to the observations in order to determine the consistency of our scenario with available evidence.

2. NUMERICAL SIMULATION

We here present the first three-dimensional radiation-hydrodynamical simulation of massive star formation that simultaneously include feedback by protostellar outflows as well as by heating of both ionizing and non-ionizing radiation. We use the adaptive-mesh code FLASH (Fryxell et al. 2000) with sink particles (Federrath et al. 2010) and our improved version of the hybrid-characteristics ray-tracing method (Rijkhorst et al. 2006; Peters et al. 2010a). For more details on our numerical technique as well as its limitations, we refer to Peters et al. (2010a).

The initial conditions are identical to our previous simulations (Peters et al. 2010a, 2010b, 2010c, 2011). We start from a 1000 M_☉ molecular cloud that has an initial temperature T = 30 K as well as a core of constant density ρ = 1.27 × 10⁻²⁰ g cm⁻³ within a radius of r = 0.5 pc and then drops as r^−3/2 out to a radius of r = 1.6 pc. The cloud initially rotates as a solid body with angular velocity ω = 1.5 × 10⁻¹⁴ s⁻¹. Magnetic fields are not included in this simulation. We resolve the collapse of the cloud on the adaptive mesh with a minimum cell size of 98 AU. Sink particles form at a cut-off density of ρ_crit = 7 × 10⁻¹⁶ g cm⁻³ and an accretion radius of r_sink = 590 AU.

The sink particles carry an intrinsic angular momentum (or spin). Their angular momentum vector changes during accretion of gas such that the total angular momentum of the system (sink particle and gas) is conserved. Thus, if L_acc is the angular momentum of the accreted material, S the intrinsic sink angular momentum before accretion and L and L' the orbital angular momentum of the sink particle before and after accretion, then the new intrinsic sink angular momentum satisfies S' = S + L − L' + L_acc. Note that L' is determined by mass, center of mass, and linear momentum conservation during accretion (Federrath et al. 2010).

We launch the protostellar outflows around the sink particles by injecting 10% of the gas that is currently being accreted into the cells surrounding the sink particle, as suggested by Königl & Pudritz (2000) as a consequence of the angular momentum relation in cold, thin disks. The material is added to the cells within a cone of height 1567 AU (or 16 grid cells) and an opening angle 15° aligned to the angular momentum vector of the sink particle. The mass loading of the cells is not homogeneous but smoothed such that it goes to zero at the boundaries of the cone both as function of radius and angle to avoid sharp density contrasts. We assume a footpoint of 0.5 AU (Ray et al. 2007; Pudritz et al. 2007) and eject the outflow material with twice the escape velocity $v_\mathrm{esc} = \sqrt{2 G M / r}$ at this radius, which is v_esc = 60 km s⁻¹ for M = 1 M_☉. The most massive star formed in the simulation reaches a mass of about M = 10 M_☉. This corresponds to a maximum outflow velocity of 380 km s⁻¹, which is on the conservative side of the observed range of velocities (e.g., Micono et al. 1998; Coppin et al. 1998). Based on observational constraints (e.g., Bacciotti et al. 2002; Bacciotti 2004), we also transfer 90% of the angular momentum of the accreted gas to the cells in the outflow cone. We stress that we do not simply set the velocities within the outflow cone to the outflow velocity, but rather add the outflow momentum to these cells. All variability in the outflow is due to variations in the accretion flow determined from the mass and angular momentum of accreted gas and to changes in the outflow axis determined from angular momentum vector of the sink particle. We impose no additional time dependence on the launching mechanism.

Our simulation is not isothermal but takes the actual heating–cooling balance between compression and radiative heating on the one hand and molecular, dust and metal line cooling on the other hand into account. Because of strong shocks, the gas in the outflows can heat up to T ≳ 10⁶ K. Such high temperatures lead to a prohibitively small hydrodynamical timestep. To circumvent this problem, we have artificially enhanced the cooling rates for temperatures in excess of T ⩾ 10^4.3 K. Material below this temperature range, such as photoionized gas, will stay completely unaffected. The only dynamical effect of the enhanced cooling rates is to prevent the shock-heated gas from reaching temperatures well beyond T ≈ 25, 000 K.

This procedure does remove a lot of thermal energy from the shock. To check whether this will have an impact on the outflow dynamics, we examine whether the outflowing gas is momentum- or energy-driven. Let v_jet be the velocity of the jet material and n_jet be the H₂ number density in the jet. The cooling time at the jet termination shock, assuming diatomic gas, is

$\begin{equation} t_c = \frac{7}{5} \frac{k_\mathrm{B} T_s}{n \Lambda (T_s)}, \end{equation} \tag{ 1 }$

where

$\begin{equation} T_s = \frac{5}{36} \frac{\mu }{k_\mathrm{B}} v_\mathrm{jet}^2 \end{equation} \tag{ 2 }$

is the shock temperature for γ = 7/5, n = 4.6 n_jet is the number density of free particles assuming fully ionized gas, and

$\begin{equation} \Lambda (T_s) \sim 10^{-22} \frac{{\rm erg\,cm}^3}{{\rm s}} \left(\frac{T_s}{10^6\,\mathrm{K}}\right)^{-0.7}. \end{equation} \tag{ 3 }$

The last relation is valid over the range of temperatures of roughly 10^4.5–10^6.5 K for solar metallicity (Mac Low & McCray 1988). Here, k_B is Boltzmann's constant, μ = 2.14m_p is the mean molecular weight with the proton mass m_p. Plugging in some typical numbers, v_jet = 100 km s⁻¹ and n_jet = 10⁻²⁰ cm⁻³, we arrive at t_c ≈ 3 × 10⁷ s. On the other hand, the crossing time of the jet is t_x = L_jet/v_jet with the jet length L_jet. For a snapshot near the beginning of the simulation, L_jet = 0.4 pc, we have t_x ≈ 10¹¹ s and thus t_x ≫ t_c. The difference between these two numbers is so large that t_x ≫ t_c still holds when slightly different representative numbers and later snapshots are considered. Hence, we can safely conclude that the outflow is physically completely in the momentum-driven regime so that the enhanced cooling does not further impact its dynamics.

We here report results from the first part of the simulation, before the stars reach masses larger than 10 M_☉. At these early times, the H ii regions only extend to a few grid cells around the sink particles. The adaptive mesh at the end of this part of the simulation has 487 million grid cells, and the simulation has consumed 1.73 million CPU hours.

3. ANALYSIS

We start our analysis with a discussion of the stellar group that forms during the simulation and the associated outflow feedback (Section 3.1). We then analyze the collective outflow that forms from the interactions of the individual stellar outflows in Section 3.2. In Section 3.3, we present synthetic observations of this outflow.

3.1. Accretion History and Outflow Feedback

The accretion history of the stellar group is shown in Figure 1. Four stars form during the simulation runtime, three of which are close to 10 M_☉. These three stars accrete at an average rate of ∼4 × 10⁻⁴ M_☉ yr⁻¹, while the first star that forms accretes at the slightly lower rate of ∼2 × 10⁻⁴ M_☉ yr⁻¹.

The momentum injected into the cloud by the protostellar outflows is shown in Figure 2. We plot both the instantaneous momentum ΔP at a single time step and the cumulative momentum P integrated over the age of the star. All four stars make a similar contribution to the total momentum of the collective outflow. The lower accretion rate and consequently lower mass of the first star results in a smaller contribution compared to the other three stars, so that the cumulative momentum output of the first star is only half of the momentum output of each of the other three stars. The fluctuations of the accretion rate result in a related variability of the instantaneous outflow momentum. This effect dominates the magnitude of the instantaneous outflow feedback, so that the outflow from a more massive star is not necessarily stronger than for a lower-mass star, although the escape velocity scales with stellar mass.

In our setup, all stars form in a rotationally flattened structure or toroid in the center of the simulation box (Peters et al. 2010a, 2010b). The stars revolve around the center of rotation while driving the outflows. Since the outflow cones of the different stars furthermore partially overlap, the momentum transport within the outflow is quite complex. The fact that one of the stars is injecting a little less momentum into the outflow than the other stars does not mean that a region in the outflow that moves slower than the rest can be identified. Furthermore, the momentum in the outflow is redistributed through shocks that form when the individual outflows collide (see Section 3.2).

Directly after their formation, the initial angular momenta S₀ of the sink particles point in the same direction as the disk angular momentum vector, which is oriented perpendicular to the disk. Later, as more gas falls onto the disk, gravitational instability sets in and leads to the formation of dense filaments. The stars are embedded in these filaments, which extend significantly above and below the disk, so the accreted gas has a preferred location with respect to the sink particles that is different from the disk midplane. Since the filaments are generally denser than the rest of the disk, the material of the filaments dominates the angular momentum balance of gas accreted onto the sink particles. Continuous accretion of gas above or below the disk plane results in a systematic trend for the angular momentum S of the sink particle. Over time, the favored accretion from a certain direction changes the direction of the sink angular momentum vectors significantly.

We quantify the changes in the sink angular momenta by looking at the angle ψ between the initial and the current angular momentum vector, $\psi = \arccos (\mathbf {j}_0 \cdot \mathbf {j})$ with j₀ = S₀/|S₀| and j = S/|S|. Since j₀ is orthogonal to the disk plane, ψ can be interpreted as an inclination angle. Figure 3 shows ψ as a function of time for all sink particles. Three of the four sink angular momentum vectors begin to wander after (t = 0.635 Myr), and one of them reaches a ψ of almost 50° by the end of the simulation.

**Figure 3.** Angle ψ between the initial and the current angular momentum vector of the sink particles as function of time. As the disk becomes gravitationally unstable and vertically extended filaments form, the angular momentum vectors become misaligned with the disk angular momentum.
Download figure:
Standard image High-resolution image

The formation of the filaments in the disk is displayed in Figures 4 and 5. The figures shows snapshots from the time when all sink particle angular momenta are still aligned perpendicular to the disk (t = 0.625 Myr), the time when the first angular momentum vectors begin moving (t = 0.631 Myr) and the very end of the simulation (t = 0.642 Myr) in face-on and edge-on slices, respectively. The correlation between the formation of filaments and the departure of the inclination angle ψ from ψ = 0° is clearly visible.

**Figure 4.** Density slices through the disk plane for three different snapshots. As more and more gas falls onto the disk, gravitational instability sets in and filaments form.
Download figure:
Standard image High-resolution image

**Figure 5.** Density slices perpendicular to the disk plane for the same snapshots as in Figure 4. One can discern the accumulation of gas above and below the disk as time proceeds.
Download figure:
Standard image High-resolution image

3.2. Evolution of the Collective Outflow

We show the time evolution of the collective outflow in Figure 6. The outflow develops immediately after the first star forms at t = 0.613 Myr and its leading edge leaves the simulation box at t = 0.635 Myr. In the middle of this time interval, at t = 0.624 Myr, the outflow only extends to a quarter of the simulation box, which shows that the outflow growth is not linear in time. The collective outflow is driven by the momentum of all the individual outflows, and this outflow feedback is highly episodic due to the coupling of the outflow momentum to the current accretion rate, whose strong variation is shown in Figure 2. We therefore do not expect a simple relation between the size of the outflow and the time since its formation.

Figure 6. Density slices through the collective outflow at four different snapshots (rows) and with three different magnification factors (columns). The masses of the stars in the group are given in the image.

(A color version of this figure and an animation of this figure are available in the online journal.)

Download figure:

Video Standard image High-resolution image

As can be seen in Figure 6, the quickly moving shock front (strong density enhancement) at the tip of the outflow is subject to instabilities that lead to the formation of finger-like structures. The gas density inside the outflow shows substantial variations between 10⁻²⁴ and 10⁻¹⁹ g cm⁻³. The sudden increase of the accretion rate onto one of the stars leads to intensified outflow feedback, which in turn triggers the formation of a shock front, which is visible as a density enhancement in the interior of the outflow. This shock then propagates through the entire outflow until it dissipates at the outflow boundary or leaves the simulation box. The movie provided in the online material shows that the variation of the outflow angle and the movement of the stars in the disk leads to a significant variation in the way the outflow material is injected near the disk plane. In particular, the collision of shocks launched from different stars and moving at different speeds or encountering each other at acute angles can be observed frequently during the entire evolution.

The collimation of the outflow decreases slightly with time. This observation can be related to the increasing growth of the inclination angle of the individual outflows (compare Figure 3). As the individual outflow inclinations rise, the effective opening angle of the collective outflow becomes larger, extending the base of the collective outflow.

The head of the collective outflow has a distinctive morphology. Figure 7 shows the time evolution of one of the outflow heads until the outflow head leaves the simulation box. The movie in the online material shows that the multiple heads continuously come and go during the evolution of the jet. Their origin becomes clear once one realizes that the velocity of the jet increases as it propagates down the density gradient of the envelope, by roughly a factor of three from early times to when the head leaves the box. The accelerating interface between the lower density jet gas and higher density envelope gas has an effective gravity with a sign opposite to that of the acceleration, so that the interface is Rayleigh–Taylor unstable. The head structure thus reflects the usual Rayleigh–Taylor bubble and spike morphology (e.g., Sharp 1984). As the minimum wavelength of this instability is determined only by viscosity, higher-resolution simulations will show finer small-scale structure until the effects of magnetic fields and physical viscosity are included and resolved. Fujita et al. (2009) show an example of this phenomenon in a blowout-induced Rayleigh–Taylor instability.

3.3. Synthetic Observations

We have generated synthetic observations of our simulated outflows to compare our results with previous simulations and observations. We present H₂ maps of the outflow in Section 3.3.1 and CO observations in Section 3.3.2. We have used our customized version of the three-dimensional adaptive-mesh radiative transfer code RADMC-3D⁸ to simulate the observations.

3.3.1. H₂ Emission Maps and Spectra

We study emission of the collective outflow in the H₂ S(1) 1–0 line at 2.1218 μm. To this end, we employ the method developed by Suttner et al. (1997). This scheme assumes statistical equilibrium, but not local thermodynamic equilibrium (LTE), for the vibrational states, and solves the full set of rate equations for the three level system consisting of the lowest energy states. To determine the H₂ density, we assume statistical equilibrium between H₂ formation on dust grains and the two temperature-dependent collisional dissociation processes H₂ + H → H + H + H and H₂ + H₂ → H₂ + H + H. We use the dust formation rate from Hollenbach & McKee (1979), and for the destruction processes we employ the interpolation formula from Glover & Mac Low (2007), which interpolates between the low-density (Lepp & Shull 1987; Shapiro & Kang 1987) and high-density (Mac Low & Shull 1986; Martin et al. 1998) limits of these rates. However, at shock fronts, although we account for collisional dissociation, we do not account for shock photodissociation, so our emission levels are likely upper limits at fast shocks such is in our jet heads.

Figure 8 shows the H₂ emission for the four snapshots displayed in Figure 6. The spatial scale of the four images is identical and corresponds to the left panels of Figure 6. It can be seen that the extent of the outflow is well represented by the density slices of Figure 6, but that multiple tips are present at the outflow head that are not visible in the selected slice plane. This is not surprising because the slice plane lies perpendicular to the disk that harbors the four stars (compare Figure 4). Since the H₂ emission is optically thin, the H₂ maps reveal the three-dimensional structure of the outflow better than the density slices. However, the emission from different regions can overlap, leading to confusion along the line of sight.

Figure 8. Intensity of the H₂ S(1) 1–0 line for the four snapshots of Figure 6. The spatial scale is the same in all four images.

(A color version of this figure and an animation of this figure are available in the online journal.)

Download figure:

Video Standard image High-resolution image

We can compare these maps to the three-dimensional outflow simulations by Völker et al. (1999). Their Hammer jet model includes precession of the outflow launching angle, a sinusoidal perturbation of the outflow velocity as function of time (pulsation), a velocity gradient across the launching area (shear) and a small opening angle at the outflow footpoint (spray). Our outflow sub-grid model also includes shear and spray as free parameters, but the precession and velocity variations of our outflows are determined self-consistently by the properties of our sink particles.

As in Völker et al. (1999), the H₂ emission comes predominantly from the outflow head. However, we also see significant emission from the bulk of the outflow because the amplitude of our velocity pulsations is not limited, while Völker et al. (1999) do not allow the outflow velocity to grow by more than a factor of two during a pulsation period. We have therefore much more shock-heated gas in the interior of our outflow. The Völker et al. (1999) Hammer model leads to regularly spaced mini-bow shocks that are aligned with the bulk flow.

Our outflows look considerably more irregular for a number of reasons. First, our precession angles grow to markedly larger values than the few degrees considered by Völker et al. (1999; see Figure 3). Second, the momentum output from each star can vary by more than an order of magnitude (see Figure 2). Third, we have four outflows interacting with each other as they propagate instead of just a single outflow expanding in a quiescent medium.

In Figure 9, we show the outflow heads from Figure 7 in H₂ emission. The head structure with its multiple tips looks very similar to the density slice in general. Shocks running through the outflow can be identified through their stronger emission, because the shocks are generally denser than the outflow interior. The brightness of the emission at the outflow tip depends on when a shock has reached this tip for the last time. This is because the shock rapidly heats up the gas, which then subsequently cools. The H₂ emission then diminishes with the gas temperature. This effect can be nicely seen in the animated version of Figure 8.

We study the kinematics for a simulation snapshot at t = 0.630 Myr, when the outflow is still completely contained inside the simulation box. To compare the kinematics of the collective outflow to observations, we have generated maps at an inclination of 45° and 90° at a spatial resolution of 191 AU and a spectral resolution of 9.82 km s⁻¹ (63 channels). This corresponds to the observing conditions of Hiriart et al. (2004). We show channel maps of the observed outflow in Figures 10 and 11 for the two inclination angles. As expected, in the edge-on maps the emission is strongest near the rest velocity, while the 45° maps have larger emission at finite velocities. The finger-like structure of the outflow heads is visible in almost all channels.

**Figure 10.** Edge-on channel maps of the H₂ S(1) 1–0 line for the collective outflow at t = 0.630 Myr. Spectra through the marked points are shown in Figure 14.
Download figure:
Standard image High-resolution image

**Figure 11.** Channel maps at an inclination of 45° of the H₂ S(1) 1–0 line for the collective outflow at t = 0.630 Myr. Spectra through the marked points are shown in Figure 15.
Download figure:
Standard image High-resolution image

We present a position–velocity (PV) diagram for the 90° map in Figure 12 and for the 45° map in Figure 13. The slices were taken perpendicular to the disk through the center of the group, corresponding to the horizontal midline in Figures 10 and 11 channel maps. The emission covers the PV plane only fragmentarily, which may hint at the different components of the collective outflow. However, it is impossible to identify the four individual outflows because of the strong interaction between them. The radial spikes visible in the 45° PV diagrams are typical of pulsation and density enhancements in the outflow (Smith et al. 1997; Völker et al. 1999), in our case driven by the variations in accretion rate.

Spectra for positions marked in the PV diagrams for the 90° view (in Figure 10) and the 45° view (in Figure 11) are displayed in Figures 14 and 15, respectively. The observation at 45° shows more complex spectra since the line-of-sight velocities are much greater in this case. We find a variety of spectral forms, ranging from a sharp single peak through broad single peaks to multiply peaked spectral profiles. The spectrum is a strong function of position, and there is no apparent systematic structure, except that the spectrum becomes more complex in regions where shock fronts overlap, which generally corresponds to the brightest regions in Figure 11.

**Figure 14.** Spectra taken through the three positions in the outflow viewed from 90° indicated in Figure 10. The position of the spectra and the relative maximum intensities are indicated.
Download figure:
Standard image High-resolution image

**Figure 15.** Spectra taken through six positions in the outflow viewed from 45° indicated in Figure 11. The position of the spectra and the relative maximum intensities are indicated.
Download figure:
Standard image High-resolution image

3.3.2. CO Observations of the Collective Outflow

We model emission from the J = 2–1 transition of CO assuming LTE, with the Einstein coefficient taken from the Leiden Atomic and Molecular Database (Schöier et al. 2005), and a constant ratio of H₂ to CO abundance [H₂]/[CO] = 1 × 10⁻⁴. ALMA observation simulations were carried out using CASA (McMullin et al. 2007), version 4.1. The intensities of the output image cubes from RADMC-3D were converted to Jansky units by assuming a source distance of 3 kpc. The source coordinates were chosen for comparison to a southern hemisphere analog of DR 21 (with a declination of −42°); however, all coordinates in the resulting figures are given in spatial offsets from the phase center. A southern declination was necessary for properly sampling the UV plane.

We simulated observing CO (J = 2–1) at a spectral and spatial resolution of 0.5 km s⁻¹ and 0 farcs 22 using one of the full ALMA configurations bundled into CASA (alma.out08.cfg). The simulations were run assuming four hours of integration, and a precipitable water vapor column of 1 mm. This is consistent with good weather observing for Band 6 (230 GHz) ALMA observations, allowing an accurate representation of atmospheric noise in the visibilities. The online sensitivity calculator estimates a rms noise level of 4 mJy beam⁻¹, and we used a threshold slightly higher than this (6 mJy beam⁻¹) as our 1σ level when cleaning. Cleaning was done non-interactively, without the use of clean boxes. We used natural weighting, and our final image cubes spanned 384 × 384 × 256, where the last dimension indicates the number of spectral channels.

We show the first moment (intensity-weighted velocity) maps of the snapshots from Figure 6 at an inclination of 30° in Figure 16. Note that the size of the mapped region grows with the outflow scale. Also notice that the outflow has already started leaving the simulation box in the third snapshot, and that a significant portion of it is missing in the fourth. The inferred mass and kinematics of the collective outflow in these cases is, of course, only a lower bound on the true values.

The CO maps prominently exhibit the outflow heads with multiple tips in each of them. The maximum outflow velocity grows with time, but since the maximum value in the first moment maps is attained in the heads, this trend seems to break with the fourth snapshot. In actual fact, the outflow velocity is nearly constant along the outflow. It is the increased density near the outflow boundary that, via the intensity-weighting in the first moment maps, leads to the impression that the velocity in the head was higher. Note that no property of the outflow morphology suggests that the outflow might be driven by several protostars. Although the head does have multiple tips, they come and go in a way consistent with Rayleigh–Taylor instability (see Section 3.2).

A PV slice through the third outflow from Figure 16 is shown in Figure 17. One can distinguish multiple components, but again it is not possible to clearly associate them with the four protostars. The PV slice shows that the emission from the collective outflow can be decomposed into three parts: (1) small-scale low-velocity gas, (2) large-scale high-velocity gas, and (3) large-scale medium-velocity gas. Component (1) is produced by the disk-like structure that contains the four protostars. Components (2) and (3) trace different pieces of the outflow, which are possibly produced by different protostars. The CO bullets correspond to shocks running through the outflow and are also visible in Figure 16.

The image cubes for each snapshot were analyzed to determine the outflow mass, momentum, and energy. Using a kinematic age for the outflow, we also determined an outflow mechanical luminosity and mass loss rate for each snapshot. The total outflow mass was determined by summing the mass in each velocity bin outside of the central 2 km s⁻¹. The gas mass in an individual velocity bin was determined by multiplying the average intensity by a multiplication factor based on the ambient temperature (in this case 30 K). This results in the average number of emitting particles per unit area, assuming that the emission is optically thin. This average number of emitting molecules was then multiplied by the emitting area to determine the total number of emitting molecules. This was then further multiplied by the abundance ratio for CO, and the mean molecular weight 2.3. More details on this type of calculation can be found in Klaassen & Wilson (2007). The kinematic age of each snapshot was calculated based on the extent of the outflow and the intensity-weighted velocity at the edge of the outflow (12, 30, 36, and 31 km s⁻¹ respectively for the earliest to latest snapshot). The timescale was then determined by dividing the distance travelled by the velocity at which the gas is moving. These kinematic ages allowed for the outflow mechanical luminosity and mass loss rate to be determined. To calculate the outflow momentum and energy, the mass in each bin was multiplied by the velocity of that bin, and the momenta (P = M × V) and energy (E = 0.5M × V²) were summed over all channels with velocities greater than 2 km s⁻¹.

4. OBSERVATIONAL COMPARISONS

We report results from early times in the development of the collective outflow from a massive star forming region. In particular, none of the protostars in our model has yet reached a high enough mass, and ionizing luminosity, to form an ultracompact H ii region. There are no regions passing through this short-lived phase with distances less than a few kiloparsecs, close enough for their jets to be clearly resolved and imaged. However, at least two rather closer regions, Cepheus A and DR 21, appear to be in the immediately subsequent phase, with a major outflow associated with small, highly obscured ultracompact H ii regions. We therefore focus our comparison of our results to these two regions, after beginning with a general discussion of the outflow kinematics.

4.1. Kinematics of the Collective Outflow

All of the outflow properties derived from our simulated CO observations are shown in Table 1. The kinematic ages on which the calculation of the luminosity and mass-loss rates is based are also shown. The second snapshot appears younger than the first one because its intensity-weighted velocity grows disproportionally to its size. This can be explained by the fact that shortly before snapshot two, the third and fourth star in the group have formed (compare Figure 1), which already contribute to the collective outflow at its base but have not influenced its spatial extent yet. Table 1 also contains the real outflow ages, measured from the formation of the first sink particle. All measured quantities—the mass, momentum, kinetic energy, mechanical luminosity, and mass-loss rate—match well with the observed properties of high-mass outflows (e.g., Beuther et al. 2002; Wu et al. 2004).

Table 1. Outflow Properties Determined from Synthetic CO Observations

	Outflow Lobe	M	P	E	L	$\dot{M}$	t
	Outflow Lobe	(M_☉)	(M_☉ km s⁻¹)	(10⁴⁴ erg)	(L_☉)	(10⁻³ M_☉ yr⁻¹)	(kyr)
t = 0.624 Myr	Blue	13.0	136	179	12.4	1.09	12.0
	Red	11.7	147	232	16.1	0.99
t = 0.628 Myr	Blue	23.9	352	716	77.6	3.15	7.6
	Red	18.7	303	672	72.9	2.46

t = 0.635 Myr	Blue	30.4	487	1080	72.5	2.48	12.3
	Red	22.8	380	866	58.1	1.86

t = 0.642 Myr	Blue	16.4	224	438	20.7	0.94	17.5
	Red	23.2	391	893	42.2	1.33

Notes. Outflow mass M, momentum P, kinetic energy E, mechanical luminosity L, mass-loss rate $\dot{M}$ , and outflow age t as determined from the synthetic CO observations.

Download table as: ASCII Typeset image

We can compare the outflow kinematics measured from the simulated CO emission map to the actual model values. We determine the actual values of the mass, momentum and kinetic energy by summing over all cells above the disk plane with positive vertical velocity components, and similarly below the disk plane with negative vertical velocities. The resulting values are shown in Table 2. All masses and energies agree within a factor of two with the observationally determined numbers, while for the momentum the difference can be as much as a factor of three. This level of agreement is consistent with the results of a similar study for outflows in magnetohydrodynamical simulations (Peters et al. 2014). The kinematic and real ages also agree within a factor of two.

Table 2. Outflow Properties Determined from the Simulation Data

	Outflow Lobe	M	P	E	L	$\dot{M}$	t
	Outflow Lobe	(M_☉)	(M_☉ km s⁻¹)	(10⁴⁴ erg)	(L_☉)	(10⁻³ M_☉ yr⁻¹)	(kyr)
t = 0.624 Myr	Blue	15.3	45.5	252	20.7	1.53	10.0
	Red	14.4	48.4	338	27.8	1.44
t = 0.628 Myr	Blue	16.9	105	745	43.8	1.21	14.0
	Red	12.0	109	893	52.5	0.86

t = 0.635 Myr	Blue	19.8	216	2230	90.5	0.98	20.3
	Red	14.0	211	2260	91.7	0.69

t = 0.642 Myr	Blue	21.9	230	1660	50.1	0.80	27.3
	Red	14.5	198	2050	61.7	0.53

Notes. Outflow mass M, momentum P, kinetic energy E, mechanical luminosity L, mass-loss rate $\dot{M}$ , and outflow age t as measured from the simulation data.

Download table as: ASCII Typeset image

We performed ¹³CO and C¹⁸O ALMA simulations for two of the snapshots to determine the optical depth of the simulated observations of ¹²CO. We found that the emission is not very opaque. Accounting for the opacity of the ¹²CO line changed our derived outflow properties by less than a few percent.

We can compare these results with high resolution observations of the early stages of high-mass star formation. We note that statistical studies of these regions are not at high resolution, and further high-resolution observations are required for this type of analysis.

The energetics in the collective outflow are slightly lower than those observed in the group of outflows in 20293+3952 observed by Beuther et al. (2004), although the masses are close to those observed in this source, suggesting the velocities of the gas in the simulation are lower than those observed by Beuther et al. (2004). Conversely, we find masses higher than those observed by Palau et al. (2013) toward 19520+2759 MM1, but energies similar to their values and with similar dynamical timescales.

The mass, momenta and energetics are all two to three orders of magnitude larger than those seen by Arce et al. (2013) in HH 46/47, which is expected since the stars in our simulations are substantially more massive than those powering HH 46/47. What is interesting is that in Figure 17, we see multiple outflow components, just as they do (compare their Figure 9).

Qiu et al. (2009) calculated the outflow energetics for G240.31+0.07 and all of their values are approximately a factor of five greater than those determined here. Their region already shows evidence for an H ii region, and may be in a later evolutionary stage than the region being studied here. Similarly, Zapata et al. (2009) show the collimated outflow from a cluster of forming stars in W51N at a similar evolutionary stage to that of the Qiu et al. (2009) study.

These comparisons show that our outflow has higher masses and energetics than low-mass regions studied at similar resolutions, and lower values than those expected for more evolved clusters of forming stars: those that have already formed their H ii regions.

4.2. Cepheus A

We now turn to a detailed comparison with nearby objects. Cepheus A is a massive star-forming region lying at a distance of 700 +30/−28 pc, based on Very Long Baseline Array (VLBA) parallax measurements to two of the radio continuum sources it contains (Moscadelli et al. 2009; Dzib et al. 2011). This is consistent with previous estimates of at least 620 pc (de Zeeuw et al. 1999), and possibly as far as 730 pc (Johnson 1957; Crawford & Barnes 1970).

Cep A contains young stellar objects with a total bolometric luminosity of 2.5 × 10⁴ L_☉, exceeding that of a main-sequence B0 star (Evans et al. 1981; Ellis et al. 1990), lying behind a visual extinction of 90 mag. Hughes & Wouterloot (1984) detected multiple radio sources within the region, now denoted HW 1–9, consistent with 14 B3 stars with masses of around 8 M_☉ lying within a 0.1 pc radius. A rotating, dense core has been traced in CS emission with a radius of about 0.16 pc and a dynamical mass of 330 M_☉ (Narayanan & Walker 1996). The brightest source, HW 2, has about half the total luminosity, suggesting it is a B0.5 star approaching 15 M_☉ (Hughes et al. 1995). A Keplerian disk a hundred times smaller was resolved around HW 2 using CH₃CN line emission mapped with the Submillimeter Array (Patel et al. 2005). A jet whose knots have proper motions of order 500 km s⁻¹ has been observed in radio continuum within HW 2 (Curiel et al. 2006). Maser activity reveals at least two additional centers of star formation activity within 200 AU of HW 2 (Torrelles et al. 2001; Gallimore et al. 2003). A cluster of class I and class II YSOs also lies in this location (Gutermuth et al. 2005). The surrounding gas core traced by CS has a dynamical mass of 330 M_☉ within a radius of 0.32 pc (Narayanan & Walker 1996). Rodríguez et al. (1980) initially identified the major east–west bipolar CO outflow associated with Cep A. Additional components are aligned northeast-southwest (Bally & Lane 1991). CO radial velocities as high as 70 km s⁻¹ have been observed in the central 2' (Narayanan & Walker 1996), suggesting a powerful outflow has only recently become active.

Hiriart et al. (2004) resolved the emission from the H₂ v = 1–0 S(1) emission line with the spectral parameters used in our Figures 10 and 11. Their Figure 1 shows the same characteristic clumpiness over the same velocity scale as seen in those figures, including a hint of the highly structured leading bow shock seen in our models. Our edge-on Figure 10 is the better match to the observations, as Cep A lies nearly in the plane of the sky. This morphology is much better resolved in the line image shown in Figure 1 of Cunningham et al. (2009). This clearly shows the same multiple structures resembling bow shocks in both the interior and leading edge of the jet that are seen in our Figures 8 and 9. Cunningham et al. (2009) trace linear features through these structures and attribute them to a single, precessing object with arbitrarily chosen parameters, but our models self-consistently reproduce the same morphology with interacting, time-variable jets from the multiple sources forming during gravitational collapse of a core. Hiriart et al. (2004) also show a PV diagram in their Figure 4 generally consistent with our Figure 13. In particular, we reproduce the characteristic variation of the peak velocity along the length of the jet, due to intersections with interior shocks having high transverse velocity components.

4.3. DR 21

The distance to the next nearest massive star-forming region, DR 21 in the Cygnus X GMC complex (Schneider et al. 2006) has been determined to be 1.4 kpc by VLBA maser parallax measurements (Rygl et al. 2012), confirming the 1.5 kpc estimate of Odenwald & Schwartz (1993). Like Cepheus A, DR 21 contains a massive core with multiple ultracompact and hypercompact H ii regions (Harris 1973; Dickel et al. 1986; Cyganowski et al. 2003) linked to zero-age main-sequence O stars (Roelfsema et al. 1989); water, OH, and methanol masers (Plambeck & Menten 1990; Araya et al. 2009); many embedded IR sources; and a complex molecular outflow nearly in the plane of the sky (Davis et al. 2007). Zapata et al. (2013) showed multiple outflows traced by CO, SiO and methanol associated with the DR 21(OH) hot core roughly 2 pc away from the DR 21 H ii region, and probably part of the same star formation event.

High-resolution images of DR 21 with a narrow-band H₂ v = 1–0 S(1) emission line filter by Davis & Smith (1996), Cruz-González et al. (2007), and Davis et al. (2007) revealed morphology strongly reminiscent of that described above for Cep A, including multiple interior structures resembling bow shocks, and a complex, structured jet head, as shown in Figure 2 of Davis et al. (2007).

Davis & Smith (1996) measured line profiles with 20 km s⁻¹ velocity resolution, finding both asymmetric, blue-shifted profiles, and also more symmetric profiles (see their Figures 6 and 7). They often observed multicomponent structure in both sorts of profiles. In our Figure 14, we show line profiles from our run viewed at 90° (edge on) with 9.82 km s⁻¹ resolution. These include a clearly blue-shifted profile with multiple components, and symmetric profiles with and without multiple components. This agrees well with the somewhat lower spectral resolution Davis & Smith (1996) results. Cruz-González et al. (2007) published a PV diagram (their Figure 6), again showing the characteristic variation in the velocity direction associated with clumpiness in the emission seen in our Figure 13. Our model does not reproduce the large scale structure seen by Cruz-González et al. (2007), though, perhaps because of our choice of smooth initial conditions, which reduces the amount of large-scale density structure that the jet interacts with.

4.4. Summary

We have demonstrated that an outflow driven by multiple stars formed by fragmentation of the accretion flow onto a massive star appears consistent with observations of massive stellar outflows. We have done this in two ways. First, we compared the integrated mass and momentum of our simulated outflow to typical examples of massive stellar outflows at similar early stages in their evolution, finding that the simulated values lie within the observed ranges. We note that other mechanisms we have examined, such as magnetic driving (Peters et al. 2011) or ionization (Peters et al. 2012), fail this fundamental test.

Second, we compared our models to the relatively well resolved outflows from Cep A and DR 21. These share several characteristics. They have a chaotic structure with multiple apparent bowshocks, and relatively poor collimation. In particular, they do not resemble the well-defined structure of isolated, low-mass stellar jets, with a train of Herbig–Haro objects ending in a single bow shock, often connected by a well-defined high velocity beam. We argue that these multiple bow shocks can be naturally explained by Rayleigh–Taylor instability acting on the broad head of the jet as it propagates down the density gradient of the envelope; while the poor collimation occurs as the angular momenta of the sources gradually diverge from a common direction.

The kinematics of the massive outflows confirm these morphological impressions, with multiple velocity components leading some workers to characterize the flow as turbulent (Davis & Smith 1996; Hiriart et al. 2004; Cruz-González et al. 2007). Despite these chaotic tendencies, the jets do show general collimation, suggesting origination from stars formed from material sharing a common angular momentum axis, as would be expected in a set of stars formed from any single, gravitationally collapsing region, even in the presence of a turbulent background flow. All of these morphological and kinematic characteristics are also seen in the combined outflow generated by multiple, time-variable jets in our model of massive star formation.

We thank R. Banerjee for contributions to the development of the outflow subgrid model. We further acknowledge useful discussions with J. Bally, J. Mackey, and M. Krumholz, as well as assistance from J. C. Mottram. We thank the anonymous referee for comments that helped improve the presentation and discussion of this work. T.P. acknowledges financial support through SNF grant 200020_137896 and a Forschungskredit of the University of Zürich, grant no. FK-13-112. R.S.K. acknowledges funding from the Deutsche Forschungsgemeinschaft DFG via grants KL 1358/14-1 as part of the Priority Program SPP 1573 Physics of the Interstellar Medium as well as via the Collaborative Research Project SBB 811 The Milky Way System in subprojects B1, B2, B3, B4 and B5. M.-M.M.L. acknowledges partial support of this work from NSF grant AST11-09395. He also acknowledges hospitality during preparation of this paper by NORDITA, and by the Aspen Center for Physics, supported by their NSF grant PHY10-66293. C.F. acknowledges funding by a Discovery Projects Fellowship from the Australian Research Council (grant No. DP110102191). We acknowledge computing time at the Leibniz-Rechenzentrum (LRZ) in Garching under project ID h1343, at the Swiss National Supercomputing Centre (CSCS) under project IDs s364/s417 and at Jülich Supercomputing Centre under project ID HHD14. The FLASH code was developed in part by the DOE National Nuclear Security Administration ASC- and DOE Office of Science ASCR-supported Flash Center for Computational Science at the University of Chicago.

COLLECTIVE OUTFLOW FROM A SMALL MULTIPLE STELLAR SYSTEM

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

ABSTRACT

1. INTRODUCTION

2. NUMERICAL SIMULATION

3. ANALYSIS