EXPLOSIVE OUTFLOWS POWERED BY THE DECAY OF NON-HIERARCHICAL MULTIPLE SYSTEMS OF MASSIVE STARS: ORION BN/KL

Massive stars are relatively rare and tend to form in distant, obscured, crowded, and dynamic environments. The Orion Molecular Cloud 1 (OMC1) located immediately behind the Orion Nebula is the nearest site of ongoing massive star birth with a bolometric luminosity of about 10⁵ L_☉ (Gezari et al. 1998). OMC1 contains the Becklin–Neugebauer (BN) object (Becklin & Neugebauer 1967; Scoville et al. 1983), the first massive young stellar object to be discovered in the infrared and the Kleinman–Low (KL) thermal-infrared nebula (Kleinmann & Low 1967) which contains additional highly obscured massive stars. In addition to BN, radio sources I and n are thought to be massive because they are detected at radio wavelengths (Menten & Reid 1995). Recent submillimeter interferometry has revealed the presence of an additional protostellar candidate, SMA1 about 2'' south–southwest of radio source I (Beuther et al. 2006; Beuther & Nissen 2008). Although this is the brightest source in the submillimeter, it has not been detected at other wavelengths.

The BN/KL region contains a spectacular, wide opening-angle, arcminute-scale outflow (Allen & Burton 1993) traced by molecules such as CO and NH₃ that exhibit broad (>100 km s⁻¹) emission line wings (Kwan & Scoville 1976; Wiseman & Ho 1996; Furuya & Shinnaga 2009) and high-velocity OH, H₂O, and SiO maser emission (Genzel et al. 1981; Greenhill et al. 1998). This bipolar outflow with a southeast–northwest (redshifted–blueshifted) axis contains at least 8 M_☉ of accelerated gas with a median velocity of about 20 km s⁻¹. Interferometric CO images, H₂O masers, and dense-gas tracers such as SiO reveal a smaller (8'' long) and younger (∼200 year old) outflow along a northeast–southwest axis emerging from radio source I orthogonal to the arcminute-scale CO outflow (Beuther & Nissen 2008; Plambeck et al. 2009). The momentum and kinetic energy content of these flows is at least 160 M_☉ km s⁻¹ and 4 × 10⁴⁶ erg (Snell et al. 1984) to 4 × 10⁴⁷ erg (Kwan & Scoville 1976). Zapata et al. (2009) presented a CO J = 2–1 interferometric study and found a dynamic age of about 500 years for the larger OMC1 outflow. They noted its impulsive nature, that its structure is different from accretion-disk powered flows, and that it originated several arcseconds north of the OMC1 hot core.

While most protostellar outflows appear to be driven by jets or collimated winds (Reipurth & Bally 2001; Beuther & Shepherd 2005; Arce et al. 2007; Qiu et al. 2008), the BN/KL outflow is a wide-angle explosion. Near-infrared images reveal hundreds of individual bow shocks in the lines of [Fe ii] and H₂ that emerge in nearly every direction (Allen & Burton 1993; Kaifu et al. 2000). The higher excitation species such as the 1.64 μm line of [Fe ii] trace the bow-shock tips (dubbed "bullets") while lower excitation tracers such as H₂ trace the wakes or bow-shock skirts (dubbed "fingers"). A few of these shocks protrude from the molecular cloud and can be seen at visual wavelengths in [O i], [S ii], [N ii], and Hα (Axon & Taylor 1984; O'Dell et al. 1997). The highest velocity shock, HH 210 (∼400 km s⁻¹) has been detected in X-rays (Grosso et al. 2006). Proper motions indicate a common age of about 10³ yr or less for many features (Doi et al. 2002), comparable to the 500 year age estimate of Zapata et al. (2009).

The three brightest radio-emitting stars in OMC1, sources BN, I, and n, have proper motions (motions in the plane of the sky) of 25, 15, and 26 km s⁻¹ away from a region less than 500 AU in diameter from which they were ejected about 500 years ago (Rodriguez et al. 2005; Gómez et al. 2005, 2008). Apparently, a non-hierarchical multiple star system containing at least four members experienced a dynamical interaction resulting either in the formation of a tight binary or possibly even a stellar merger (Bally & Zinnecker 2005; Bally 2008; Zapata et al 2009) whose (negative) gravitational binding energy ejected these stars from the OMC1 core. With estimated stellar masses of 8–13, 12–20, and 2–3 M_☉ for BN, I, and n, respectively, the kinetic energy of the ejected stars is 1 to 2 × 10⁴⁷ erg (Rodriguez et al. 2005; Gómez et al. 2005, 2008). Thus, the total kinetic energy of the ejected stars and the outflow is about 1.5 to 5.5 × 10⁴⁷ erg. Assuming that the outflow has suffered some dissipation of its original kinetic energy as it interacted with its surroundings, a reasonable range of initial energies required to eject the stars and outflow is 2 to 6 × 10⁴⁷ erg.

In this paper, new proper-motion measurements of H₂ are presented and combined with published data. The results indicate a 500 year age for the fastest ejecta in the OMC1 bullet and finger system and show that the point of origin of the flow coincides with the location of the ejected stars 500 years ago to within the measurement errors. Motivated by an apparent common event that ejected both the stars and the outflow, we present a plausible scenario for the recent evolution of the OMC1 region.

New images of the BN/KL outflow were obtained on 2004 November 21 with the Near-Infrared Camera and Fabry–Perot Spectrometer (NICFPS) on the ARC 3.5 m telescope at Apache Point Observatory in New Mexico. NICFPS uses a 1024 × 1024 pixel Rockwell Hawaii 1-RG HgCdTe detector. The pixel scale of this instrument is 0273 per pixel with a field of view 458 on a side. Images with 300 s exposures were obtained in the 2.122 μm S(1) line of H₂ using a narrowband filter (FWHM = 0.4% of the central wavelength). Separate sky frames were obtained using the same exposure time at a location 600'' east. A set of five dithered images were obtained both on-source and on the sky position. A median-combined set of unregistered, mode-subtracted sky frames were used to form a master sky frame that was subtracted from individual images. The reduced images were corrected for optical distortions. Field stars were used to align the frames which were median-combined to produce the image shown in Figure 1. Atmospheric seeing produced 08 full-width at half-maximum stellar images. A full description of the data acquisition, the data reduction, and analysis procedures is given in Cunningham (2006).

Zoom In Zoom Out Reset image size

Proper motions for individual features in the BN/KL outflow were determined by comparison of the 2004 epoch NICFPS image with H₂ images taken on 1992 September 13 using the IRIS instrument at the Anglo-Australian Telescope (AAT) described by Allen & Burton (1993), and the narrowband images taken on 1999 January 11 and 13 with the CISCO instrument on the Subaru telescope on Mauna Kea, Hawai'i (Kaifu et al. 2000). These images were registered to the NICFPS image using field stars to an accuracy of about 005, limited by optical distortions. The time intervals between the first image and the second and third images are 6.33 and 12.19 years.

Proper motions were determined using difference-squared cross-correlation technique of Currie et al. (1996). All images were re-sampled to the pixel scale of the 2004 epoch NICFPS image and blinked to identify moving features. Measurement boxes containing between 50 and 500 pixels were defined on the 2004 epoch image as shown in Figure 1. On the prior-epoch images, a larger measurement box with six additional pixels beyond each boundary was defined at each location. Slowly varying background emission was removed by fitting a two-dimensional cubic function and subtraction of the result. The smaller box is shifted within the larger box for all possible x and y pixel offsets and subtracted from the larger box. The square of the sum of all pixel-to-pixel differences is assigned to the pixel that corresponds to the x and y offset in an "offsets" image consisting of all possible pixel shifts. The location of the minimum in this two-dimensional array is an estimator of the x and y offset that results in the best match between the intensities in the boxes defined on the first and second epoch images. Proper motions are estimated using the centroid of a two-dimensional Gaussian fit to this "offsets" image. Field stars were used to estimate a mean proper-motion error of about 15 km s⁻¹ for bright compact sources for both the 2004–1999 and 2004–1992 comparison.

Most proper-motion vectors point directly away from the OMC1 cloud core. However, about 10% of the vectors differ by more than 30° from the expected direction of motion away from the point of origin of the outflow. The four boxes located in the lower-right portion of Figure 1 trace a highly collimated chain of H₂ bow shocks that emerge at position angle P.A. ∼320° from the OMC1-S core located about 90'' south OMC1. Another collimated outflow crosses the northwestern lobe of the OMC1 outflow at position angle P.A. ∼240°–250°. Its brightest southwest-facing bow shock is located at J(2000) = 5^h35^m116, −05°20'54''. The orientation of this shock indicates that the source is likely to be located toward the northeast in the "Integral Shaped Filament" which contains the highest column density of molecular gas (Johnstone & Bally 1999).

Figure 1 shows the H₂ proper-motion vectors. The fastest motions are located at the greatest distance from the OMC1 core while the slower motions are located closer in. Thus, to first order, the motions are consistent with velocity being proportional to distance (e.g., a "Hubble flow"), an indication of an explosion.

Figure 2 shows the proper motions of 2.122 μm H₂ features measured by Cunningham (2006, and this work), 1.644 μm [Fe ii] (Lee & Burton 2000), and Herbig–Haro objects traced by Hα, [O i], [N ii], and/or [S ii] (Doi et al. 2002) as a function of distance from the point of origin of the three high-velocity stars. H₂ tends to only trace features with relatively low proper motions below about 70 km s⁻¹. Presumably, this is a consequence of H₂ dissociation in magnetized C shocks with speeds in excess of about 70 km s⁻¹. For non-magnetized J shocks, H₂ dissociates at shock speeds of 10–30 km s⁻¹. The pair of ∼250 km s⁻¹ H₂ features located about 0.1 pc from the expansion center in Figure 2 may trace reverse shocks moving into dense molecular clumps that have high velocities. Hα excitation requires shock speeds of at least 60 km s⁻¹ while the λ 5007 A line of [O iii] requires shocks with speed in excess of at least 150 km s⁻¹. Most visual-wavelength-forbidden transitions such as the λ 6717/6731 [S ii] lines trace the cooling layers behind fast shocks. [Fe ii] and the visual-wavelength tracers are commonly associated with features having velocities up to 425 km s⁻¹.

Zoom In Zoom Out Reset image size

Figure 2 demonstrates that some proper-motion features have apparent ages between 500 and 1000 years. Dynamical ages are determined from the distance traveled divided by the magnitude of the proper motion. Unlike age determination from radial velocities, knowledge of the inclination of the flow vector is not required. For example, HH 210, the fastest HH object associated with the BN/KL outflow, has a dynamic age of 605 years assuming no deceleration. Given that most shocks are likely to have been decelerated by interactions with their environment, Figure 2 shows that the ages of many features in the BN/KL outflow are consistent with ejection about 500 years ago, similar to the time elapsed since the ejection of the three high-velocity stars, BN, I, and n.

Figure 3 shows the apparent point of origin of the H₂ proper motions in the OMC1 outflow. The left panel is made by backward tracing of the proper-motion vectors, effectively marking their trajectories in the absence of deflections. The peak in the density of intersecting trajectories provides the best estimator for the point of origin of the measured motions. The "intensity" of each trail is weighted by the speed. Thus, large motions, which can be measured more precisely than small motions, have a greater weight. The image in the right panel has been boxcar-smoothed by 31 pixels (25'') to aid visualization of the intersection point. The point of origin of the proper motions is located at J(2000) = 5^h35^m145, −5°22'23'', with a 1σ uncertainty of about 5''. This location is offset by about 5'' to the northeast of the apparent ejection center of the massive stars and 9'' north of the current location of radio source I. Zapata et al. (2010) found an expansion center based on interferometric CO J = 2–1 observations at J(2000) = 5^h35^m1437, −5°22'279, about 5'' due south the expansion center derived from our proper motions, but consistent within the uncertainty estimates.

Zoom In Zoom Out Reset image size

High-velocity runaway stars are relatively common among massive O and early B stars but rare among lower-mass stars (Gies & Bolton 1986; Gies 1987). While approximately 30% of O stars are runaways with velocities greater than 30 km s⁻¹ with respect to their birth clusters and associations, less than 2% of B and A stars have such velocities (Stone 1991; Schilbach & Roser 2008; Zinnecker & Yorke 2007). Many (possibly most) massive runaway stars were ejected by three-body encounters between a single star and a massive binary (Gvaramadze & Gualandris 2011), four-body encounters between a pair of massive binaries (Poveda et al. 1967), or the decay of non-hierarchical multiple star systems (van Albada 1968; Valtonen & Mikkola 1991; Aarseth et al. 1994; Sterzik & Durisen 1995, 1998). Such three- or four-body encounters are most likely to occur in very dense, young star clusters. The runaway stars AE Aur and μ Col, ejected with speeds of over 100 km s⁻¹ about 2.5 × 10⁶ yr ago (Hoogerwerf et al. 2001) from the NGC 1980 cluster located about 30' south of the Orion Nebula, provide a prototypical example. This event produced the massive colliding-wind X-ray binary ι Ori (Gualandris et al. 2004).

Non-hierarchical multiple star systems, triple or higher-order multiples in which the time-averaged distances between components are similar are dynamically unstable. Numerical experiments show that non-hierarchical systems tend to rearrange into hierarchical configurations consisting of a compact binary and either ejected stars or ones in elongated orbits with a much larger semi-major axis than the binary. Non-hierarchical systems tend to decay in about 100 crossing times; the crossing time is roughly given (in years) by t_cr ∼ 0.17(R³/M)^1/2 where R is the characteristic length scale of the system in AU and M is total system mass in solar masses (Anosova 1986; Sterzik & Durisen 1995). Many low-mass pre-main-sequence stars are born in such non-hierarchical multiples (Reipurth 2000). Among massive stars, Trapezium-like groups tend to be non-hierarchical and subject to dynamical decay.

Both the ejection of massive stars and the launch of the OMC1 BN/KL outflow may be the result of the dynamic interaction and rearrangement of a system of massive stars in OMC1 about 500 years ago. The kinetic energy originates from the release of gravitational potential energy accompanying the dynamic formation of a compact binary. According to the virial theorem, one-half of the (negative) potential energy of the binary goes into orbital motion of its constituent stars. The other half is removed from the system and is available to accelerate the binary, other stars involved in the interaction, and surrounding gas. The most massive component of the system, which is the slowest moving object, is likely to be the binary. The binding energy of a binary consisting of masses m₁ and m₂ with an orbital semi-major axis a, E_b = −Gm₁m₂/2a, must be twice the (negative) sum of the observed kinetic energies in the outflow and ejected stars.

Assuming that source I consists of two 10 M_☉ stars, to release 2 × 10⁴⁷ (or 6 × 10⁴⁷) erg of free kinetic energy, its members must have a semi-major axis less than 2.2 (or 0.74) AU with an orbital period shorter than 2 (or 0.4) years. The perihelion velocity of the stars must be at least 33 (or 55) km s⁻¹ and the binding energy of the binary must be less than −4 × 10⁴⁷ (or −12 × 10⁴⁷) erg. If either BN or n is the binary, the binary separation must be much smaller to provide the energy for the ejection of the three stars plus outflow. Stellar ejection removed a total of about 22–36 M_☉ of stars and least 8 M_☉ of gas. Thus, this event removed 30–50 M_☉ from the OMC1 core. The point of origin of the ejected stars is located at the center of the BN/KL outflow about 4'' (∼1700 AU in projection) northwest of the present location of radio source I.

Momentum conservation requires that the vector sum of the stellar and outflow momenta sum to zero in the rest frame of the parent cloud. Thus, the runaway single stars and the dynamically formed binary are most likely to be ejected from their star-forming clump. Removal of the stars implies that the envelope whose motion was dominated by the gravity of the stellar system before ejection can become unbound after ejection. Multi-body encounters tend to eject the outer portions of pre-existing circumstellar disks at radii larger than about 1/3–1/2 of the periastron separation. Thus, the outer parts of disks around the stars that form the binary and the envelope surrounding the disrupted cluster will be ejected with roughly their respective orbit speeds.

Three distinct mechanisms may be involved in releasing gas from the gravitational potential well. First, the disruption and ejection of the outer parts of circumstellar disks during the final dynamic encounter leading the binary formation. Second, recoil of the pre-existing circumcluster envelope following stellar ejection. Finally, the release of stored magnetic stress. Each process is a consequence of the removal of gravitational attraction on a timescale given by the relevant crossing time. Disk disruption occurs within a timescale t_cross ≈ r/V≪ 1 year at r = 1 AU for V∼ 100 km s⁻¹, the Kepler speed around a 10 M_☉ star. At r = 100 AU, stellar ejection with a mean speed of V = 20 km s⁻¹ implies that the envelope gas would be unbound several decades. In both cases, the expulsion velocity would be comparable to the Kepler speed prior to ejection, or about 20 km s⁻¹ for an envelope initially orbiting 40 M_☉ of stars at r = 100 AU to over 100 km s⁻¹ for parts of disks ejected from within 1 AU of a 10 M_☉ star. These mechanisms are explored in more detail below.

3.1. Formation of a Non-hierarchical Massive Multiple in OMC1

OMC1 spawned at least 20 to over 40 M_☉ of stars. Assuming a typical star formation efficiency of 25%–50% typical of massive cluster forming clumps, the OMC1 proto-cluster clump must have initially contained at least M_clump ∼ 100 M_☉ of gas. As a fiducial example, consider a clump with an "isothermal" density distribution ρ(r) = ρ₀(r/r₀)⁻² for radius r less than the clump outer radius R. A clump with mass M_clump = 100 M_☉ inside R = 0.1 pc implies mean orbital velocity (which is independent of r), V ≈ (GM_clump/R)^1/2 = (4πGρ₀)^1/2 = 2.1 km s⁻¹, where ρ₀ ≈ 5.0 × 10¹⁶ g cm⁻³ would be the density at r₀ = 1 cm if the power-law distribution were extended to this radius. In this example, the density at r = 5000 AU from the center is n(H₂) ≈ 2 × 10⁶ cm⁻³ and the crossing time is t_cross ∼ R/V < 5 × 10⁴ yr. The most massive protostellar objects are likely to either form in the center where the density and pressure are the greatest, or if they form away from the center, to rapidly migrate to the center of the potential due to orbit decay (Zinnecker 1982; Bonnell et al. 2001).

Orbits can decay due to Bondi–Hoyle (BH) accretion as they move though clump gas and dynamic friction caused by the gravitational pull of the wake produced by converging streamlines behind the YSO. The BH accretion rate is roughly , where C_s is the effective sound speed (which includes the effects of turbulent motion and magnetic support) in the gas (usually <1 km s⁻¹). A star moving with a speed V = 2.1 km s⁻¹ 5000 AU from the center of the above gas distribution would accrete at a rate  M_☉ yr⁻¹ for a 1 M_☉ YSO and  M_☉ yr⁻¹ for a 5 M_☉ YSO. As illustrated by the simulation described in the next paragraph, the most massive objects grow most rapidly, doubling their mass and migrating to the center of the cloud within few × 10⁵ yr. Fragmentation (McKee & Tan 2002, 2003; Krumholz 2006) or the capture of sibling stars (Moeckel & Bally 2006, 2007a, 2007b) may produce multiple systems. These processes could have resulted in the assembly of a sub-cluster of massive stars and binaries in the center of the core.

Figure 4 shows a sample result of a simplified numerical model illustrating the combined effects of BH accretion onto protostellar seeds, orbit decay, formation of a non-hierarchical system of massive stars, and dynamical ejection. This simulation starts with a 50 M_☉ clump of gas in a Plummer potential given by , where M is the total mass and a = 0.025 pc is the core radius. (This potential is used to avoid the singularity at the center of an isothermal sphere.) The clump is non-rotating with an internal velocity dispersion (effective sound speed) given by the local virial velocity so that, to first order, the clump is stable to global gravitational collapse. At the start of the simulation, 20 protostellar seeds, each having an initial mass of 0.5 M_☉, are distributed randomly in the clump. Their velocity distribution is virialized in the combined gravitational potential of the stars and gas and their velocity-vector orientations are random. The simulation was run 100 times, each with a different randomly chosen set of initial locations and velocities for the seeds.

Zoom In Zoom Out Reset image size

Stellar motions are calculated from the gravitational forces exerted by other stars and the static gas cloud, modified by the effects of BH accretion. The inter-stellar forces are directly calculated, and the orbits are integrated with a global but variable time step using a (7,8) order Runge–Kutta pair (Prince & Dormand 1981). The code is a standard n-body integrator with the addition of a static potential, identified with the natal gas, that can be accreted onto protostellar seeds. The code follows the stellar mass growth and orbit evolution in the combined gravitational potential of the clump plus embedded stars. Similar approaches have been used by Moeckel & Clarke (2010) and Baumgardt & Klessen (2010). The protostellar seeds experience BH accretion from the clump and grow in mass with the accretion radius set to the minimum of the BH radius or one-third the distance to the closest neighboring star. To keep the computations simple, the response of the gas to the passage of the stars is not modeled. This simplification is reasonable because the zone influenced by the passage of a protostellar seed is restricted to its gravitational radius, r_G = G_*M/V²_* (∼200 AU for M_* = 1 M_☉ and V_* = 2 km s⁻¹). Random motions in the clump will tend to fill in the cavity formed by the passage of a star on a timescale τ ∼ r_G/C_s⩾ 10³ yr but much less than the 10⁵ yr evolution timescale of the protostars and clump. Thus, it is assumed that on average, the envelope remains spherical and smooth. However, as gas is accreted onto the stars the mass of the gas cloud, and thus the contribution of the gas to the gravitational potential, is reduced accordingly.

Protostars located in the densest portion of the clump and/or those that by chance have the lowest relative velocity with respect to the gas tend to grow in mass most rapidly. In this simulation, accreted gas lowers the orbital angular momentum of the stars about the clump center since the clump is assumed to be non-rotating. Thus, as they grow in mass, their orbits decay and the most massive objects migrate to the clump center. The most rapidly growing protostars form a non-hierarchical system near the center. As they come to dominate the mass in this region, mutual gravitational interactions reconfigure the system into a hierarchy consisting of one or more tightly bound compact binaries plus ejected high-velocity stars. Multiple runs of the simulation with different randomly chosen initial locations and velocities almost always produce mass-segregated, non-hierarchical systems of massive stars at their centers. This model illustrates one plausible scenario for the formation of a non-hierarchical system of massive stars and their subsequent ejection. Three essential features are missing from this simplified model; the dynamic ejection of gas from the central region by the orbital motion of the massive stars, the presence of accretion disks, and the streamers that transport mass from the inner boundary of the envelope onto these disks. Previous simulations of binary star formation have shown these features to be present (e.g., Lubow & Artymowicz 1996).

In a multiple star system, matter that enters the region where the stars orbit each other tends to be expelled by gravitational torques. As a result, gas bound to the system is likely to be organized into two components: an extended outer envelope having an inner boundary with a radius somewhat larger than the semi-major axis of the largest stellar orbit, and inner disks surrounding the stars with outer radii several times smaller than the periastron separations (Artymowicz & Lubow 1994; Günther & Kley 2002). For pre-decay interstellar separations of 100 AU, matter located within ∼300 AU of the cluster would either be accreted onto disks with outer radii less than about 30 AU or be expelled to beyond 300 AU. Gas in an envelope bound to a 40 M_☉ cluster of stars with an inner radius of 300 AU would have Kepler speeds less than 11 km s⁻¹. The outer boundary where the cluster has significant influence on the envelope can be defined by the radius where the gravitational influence of the cluster falls below the velocity dispersion, about 2–3 km s⁻¹ for OMC1. Disks with 30 AU outer radii orbiting 10 M_☉ stars would have Kepler speeds greater than 17 km s⁻¹. Since non-hierarchical multiples likely have chaotic orbits, disks may be truncated at smaller radii. Gas falling from the outer envelope onto individual circumstellar disks may form transient streams in the region of avoidance.

Figure 5 illustrates a possible scenario for the evolution of the OMC1 cloud core. Figure 5(a) shows the formation of several protostellar seeds destined to become massive stars as they continue to accrete from the clump and sink to its center to form a non-hierarchical multiple system. Figure 5(b) shows the multiple system after it has cleared the orbit zone. Tidal streams can continue to feed the growth of compact circumstellar disks with outer radii smaller than the typical periastron separation. Figure 5(c) shows the aftermath of the dynamic interaction that ejected the stars from the clump.

Zoom In Zoom Out Reset image size

3.2. Recoil of the Envelope: Slow Ejecta

3.3. Disruption of Disks: Fast Ejecta

3.4. The Contribution of Stored Magnetic Stress

The birth and subsequent orbital motion of a non-hierarchical cluster of four or five massive stars in OMC1 over 10⁵ or 10⁶ yr may have led to the amplification of ambient magnetic fields in both the envelope surrounding the cluster and in individual circumstellar disks. Fields can grow until either magnetic reconnection or resistive dissipation in the charged component of the weakly ionized plasma come to equilibrium with the dynamo (Lynden-Bell 1996; Wheeler et al. 2002; Vasil & Brummell 2008, 2009). The orbital motion of the protostars in the cluster gravitational potential and the orbital motion of gas in individual circumstellar disks may each contribute to the amplification of the magnetic fields (Käpylä et al. 2008).

When equipartition is reached, B² ≈ 8πρ(r)c²_s, where ρ(r) is the local gas density and c_s is the effective sound speed including turbulent motions. For the envelope having the power-law density profile described in Appendix A, the equipartition magnetic energy density at radius r is given by _Be = B²(r)/8π = ρ(r)c²_s where c_s is the local effective sound speed (including thermal and turbulent motion) and the subscript e refers to the envelope. Under this assumption, the magnetic energy in a shell at radius r and thickness dr is

which can be integrated from the inner edge of the envelope to a radius r to give

Setting r = r_out gives E_Be = M_ec²_s. For M_e = 60 M_☉ and c_s = 2 km s⁻¹ this results in E_Be ≈ 5 × 10⁴⁵ erg for the equipartition magnetic energy that could be stored in the envelope.

If magnetic stress reaches equipartition with the effective sound speed within circumstellar disks, a similar calculation yields E_Bd = M_dc²_d where the subscript d refers to the disk quantities. The effective sound speeds in the disks are likely to be considerably larger due to proximity to the massive stars. Using c_s = 10 km s⁻¹ gives E_Bd = 2.0 × 10⁴⁵M₁c²₁₀ erg for the equipartition magnetic energy that could be stored in a disk. Here, M₁ is in units of 1 M_☉ and c₁₀ is in units of 10 km s⁻¹.

Magnetic fields may become stronger than equipartition with the local pressure in the disk or envelope. Keplerian shear can wind up the field so that its energy approaches the gravitational potential, a circumstance discussed in the context of circumstellar disks by Shu et al. (2008). The condition that Alfvén speed (or effective sound speed) is equal to the gravitational escape speed sets an upper bound on the field strength. Thus, V²_A = B²/4πρ ≈ GM/r, where M is the mass enclosed inside radius r. In this case, the field could have stored more than 10⁴⁷ erg of magnetic energy in the volume inhabited by the cluster. The removal of about 40 M_☉ of stars would have left the circumcluster gas severely magnetically overpressured. The release of magnetic energy, in effect a "magnetic bomb" (Matt et al. 2004, 2006), may have been a major source of energy for launching the OMC1 outflow. The strongest fields are likely to be generated by shear dynamos in the innermost parts of circumstellar disks surrounding individual massive stars. Future Zeeman measurements of clumps containing sub-clusters of massive protostars using the 24–26 μm iron lines with SOFIA may determine if energetically significant, gauss-strength fields can be generated in such environments (R. Crutcher 2010, private communication).

The BN/KL outflow is velocity-segregated with relatively slow gas confined to a few tens of arcseconds from radio source I and the fastest ejecta, located near the tips of the [Fe ii] and H₂ fingers about 120'' northwest and southeast of the core region. Proper-motion measurements of visual-wavelength HH objects (such as HH 205 through 210), near-infrared emission from shock-excited [Fe ii] and H₂ (Figures 1 and 2), and the interferometric CO J = 2–1 study of Zapata et al. (2010) indicate that the BN/KL outflow was generated by an explosive event about 500 years ago. This age coincides with the dynamical re-arrangement of the massive stars in OMC1 during which radio sources I, n, and BN were ejected (Gómez et al. 2005, 2008). The slowest star, radio source I, is probably a compact binary consisting of roughly 10 M_☉ stars with a separation less than about 2 AU. Future long-baseline radio interferometry, or precision radial velocity measurements in the infrared, may determine if this conjecture is true. It is proposed that the formation of this binary and the consequent release of 2 to 6 × 10⁴⁷ erg of gravitational potential energy powered both the motion of the stars and the supersonic expulsion of gas in the BN/KL outflow. The slowest ejecta may have been produced at large radii as the pre-ejection orbits of clumps in the envelope became unbound following stellar ejection. The fastest ejecta may have originated from near the center of the potential well vacated by the stars. This situation is similar to that postulated by McCaughrean & Mac Low (1997) and Stone et al. (1995), where the H₂ fingers are generated by the interaction between a slow, outer wind and a fast, later ejected flow. The estimates based on simple models of envelopes and circumstellar disks with power-law density distributions indicate that the gravitational energy stored in orbital motion prior to ejection can explain the energetics of the outflow following stellar ejection. However, numerical modeling of realistic configurations is needed to determine if this hypothesis is valid.

The explosive OMC1 outflow consists of multiple fingers and wakes of H₂ emission produced by the fastest ejecta plowing through the OMC1 core. The volume-filling factor of the cavities created by this ejecta is likely to be lower than the enclosed volume, especially toward the southeast. Observations show that a reservoir of dense gas is located in the "hot core" 5''–20'' southeast of the site of stellar ejection. Zapata et al. (2010) propose that the "hot core" may be blocking portions of the outflow in this direction. The bulk of the H₂ emission is located to the northwest, consistent with more of the original core being blown-out in this direction. The ejected stars moving south likely plow through the remaining dense material.

4.1. Post-ejection Disks and Outflow Orientations

4.2. Disks Accreted from the Surrounding Medium?

In the proposed scenario, the currently observed disks either trace (1) material that has expanded from within a few AU of each star, (2) disk debris ejected at less than escape speed that fell back to re-form a new disk in the last 500 years, or (3) was captured by BH accretion as the ejected stars plowed through remaining dense gas such as the "hot core" in OMC1.

A star with mass M_* moving with velocity V through a cloud can accrete material lying within a cylinder defined by the gravitational radius r_G = GM_*/V² which is measured orthogonal to the stellar trajectory. After traversing a distance L the amount of mass captured is about

where μ is the mean molecular weight, m_H is the mass of hydrogen, n₅ is the number density of H₂ in units of 10⁵ cm⁻³, L₃ is the path traversed by the star in units of 10³ AU, M₁₀ is the mass of the star in units of 10 M_☉, and V₁₀ is the velocity of the star through the medium in units of 10 km s⁻¹. Throop & Bally (2008) explored some consequences of BH accretion onto low-mass protostars as they moved through dense clouds. Moeckel & Throop (2009) modeled this late-phase BH accretion from a uniform density cloud onto a young star surrounded by an accretion disk using a hydrodynamic code with an isothermal equation of state but did not consider velocity or density gradients.

BH accretion tends to be subject to a "flip-flop" instability that leads to highly variable accretion rates and fluctuating orientations of the angular momentum vector of the flow (Shima et al. 1998; Fryxell & Taam 1988; Ruffert 1997, 1999; Krumholz et al. 2006). However, gradients in either velocity or density orthogonal to the stellar velocity give the flow vorticity with respect to the accretor. The resulting angular momentum can curtail accretion directly onto the moving mass because the flow tends to form a disk-like structure. Ruffert (1997) presented three-dimensional simulations of BH accretion from a medium having a velocity gradient while Ruffert (1999) considered flows from media with density gradients. Krumholz et al. (2005) modeled accretion from a medium having vorticity. These studies show that BH accretion from a medium with a small velocity or density gradient is still subject to the "flip-flop" instability. However, large gradients tend to produce accretion disks with a relatively stable angular momentum vector orthogonal to the star's velocity (Ruffert 1999).

The BN object is moving away from the hot core toward the northwest and has traversed about 4000 AU since ejection. If the average density of the medium through which it moved had n(H₂)∼ 10⁵ cm⁻³, a 30 km s⁻¹ velocity implies that it could have swept up about M_BH = 1.4 × 10⁻⁶ M_☉. If source n has a mass of 5 M_☉, and traversed 2000 AU through a medium with a density of 10⁵ cm⁻³ with a velocity of 30 km s⁻¹, it could have BH-accreted M_BH = 1.2 × 10⁻⁷ M_☉.

Radio source I is moving toward the dense "hot core." Assuming that it has a mass 20 M_☉, is moving through a medium with n(H₂) = 10⁶ cm⁻³ with a velocity of 10 km s⁻¹, and has traversed a distance on order 1000 AU, it could have swept up as much as M_BH = 8 × 10⁻⁴ M_☉. For source I, assuming a mass M_I = 20 M_☉, the time to cross the gravitational radius, GM_I/V²_I∼ 100 AU at velocity V_I = 13 km s⁻¹ is less than 40 years. The Kepler orbit time at a distance of 50 AU from a 20 M_☉ compact binary is about 80 years. Thus, any material accreted within the last several hundred years would have orbited at least three times at 50 AU. Thus, enough time has elapsed for accreted matter to have damped much of its motion orthogonal to its plane of symmetry and would resemble a disk.

Both BN and radio source I are moving close to the plane of the sky. If the accreted disk's angular momentum vector lies in a plane orthogonal to the stellar velocity vector, the apparent major axis of such a disk as projected onto the plane of the sky will be aligned with the proper-motion vector. Observations show that the major axes of both the 7 mm continuum emission (Reid et al. 2007) and the SiO maser emission (Goddi et al. 2009) from radio source I (which are interpreted as a nearly edge-on disk), and the putative disk associated with BN (Jiang et al. 2005) are parallel to their proper motions. These disk orientations are consistent with BH accretion onto the ejected stars from a medium having a significant density or velocity gradient orthogonal to the stellar motion.

4.3. Is the Orion BN/KL Outflow Unique?

The main conclusions of this study are the following.

• 1.  
  New proper-motion measurements show that some of the fastest and most distant H₂ knots in the OMC1 BN/KL outflow have dynamical ages of 500 years, similar to the 500 year interval since dynamical decay ejected radio sources I, BN, and n from a region located only a few arcseconds from the center of the outflow. Thus, it is likely that the outflow was produced by the same event that ejected the massive stars BN, I, and n from the core.
• 2.  
  To within an order of magnitude, the kinetic energies of the outflow and ejected stars are comparable.
• 3.  
  A scenario is explored in which the explosive OMC1 outflow is produced by the dynamical ejection of massive stars. Within the last few hundred thousand years, a compact non-hierarchical multiple system of massive stars or a pair of massive binaries must have formed in the OMC1 cloud core. A highly idealized numerical model is presented in which protostellar seeds grow by BH accretion from the parent clump. Seeds which by chance wander into the densest part of the core or have the lowest relative velocity with respect to the surrounding gas experience the highest accretion rates. They become massive, experience orbit decay, and migrate to the center of the model where they form a dynamically unstable multiple system.
• 4.  
  The kinetic energy of the stars and the outflow could have been generated by release of gravitational binding energy accompanying the formation of a compact binary, most likely radio source I. If source I consists of a pair of 10 M_☉ stars, then to produce the observed kinetic energy of the ejected stars and the outflow, between 2 and 6 × 10⁴⁶ erg, the mean separation of the binary must be between 2.2 and 0.7 AU.
• 5.  
  In the proposed scenario, the outflow is driven by the energy and momentum stored in orbital motion prior to dynamical decay. Previous studies of gas flows in multiple-star systems show that orbiting material tends to be organized in a hierarchy consisting of tightly bound circumstellar disks with outer radii less than about 1/3 of the periastron separation and relatively loosely bound envelopes with inner radii several times the apastron separation. Such a scenario is envisaged to have been present in OMC1 prior to the ejection of its massive stars.
• 6.  
  Stellar ejection may contribute to the formation of an explosive outflow by three mechanisms. Recoil of the envelope: the removal of stellar mass from the center of the envelope would result in the conversion of its orbital motion into linear motion. Disruption and ejection of circumstellar disks: the final penetrating encounter that formed a compact binary would eject the inner parts of pre-existing disks to produce the fastest ejecta. Release of magnetic stress: magnetic energy potentially produced by shear-dynamo action in both circumstellar disks and the envelope might boost the velocities of the ejecta. In this scenario, the fastest ejecta is launched first from deep inside the gravitational potential of the decaying cluster. This material must plow through the slower moving and later ejected envelope. Such interactions are prone to Rayleigh–Taylor-type instabilities as shown by models of fast winds slamming into slower, previously ejected winds. Order-of-magnitude energy estimates show that all three mechanisms are plausible, and that all three may contribute to the energy budget.
• 7.  
  The ejected high-velocity stars may have accreted new circumstellar material as they traversed the dense OMC1 core. In the presence of gradients in the medium, the angular momentum vectors of accreting gas tend to be orthogonal to the stellar velocity vectors. Therefore, for stellar motions close to the plane of the sky, disk major axes will be aligned (in projection) with the stellar proper-motion vectors, as observed for both radio source I and BN. Consequently, any recent outflow activity produced by the ejected stars will have axes orthogonal to the stellar proper-motion vector. The youngest ejecta from radio source I as traced by H₂O and SiO masers and CO has an axis orthogonal to its proper-motion vector. The axis of symmetry of the near-infrared bipolar reflection nebula associated with BN is also orthogonal to its proper motion.

Future numerical hydrodynamic or magnetohydrodynamic modeling of a massive star-forming clump in which a non-hierarchical system experiences a dynamical decay are needed to determine if the proposed scenario can indeed result in the generation of a powerful, explosive outflow such as is observed in Orion.

This work was supported by NSF grant AST0407356 and the CU Center for Astrobiology funded by NASA under Cooperative Agreement No. NNA04CC11A issued by the Office of Space Science. This paper is partially based on observations obtained with the Apache Point Observatory 3.5 m telescope, which is owned and operated by the Astrophysical Research Consortium. We thank Bruce Elmegreen, Hans Zinnecker, and Bo Reipurth for their insights and comments. We thank the referee, Paul Ho, for a thorough reading and very helpful criticism that greatly improved the manuscript.

Note added in proof. Since submission of this manuscript, a paper has appeared presenting new observations of the proper motions of radio sources BN, I, and n along with a numerical study of the dynamical decay in OMC1 and its implications for the origin of the outflow (Goddi et al. 2010).

Assume that a gaseous envelope surrounds a compact cluster of stars with mass M_cl, and that the motion of the gas is dominated by gravitational forces. The stars are all interior to the inner spherical boundary at radius r_in, and the average gas density is given by a power law so that ρ(r) = ρ₀r^−α out to an exterior radius r_out. If the total mass of the envelope is M_e, the normalization constant is

The potential energy of the envelope is found from

where M_tot = M_cl + M_e(r) is the total mass inside radius r before ejection of the cluster of massive stars. The first term is associated with the central cluster mass, while the second term is the potential associated with the gas envelope itself M_e(r) being the envelope mass inside radius (r). Taking the pre-stellar-ejection envelope to be virialized, we have −W_init = 2T_init. Because the observed stellar ejection velocities are much faster than the Kepler velocities in the envelope, stellar ejection happens on a timescale fast compared to the speed at which the gas can redistribute itself. Immediately after ejection of the stars, the potential energy of the envelope will be given solely by the second term, which remains the same, and the kinetic energy will still be given by T_init. If the post-ejection envelope were virialized, it would have −W_final = 2T_final,vir; the excess kinetic energy is then given by

The quantity |ΔW| is simply the term in Equation (A1) associated with the cluster mass, which integrates to give

To estimate the kinetic energy stored in the gravitationally bound motions of a circumstellar disks surrounding a star of mass M_⋆, assume that between the inner edge of the disk, r_d,in, and the outer edge at radius r_d,out, the average gas surface density is represented by a power law on the radius of the form Σ(r) = Σ₀r^−β. With a total disk mass M_d, the normalization is given by

Assuming a disk in near-Keplerian rotation (thus with a disk mass small compared to the stellar mass), the potential energy of the system is to a good approximation just due to the stellar potential,

With Keplerian motion we have −W_d,init = 2T_d,init. Immediately after ejection of the star, which we assume happens on a timescale fast compared to the gas-redistribution timescale, the potential energy of the disk will be due only to its self-potential, which we have neglected due to its assumed small contribution to the total. The excess kinetic energy is thus

Since the stellar potential was dominant, the quantity |ΔW_d| is simply the potential energy due to the star from Equation (B2), which integrates out to