OUTFLOW FEEDBACK REGULATED MASSIVE STAR FORMATION IN PARSEC-SCALE CLUSTER-FORMING CLUMPS

Our simulations are performed using the parallel AMR code, Enzo, with a newly added solver for unsplit conservative hydrodynamics and MHD (Wang & Abel 2009). As usual, at the heart of the MHD solver lies the enforcement of the divergence-free condition ∇ · B = 0. There are three classes of method for this purpose: constraint transport (CT; Evans & Hawley 1988), projection method, and hyperbolic clean (see, e.g., Tóth 2000 for a review). We have evaluated all these methods and found that divergence cleaning is most suitable for the problem at hand, which involves rather extreme variation in physical quantities (including the field strength) due to both gravitational collapse and protostellar outflows. The fast and robust method of hyperbolic clean of Dedner et al. (2002) is adopted. It enables us to follow the formation of hundreds of stars over several dynamical times. We monitor the value of ∇ · B to ensure that it is small in all of our MHD simulations.

Truelove et al. (1997) pointed out that in order to avoid artificial fragmentation, one need to use at least four cells to resolve the Jeans length. For a given cell size at the finest level of refinement, the condition translates into a maximum density, which we will call the "Jeans density," above which artificial fragmentation may happen. Since it is difficult in large-scale simulations such as ours to resolve the Jeans length all the way to the stellar density, sink particles are used to handle the gravitational collapse beyond the Jeans density. This approach has previously been successfully used in both SPH (e.g., Bate et al. 1995) and grid-based simulations of star formation (e.g., Krumholz et al. 2004; Offner et al. 2009).

Our procedure for implementing the sink particles is as follows. When a cell violates the Truelove criterion on the highest refinement level, a sink particle is inserted at the center of that cell. The initial mass of the sink particle is calculated such that after subtracting the sink mass, the cell will be at the Jeans density. The initial velocity of the sink particle is calculated using momentum conservation.

After creation, the sink particle accretes gas from its host cell according to a prescription inspired by the BH formula (Ruffert 1994),

$$\begin{equation} \dot{M}_{{\rm BH}} = 4\pi \rho _\infty r_{{\rm BH}}^2\sqrt{1.2544c_{\infty }^2 +v_\infty ^2}, \end{equation} \tag{ 1 }$$

where r_BH is the Bondi radius, ρ_∞, v_∞, and c_∞ are the gas density, velocity, and sound speed of the (uniform) medium far from the point mass.

Even though Equation (1) is not strictly applicable to the highly turbulent medium under investigation, we will follow Krumholz et al. (2004) and use it as a guide to construct a prescription for sink particle accretion. Since the cloud material is assumed isothermal, it is natural to set c_∞ to the isothermal sound speed c_s. For v_∞, we adopt the simple prescription v_∞ = v_cell − v_sink, where v_cell is the flow velocity in the cell and v_sink the velocity of the sink particle. We compute the Bondi radius in Equation (1) using r_BH = GM_sink/(c²_s + v²_∞), where M_sink is the sink mass. If r_BH is smaller than the cell length Δx, the cell density ρ_cell is used for ρ_∞. Otherwise, an extrapolation assuming an r^−1.5 density profile is used. In other words, we set ρ_∞ = ρ_cell min[1.0, (Δx/r_BH)^1.5], where Δx = 200 AU is the size of our finest cells. The amount of gas accreted during a single time step Δt is then $\dot{M}_{{\rm BH}}\Delta t$. The amount may or may not match the true value exactly at any given time, given the crudeness of the above prescription. This deficiency is corrected in a second step, through sink particle merging.

Sink particle merging is needed not only to ensure the correctness of the mass accretion algorithm, but also to save computation time when many sink particles are created. It is a tricky problem, because we want to eliminate, on the one hand, the artificial particles created by the mismatch between the true accretion rate and that given by Equation (1), and to preserve, on the other hand, the legitimate stellar seeds. It requires sub-grid information that is not available in the simulation. At this early stage in the development of the sink particle technique, our design principle of the merging recipe is to make it as simple as possible and with as few parameters as possible. We have experimented with a number of recipes, and settled on one with two parameters: a merging mass M_merg and a merging distance l_merg. We divide the sink particles into two groups: the "small" ("big") particles have masses smaller (larger) than M_merg. At every time step, the merging is done in two steps. First, for each small particle, we search within the merging distance for its nearest big particle. If the search is successful, we merge this particle with the located big particle. Second, for all small particles that remain after the first step, we group them using a Friend-of-Friend algorithm with the merging distance l_merg as the linking length (Davis et al. 1985). In this work, we use M_merg = 0.01 M_☉ and l_merg = 5Δx, which is about 1000 AU. Experimentation shows that increasing M_merg to 0.1 M_☉ or l_merg to 10Δx does not change the star formation rate or the stellar mass distribution significantly. Lowering the merging distance to 3Δx has little effect on the total star formation rate, but can change the mass of a star by up to 50%, probably because the flow pattern within a few cells of a sink particle is not well resolved. The flow so close to a forming star may be strongly affected by radiative heating (Offner et al. 2009; Bate 2009; Smith et al. 2009), which is not included in our simulations. For these reasons, we believe that the properties of circumstellar disks and binaries may not be reliably captured by our simulation; they are not discussed further in the paper.

With the inclusion of sink particle merging, the crude prescription for sink particle accretion based on the BH formula (Equation (1)) does not appear to pose any serious problem. We have experimented with changing the parameters in the formula and did not see any significant difference in either the stellar mass accretion rate or the final stellar mass distribution. As stressed by Krumholz et al. (2004), the reason is probably that, when the BH formula underestimates the true accretion rate, the sink particle creation routine would transform the unaccreted gas into small sink particles, which would merge quickly with the original particle and restore the correct rate. On the other hand, when the BH formula overestimates the true rate, which happens rarely and only when the sink particle is massive, the accretion flow is typically supersonic near the particle, and the over-accretion does not affect the flow further away. In particular, we have repeated the test of the collapse of an singular isothermal sphere (SIS) described in section Section 3.1 of Krumholz et al. (2004), and found that the numerically obtained accretion rate agrees with the analytic value to within a few percent. To avoid numerical instabilities, we restrict the maximum amount of the gas accreted in a single time step to be less than 25% of the cell mass.

Finally, the gravity of the sink particles is calculated using the standard particle-mesh method, and their positions and velocities are updated using the leapfrog method.

Protostellar outflows are a crucial element of our model. In large-scale global simulations such as ours, it is unfortunately impossible to follow their production from first principles. An observationally motivated prescription is used instead. We assume that the protostellar outflow momentum injection rate is proportional to the stellar mass accretion rate, with a mass dependent proportionality constant P_* = P₀(M_*/M₀)^1/2, where P₀ = 16 km s⁻¹ and M₀ = 1M_☉. The mass dependence is chosen to reflect the fact that the outflows from high mass stars tend to be more powerful even after correcting for their larger masses (Richer et al. 2000, D. Shepherd 2007, private communication). We represent the outflow by injecting a momentum ΔP = P_*ΔM into the surrounding gas every time the sink particle has gained a mass increment ΔM. Ideally, one would like to make ΔM as small as possible, to render the injection continuous. However, we find continuous injection is too time consuming, because of the small Courant time step associated with the fast outflow speed and small grid size near the sink particles. Furthermore, if ΔM is much less than the mass of an injecting cell, the injected momentum would have little dynamical effect on the cell. These considerations led us to adopt ΔM = 0.1 M_☉. We have experimented with ΔM = 0.2 M_☉, and it did not change the results significantly. To avoid too high an outflow speed when a surrounding cell happens to have a low density, we take 10% ΔM out of the sink particle, and divide it evenly between the injecting cells. To be conservative, we represent the outflow from each accreting protostar by a bipolar jet, with two oppositely directing jet beams. The outflow in each beam is injected in the same direction; a wider injection angle should lead to a more efficient turbulence driving. For simplicity, the jet axis is chosen to be the local magnetic field direction of the host cell of the sink particle when the jet is first injected (and is kept fixed at all times), and its width is set to five finest cells.

We consider a moderately condensed clump of molecular cloud with an initial density profile of

$$\begin{equation} \rho (r) = {\rho _c\over 1+(r/r_c)^2}, \end{equation} \tag{ 2 }$$

where r_c = L/6 is the radius of the central plateau region and L = 2 pc is the length of the simulation box. We adopt a central density ρ_c = 1.0 × 10⁻¹⁹ g cm⁻³, which corresponds to a central free-fall time $t_{\rm ff,c}= \sqrt{3\pi /(32G\rho _c)}=0.21$ Myr. It yields a clump mass within a sphere of 1 pc in radius of M_cp = 1215 M_☉. The average density within that sphere is ${\bar{\rho }}= 1.96\times 10^{-20}$ g cm⁻³, corresponding to an average free-fall time of t_ff,cp = 0.48 Myr, which we will identify with the clump free-fall time. The ratio of the mass and free-fall time yields a characteristic mass accretion rate ${\dot{M}}_{\rm ff, cp}=2.55 \times 10^{-3}$ M_☉ yr⁻¹ for our model clump. Since a square box is used, the total mass in the simulation domain is somewhat larger than that within the sphere, with M_tot = 1641 M_☉. Our choice of cloud density, mass and mass distribution is motivated by observations of infrared dark clouds (e.g., Butler & Tan 2009), which are thought to be future sites of cluster and massive star formation; they tend to be centrally condensed. Centrally condensed density profiles are also inferred for dense clumps before the onset of massive star formation from dust continuum emission (e.g., Beuther et al. 2002).

The equation of state is isothermal, with a sound speed c_s = 0.265 km s⁻¹, which corresponds to a temperature of T = 20 K for a mean molecular weight μ = 2.33.

With a density profile given by Equation (2), our model clump has a gravitational potential energy E_G = 0.7GM²_cp/R, which yields a virial parameter α_vir = 2E_K/E_G = Rσ²/(0.7GM_cp), where σ is the three-dimensional rms velocity. In terms of the three-dimensional Mach number $\mathcal {M}\equiv \sigma /c_s$, we have

$$\begin{equation} \alpha _{\rm vir} = 1.5 \left({\mathcal {M}\over 9}\right). \end{equation} \tag{ 3 }$$

The initial Mach number is set to $\mathcal {M}=9$. Following the standard treatment in supersonic turbulence simulations (e.g., Mac Low et al. 1998), we impose a velocity field of power spectrum v²_k ∝ k⁻⁴ with a minimum wavenumber 2 and maximum wavenumber 10 in units of inverse box length.

In models with a magnetic field, the field is assumed to be initially uniform in the z-direction, with a strength given by

$$\begin{equation} B_0 = 1.0\times 10^{-4} \left({\lambda _{\rm tot}\over 1.4}\right)^{-1} \ \ {\rm gauss}, \end{equation} \tag{ 4 }$$

where λ_tot ≡ 2πG^1/2M_tot/(πB₀L²/4) = 1.4 is the mass-to-flux ratio in units of the critical value. Because the adopted density profile is centrally condensed, the central part of the clump within r_c has a larger mass-to-flux ratio λ_c = 2.5λ_tot = 3.3. Thus, the envelope is more magnetized than the central part. Our adopted degree of magnetization is in the range inferred by Falgarone et al. (2008).

We note that our simulations can be scaled to represent clumps with different sets of dimensional quantities. For example, if we keep the clump mass the same but increase the temperature by a factor of 2, from 20 K to 40 K, the clump dimension would shrink by a factor of 2, from 2 pc to 1 pc. The density would increase correspondingly, by a factor of 8, which would lead to a decrease in free-fall time by a factor of $2\sqrt{2}$ and an increase in mass accretion rate by the same factor.

For all simulations, the top grid has resolution 128³ and the maximum refinement level is 4, which corresponds to a highest resolution of 200 AU. With this resolution and our sink particle algorithms, the risk of over producing low-mass stars is reduced in our model even though we do not include radiative heating, which happens predominantly on the disk scales (Krumholz et al. 2007; Bate 2009). We will not be confident about the mass spectrum at the lowest mass end. However, the formation of massive stars, the main focus of this paper, is adequately resolved.

Ideally one would want to start the simulation from diffuse atomic gas and follow the build up of the turbulent molecular cloud including the relevant physics. This would eliminate potential concerns about boundary effects. However, for our numerical experiments that focus on small parsec-scale clumps, specific choices about the boundary conditions have to be made. After experimentation, we settled on periodic boundary conditions for the hydro quantities to conserve the mass in the box over the simulation time. Periodic boundary conditions are also chosen for the gravity to reflect the expectation that clumps are probably not formed in isolation. They prevent the clump from developing large infall speed (and associated low density and high Alfvén speed) near the edge. In the model with protostellar outflows, to prevent the high-speed outflows from re-entering the other side of the boundary and thus artificially increasing their effects, we reduce their speeds by a factor of 10 when they reach the boundary.

We consider four models of increasing complexity, starting with a base model (BASE) that is initially neither turbulent nor magnetized. We then add to this base model, one by one, an initial turbulence (HD), magnetic fields (MHD), and outflows (WIND). The models are summarized in Table 1. In the next two sections (Sections 3 and 4), we will focus on presenting and analyzing the numerical results. We postpone an in-depth discussion of the connection between our work and the two widely discussed scenarios of massive star formation—turbulent core (McKee & Tan 2003) and competitive accretion (Bonnell et al. 2004)—to Section 5.

Massive stars are formed in all of the four simulations listed in Table 1. In the limit of no turbulence (model BASE), all of the accreted mass goes to a single object at the center. At the end of the simulation (t = 2t_ff,c = 0.42 Myr), the object has grown to more than 400 M_☉, which is clearly unphysical. The addition of an initial turbulence in model HD enabled the formation of eight massive sink particles by t = 3t_ff,c = 0.63 Myr. The time evolution of the positions of these particles are shown in Figure 3. There are two features worth noting. First, only five lines are clearly distinguishable in the plot. This is because the first four of the eight stars are formed close together both in space and in time; they stay close together at all times, so that their trajectories on the plot are indistinguishable. For our main purpose of evaluating mass accretion rate, we will treat them as a single object. Second, the remaining four massive stars are formed at different times and locations. Two of them join the first object at later times to form a tight group, while the other two merge into another system. We refer to the merged systems as group A and B, respectively, and the massive stars formed at different times in each group as different generations. For example, A1 is the first generation of massive stars to form in group A. We will refer to all massive stars in a single generation of a group as "a massive object" (because they form at essentially the same time and same location, and the partition of mass between the stars depends somewhat on the sink particle treatment). Some properties of the massive objects are listed in Tables 2–4 at three different times t = 0.63, 0.84, and 1.16 Myr.

Zoom In Zoom Out Reset image size

There are drastic differences between the massive objects formed in different simulations. At a relatively early time t = 0.63 Myr (see Table 2), five massive objects (eight massive stars) form in the HD model, with a total mass of 131 M_☉, which is 37% of the stellar mass. The presence of a magnetic field in model MHD reduces the number of massive objects (stars) to three (four), and their total mass to 73.6 M_☉, although the mass fraction of the massive stars remains similar (45%). When outflows are turned on (model WIND), the formation of massive stars is suppressed up to this time, with no stars more massive than 10 M_☉.

Massive stars do form at later times in model WIND. At t = 0.84 Myr, there are two massive stars (16.9 and 10.5 M_☉ each) formed, with a total mass of 27.4 M_☉, which is much less than the total mass of massive stars formed at the same time in model MHD (183 M_☉). The difference re-enforces the notion that the outflows retard massive star formation (see also Dale & Bonnell 2008). The massive stars do continue to grow with time, however. Their masses reach 46.4 and 16.7 M_☉, respectively, near the end of the simulation at t = 1.16 Myr (about 2.4 global free-fall times). They are each joined by another massive star, forming two loose groups. Together, the four massive stars contain 86.9 M_☉, accounting for 35% of all the stellar mass in the cluster.

How do massive stars form in our simulations? To address this question, we will concentrate on the most massive object formed in each simulation, starting with the simpler, HD model that has turbulence but neither magnetic fields nor outflow feedback.

4.2.1. Model HD

The most massive object, HD-A1 in Table 2, formally reached a mass of 61.1 M_☉ at the end of the simulation (t = 3t_ff,c = 0.63 Myr). Its seed is among the first to form (at a time t ∼ 0.2 Myr), as found by Bonnell et al. (2004) in their hydro simulation. The average mass accretion rate is therefore ${\dot{M}}_{\rm ave}\approx 1.5\times 10^{-4}$ M_☉ yr⁻¹. The instantaneous accretion rate for this object is shown in Figure 4. It quickly rises to a plateau of ∼10⁻⁴ M_☉ yr⁻¹, before jumping at t ∼ 0.5 Myr to a second plateau that is ∼2.5 higher. About equal amounts of mass (∼30 M_☉) are accreted in the two plateau phases. During the first plateau phase, the stellar mass increased by an order of magnitude, while the mass accretion rate stayed more or less constant. This behavior is characteristic of core collapse (such as the collapse of the SIS where the mass accretion rate is a constant; Shu 1977) rather than the BH type accretion which depends sensitively on the stellar mass. The rapid increase in mass accretion rate around ∼0.5 Myr is not triggered by any sudden increase in stellar mass, which again points to an external control of the mass accretion rate; it is caused by the onset of global collapse of the clump.

Zoom In Zoom Out Reset image size

To examine the nature of the accreting flow before the global collapse in more detail, we plot in Figure 5 the mass distribution, velocity field, and gravitation potential at t = 1.5t_ff,c = 0.32 Myr, when the eventual massive object A1 is still in an early stage of accretion. It is embedded in a dense filament, which is produced by converging flows in the turbulence (see also Banerjee et al. 2006). Here, the inclusion of self-gravity in the calculation is crucial. It allows the converged material to stay together, rather than re-bouncing and dispersing quickly; in other words, the self-gravity softens the high-speed collision that creates the filament.

Zoom In Zoom Out Reset image size

The filament feeds the object at a high rate of ∼10⁻⁴M_☉ yr⁻¹, some 20 times the standard Shu's rate for an SIS. As stressed by Banerjee et al. (2006), the high accretion rate is due to the rapid formation of non-equilibrium structures (dense filaments) that are fed continuously by converging flows (and are thus very different from the equilibrium SIS). Their existence does not depend on the stellar objects embedded in them; rather, the objects are produced by the runaway collapse of the densest parts of the filaments, often at the intersections with other filaments. In this sense, the formation of the most massive object in this early phase is not due to competitive accretion, since the gravity of the object is not the primary driver of the mass accumulation in the filaments; there are few objects to compete against at this stage of star formation in any case. It may be viewed as a special version of core collapse, in which the "core" is an over-dense, highly out-of-equilibrium filament that existed before the object and that continues to grow out of the converging flow even as part of it disappears into the collapsed object. In this picture, the massive star formation is simply part of the rapid "core" formation process. The use of sink particles in our simulation enables us to go beyond the calculations of Banerjee et al. (which stopped before significant mass accretion onto the stellar object) and follow the mass accretion for a long time, and to reveal a second phase of even more rapid accretion, driven by global collapse.

The global collapse at the end of the pure hydro simulation (t = 0.63 Myr) was already shown in Figure 2. It is further illustrated in Figure 5, which shows the distribution of mass, velocity field, and contours of gravitation potential in a sub-region centered on the most massive object, with the positions of stars superposed. It is clear that the global collapse has supplied the central region near the bottom of the gravitational potential well with plenty of dense material, which fuels both the enormous rate of total stellar mass accretion (close to the characteristic free-fall rate, see Figure 1) and the high accretion rate of the most massive object, which is but one of many stars in the region. In this crowded region, competitive rather than core accretion is likely at work. The reason is that the stars in this region are not enveloped by permanent host cores, because the dense gas in the region is constantly drained onto stars and constantly replenished by the global collapse. The most massive object is located near the minimum of the global potential well and accretes the global collapse-fed material at a rate higher than any other object. Nevertheless, its accretion rate is only ∼10% of the total rate. In other words, the vast majority of the globally collapsing flow is diverted to stars other than the most massive object. Its accretion rate is large (∼2.5 × 10⁻⁴ M_☉ yr⁻¹) during this late phase only because of an even larger global collapse rate.

A common theme of the early and late phases of rapid accretion is that the dense gas that feeds the most massive object at a high rate is gathered by an agent external to the object. It is the converging flow set up by the initial turbulence in the early phase and the globally collapsing flow in the late phase. In this aspect, the formation mechanism is closer to core collapse than to competitive accretion. One may plausibly identify the dense filament in the early phase and the dense region at the bottom of the global gravitational potential well in the late phase as a McKee–Tan core (McKee & Tan 2003). However, the "cores" so identified are transient objects that are not in equilibrium. They are evolving constantly, with mass growing (from converging or collapsing flow) and depleting (into one or more collapsed objects) at the same time. The replenishment of dense gas in the "core" is the key to the long duration of high mass accretion rate for the most massive object, which is many times the free-fall time of the small dense "core" and comparable to the global collapse time of the clump as a whole. In other words, there is no hint of any drop in the mass accretion rate after a finite reservoir of core material gets depleted over a (short) core free-fall time. Indeed, the nearly constant or increasing accretion rate for the most massive object highlights the important issue of when and how the mass accretion is terminated; there is an urgent need to tackle this issue if the final stellar mass is to be determined. We will return to a more detailed discussion of the nature of the massive star's high accretion rate in the absence of magnetic fields and outflow feedback in Section 5.1.

4.2.2. Model MHD

The effect of the magnetic field in model MHD on the mass accretion rate onto the most massive object is relatively modest (see Figure 4). As in the HD case, there are two plateau phases. The first starts somewhat later than that in the HD case, indicating that the initiation of star formation in general and massive star formation in particular is somewhat delayed, as expected, because of magnetic cushion of turbulent compression. Once started, the accretion rate onto the most massive object is, in fact, slightly higher in the MHD case than in the HD case, indicating that the magnetic field does not significantly impede the initial phase of rapid mass accretion, which is enabled by the dense structures formed by turbulent compression. Indeed, dense structure formation is aided by a strong magnetic field, which forces the turbulent flows to collide along the field lines (see, e.g., also Tilley & Pudritz 2007; Price & Bate 2008). When viewed in three-dimensional, the dense structure resembles a warped, fragmented disk, with transient spirals that come and go. The spirals indicate significant levels of rotation on relatively small scales, which may hinder the stellar mass accretion due to centrifugal barriers.

The centrifugal barrier may be weakened or even removed by magnetic braking. Evidence for the braking is shown in Figure 6, which shows a magnetic bubble driven from the central region, where active mass accretion is occurring. Magnetic braking driven bubbles provide a form of energy feedback that should accompany any star formed in a magnetized cloud (e.g., Draine 1980; Tomisaka 1998; Banerjee & Pudritz 2006; Mellon & Li 2008; Hennebelle & Fromang 2008) but is absent from purely hydrodynamic models. We expect this form of feedback to become more important with a higher spatial resolution, since the rotating material closer to a star would be able to wrap the field lines more quickly (until the flux freezing under the ideal MHD assumption breaks down). Indeed, the protostellar outflow feedback may be viewed as a special form of this magnetic bubble-driven feedback, since the outflow is thought to be driven by the field lines that are forced to rotate rapidly by the circumstellar disks (Konigl & Pudritz 2000; Shu et al. 2000). Besides driving the magnetic bubble, the field can significantly modify the formation, evolution, and fragmentation of the disks around accreting massive stars compared to the pure hydro case, although higher resolution simulations focusing on the formation of individual massive stars are needed to quantify the effects.

Zoom In Zoom Out Reset image size

Despite the feedback through the magnetic bubbles, large-scale collapse does occur at late times. When enough mass has been accumulated near the bottom of the gravitational potential well, the self-gravity of the accumulated gas overwhelms the magnetic tension forces, leading to a rapid cross-field collapse and an elevated, second phase of rapid mass accretion. The mass accretion is facilitated by the magnetic braking, which enables the high angular momentum material to sink closer to the center than it would be in the absence of the braking. The collapse is not as global as in the HD case, however, as shown in Figure 7. At the end of the MHD simulation (t = 0.84 Myr), the bulk of the clump material infalls toward the bottom of the potential well at a speed of order twice the sound speed or less, unlike in the HD case, where the infall is faster and more widespread. The collapse remains significantly retarded by the magnetic tension except deep in the potential well.

Zoom In Zoom Out Reset image size

4.2.3. Model WIND

The combined effect of outflows and magnetic fields on the formation of the most massive object is much greater than that of magnetic fields alone. From Figure 4, we find that, even though a massive star of more than 45 M_☉ was eventually formed at the end of the WIND simulation, it took nearly 10⁶ yr to accrete the mass, with an average rate of only ∼5 × 10⁻⁵ M_☉ yr⁻¹. The accretion rate at early times is even lower, with a value of ∼2 × 10⁻⁵ M_☉ yr⁻¹ during the initial third of the time; it increases to ∼4 × 10⁻⁵ M_☉ yr⁻¹ during the second third. These rates are about a factor of 5 lower than those in the HD and MHD models at comparable times. The question is: why does the most massive object accrete at such a low rate?

The basic reason is of course the outflows in the WIND model, which are much more powerful than the magnetic braking driven bubbles that exist in the MHD case. The outflows modify both the mass distribution and velocity field, and thus the outcome of the gravitational dynamics. As a result, the most massive object forms at a somewhat different time and location. Despite the additional physics in both MHD and WIND models, it is initially still embedded in a dense filament, as it is in the HD case. As discussed earlier, the filament is the key to the first, turbulence compression-induced phase of rapid mass accretion in the HD and MHD models. In the WIND model, accretion from the filament is reduced in two ways. First, there are several stars formed in the filament. They all emit outflows which tend to break the filament into smaller segments. Second, the outflows emerge more easily perpendicular to (rather than along) the filament. Once they break out, they blow against the converging flow and slow down the mass accumulation that feeds the growth of the filament in the first place.

As the outflows propagate to large distances, a fraction of their momentum is lost through the computation boundary. The remaining fraction is deposited in the lower density region that surrounds the denser region of active star formation near the bottom of the potential well. The momentum deposition allows the bulk of the clump material to be supported against rapid global collapse. The absence of a global collapse is illustrated in Figure 8, which shows a rather chaotic velocity field that involves both infall and outflow. It is further quantified in Figure 7, which shows that the mass-weighted, azimuthally averaged, infall speed toward the most massive object is subsonic over most of the volume,⁵ except within a radius of about 0.2 pc, where the speed becomes two to three times the sound speed. The slower average infall leads to a smaller amount of dense gas accumulating near the bottom of the potential well. The accumulated dense gas is also more fragmented because of interaction with outflows. Both factors limit the mass accretion rate onto the most massive object, especially at late times.

Zoom In Zoom Out Reset image size

As in the HD and MHD cases, the most massive object is not formed out of a pre-existing dense core. The case against the pre-existing dense core picture is stronger in the WIND case, because it takes longer (nearly two global free-fall times) to accumulate the final mass, whereas a pre-existing dense core should collapse and exhaust its mass in a local (core) free-fall time, which should be a very small fraction of the global free-fall time. The slow accretion rate and long accretion time are clear evidence that the formation of the most massive object is controlled by the global clump dynamics in our WIND simulation.

Massive stars form quickly in our simulated parsec-scale clump in the absence of any magnetic field or outflow feedback (see Model HD). They are fueled by high mass accretion rates from either the dense filaments that are formed by turbulent compression or the clump-wide collapse due to the global turbulence dissipation. In both cases, the dense material that is depleted onto the massive star is constantly replenished. The replenishment is illustrated in the left panel of Figure 9, where the mass of the most massive object is plotted as a function of time, along with the mass of the gas within a "core" of 0.1 pc in diameter around the object; the size was chosen to coincide with the fiducial value that McKee & Tan (2003) used to define their "turbulent cores." It is clear that the "core" so defined has an initial mass that is well below the final mass of the most massive object, and thus cannot supply all of the mass of the object. As the object gains mass, the mass of the "core" stays roughly constant or even increases (rather than decreases), which demands that the core mass be replenished or fed from the surrounding clump. We are thus motivated to term this model the "clump-fed massive star formation" model.

Zoom In Zoom Out Reset image size

Our "clump-fed massive star formation" model (CF model for short) contains elements of the two widely discussed scenarios for massive star formation in the literature: the turbulent core model of McKee & Tan (2003) and competitive accretion model of Bonnell et al. (2004). It has in common with the turbulent core model that the mass accretion rate onto the massive star is not primarily determined by the star itself, but rather by the properties of the pre-existing gas that produces the seed of the star in the first place and that continues to feed the stellar growth; in other words, if the massive star were to be removed prematurely, another seed would be produced in its place and grow to a high mass. Indeed, this phenomenon was observed in the simulation of Smith et al. (2009), where the initial sink particle which was preferentially growing due to infall was displaced by a dynamical interaction before it became too massive; in its place, a new pre-stellar core started to accrete the mass and became the most massive sink (R. Smith 2009, private communication) It is primarily the properties of the fueling gas that determine the stellar mass accretion rate and thus the stellar mass. McKee & Tan envisioned this gas to be a pre-existing massive turbulent core. It collapses to form the massive star, with a mass accretion rate that depends on the core structure. For example, if the turbulent core has a density profile of r^−1.5, the mass accretion rate would increase linearly with time t (McKee & Tan 2003). In our CF model, the pre-existing gas is the cluster-forming clump, which produces the transient dense material that feeds the massive star at a high rate through two mechanisms: (1) the collapse of the dense filaments produced by turbulent compression, as already emphasized by Banerjee et al. (2006) and (2) global collapse, driven by the turbulence decay in our simulations, but can in principle be caused by an external compression as well. The first mechanism depends on the detailed properties of the initial turbulence in the clump, which are not well constrained observationally. The second mechanism depends on the dissipation of (supersonic) turbulence that is inevitable, and should be more robust.

Our unregulated CF model has in common with the competitive accretion model that a mass star can form even in the absence of a pre-existing turbulence-supported massive core and that the global gravitational potential and dynamics of the clump play an important, even the dominant, role in massive star formation. The latter is especially true at late times, when the global collapse feeds mass at a high rate to the compact region near the bottom of the potential well, where a large number of stars (including massive stars) are already present. In the absence of stellar feedback, competition for this clump-fed material in the crowded region is unavoidable. The most massive object tends to grow the fastest, mainly because it typically locates closest to the center of the potential, as emphasized by Bonnell et al. (2007). Although the numerical results of our grid-based AMR hydro simulations are in broad agreement with those of SPH simulations of Bonnell et al. (2004), we differ from their interpretation of the results regarding the formed stars, particularly as it relates to massive star formation. When a massive star is formed, the vast majority of the clump mass remains in the gas. The high accretion rates of the massive stars at late times of our hydro simulation derive directly from the global gas collapse. Even if there were no stars near the center of the collapse (or more likely the mass accretions onto individual stars are terminated by the stellar feedback), (other) massive stars can still form from scratch, as long as the collapse delivers mass to the center at a high enough rate. So we argue that it is predominantly the structure and dynamics of the gaseous component that set the relevant physics in forming the massive stars rather than the properties of the stars made previously. A similar point was made independently by Smith et al. (2009).

A moderately strong magnetic field (corresponding to the observationally inferred dimensionless mass-to-flux ratio of a few) does not qualitatively change the clump-fed picture for massive star formation by itself. As discussed in Section 4.2.2, the magnetic field has relatively little effect on the high mass accretion rate onto the most massive object induced by the initial turbulent compression and filament formation. Indeed, it is conducive to filament formation by guiding the converging flows to collide along the field lines. The magnetic field does slow down the global collapse-induced rapid accretion by a modest factor of a few. The basic ingredients of the CF model remain, however, as can be seen from the middle panel of Figure 9. Specifically, the initial core mass is still much less the final mass of the most massive object, and most of the mass that goes into the object still needs to be replenished or fed from the surrounding clump. We conclude that the feeding processes for massive star formation are only weakly regulated by the magnetic field, perhaps because the field is only moderately strong (i.e., the clump is moderately supercritical), and it resists the feeding passively rather than actively (unlike the outflows, see below). The magnetic field does change the dynamical coupling between different parts of the clump, an important aspect of cluster and massive star formation that we plan to return to in a future investigation.

Massive star formation through rapid accretion of the mass that is fed from outside the small region surrounding the forming star is particularly vulnerable to outflow feedback. This is because the feeding is directly opposed by the feedback. In the case of the rapid accretion fed by turbulent compression and filament formation, the filament is quickly chopped up into small segments by the outflows driven by the stars formed in it. The outflows also slow down any further mass accumulation in the filament after star formation is initiated. As a result, this mode of turbulent compression-fed massive star formation is strongly regulated, perhaps even choked off completely, by outflow feedback.

In the case of the rapid accretion fed by global collapse, the infall is countered by the outflows on all scales, especially on the global clump scale. The additional global support provided by the outflow feedback can reduce the total star formation rate by a large factor. This reduction affects the formation of all cluster members, especially the massive stars. This is because the massive stars receive a larger fraction of the collapse-fed material, and are thus more sensitive to the change of global clump dynamics. They also tend to complete their formation toward the end of cluster formation (even though their seeds tend to be among the first objects to form), making them more prone to the accumulative influence of multiple generations of outflows that precede their eventual formation.⁶ Since the outflows are believed to be driven by the release of gravitational binding energy from mass accretion in one form or another (Konigl & Pudritz 2000; Shu et al. 2000), and most of the accreted mass goes to low-mass (rather than high mass) stars for a Salpeter-like initial mass function (IMF), a large fraction if not the bulk of the outflow feedback should come from the low-mass stars. In this sense, the formation of low-mass stars in a dense clump can profoundly influence the formation of massive stars in the same clump, through their feedback on the clump dynamics. Whether massive stars eventually form or not in a dense clump depends on the extent to which all the outflows in the clump collectively regulate the global collapse and slow down the star formation. If the total star formation rate is reduced, for example, below the fiducial minimum rate for massive star formation, 10⁻⁴ M_☉ yr⁻¹, no massive stars would form at all. In the opposite extreme where the outflows are weak and the feedback is not strong enough to reverse the global collapse, stars (especially massive stars) would form quickly, as in our pure hydro simulation.

The degree of outflow regulation will depend on the properties of the outflows, including their strengths and degrees of collimation, both of which are somewhat uncertain. What we have demonstrated explicitly through numerical simulation is that for well-collimated jets of reasonable strength the outflow feedback can prevent rapid global collapse, and keep the total star formation rate an order of magnitude below the characteristic free-fall rate. Massive stars do eventually form in our simulation that includes outflow feedback. It demonstrates that outflows can strongly regulate massive star formation, but do not necessarily quench it completely, especially in dense massive clumps with a high characteristic free-fall rate ${\dot{M}}_{\rm ff, cp}$. For such clumps, even when the bulk of the clump material is supported, a small fraction can still percolate down the global gravitational potential well, feeding the formation of massive stars near the center at a high enough rate.

The fact that our outflow-regulated massive star formation remains clump-fed rather than core-fed is illustrated in the right panel of Figure 9. As in the HD and MHD models that do not have outflow feedback, the mass of the 0.1 pc-sized "core" that surrounds the most massive object in the WIND model is initially smaller than the final mass of the object, and does not decrease monotonically as the object accretes. Again, the bulk of the accreted stellar material must come from outside the "core," which is an essential feature of our new scenario of ORCF massive star formation. In this picture, the parsec-scale cluster-forming clump, rather than a 0.1 pc-sized turbulent core, is the basic unit for massive star formation. Indeed, the 0.1 pc-sized "core" does not appear to be a dynamically distinct region in our simulation (see Figure 8); it simply marks the inner part of the clump-wide potential well that extends from near the most massive object to large, parsec-scale distances.

As pointed out by Bonnell et al. (2007), a potential drawback of the turbulent core model is that massive, turbulence-supported, cores tend to fragment into many stars rather than collapse monolithically, although radiative feedback from the rapidly accreting massive stars can reduce the level of fragmentation (Krumholz & McKee 2008). Fragmentation is expected to be less of a problem in our clump-fed picture of massive star formation, since the material near the forming massive stars does not have to be supported by a strong turbulence; it is typically in a state of rapid collapse that feeds the growing massive stars at a high rate, even though the clump as a whole may remain supported by (possibly outflow-driven) turbulence and, perhaps to a lesser extent, magnetic fields. Fragmentation on the clump scale does occur. It produces a cluster of lower mass stars that accompany the massive stars. Our ORCF model can thus be viewed as a "turbulent clump model" for massive star formation, with the clump turbulence regulated by the outflow feedback from the cluster members. If the cluster-forming parsec-scale dense clumps are indeed the basic units of massive star formation, how they form in the first place becomes an important problem that deserves close theoretical and observational attention.

The most severe limitation of the current work is perhaps the neglect of radiative feedback. Radiative heating changes the clump fragmentation behavior, especially close to the forming stars (Krumholz et al. 2007; Bate 2009). This effect was mimicked to some extent by the relatively large sink particle merging distance adopted in our simulations (with a length of five cells or 10³ AU), which suppresses fragmentation on the small (mostly disk) scale. Furthermore, if our ORCF scenario is correct, the formation of massive stars may be more sensitive to the global clump dynamics (which are less affected by the radiative heating) than the gas properties close to the stars. Nevertheless, we believe that the formation of massive stars will benefit from the suppression of fragmentation by radiative heating, especially near the bottom of the gravitational potential well of the clump, where the thermal Jeans mass is formally smaller than the typical stellar mass. On small scales, radiative pressure may slow down the mass accretion onto individual massive stars somewhat, although the accretion of magnetized material may be enhanced to some extent by non-ideal MHD effects, such as ambipolar diffusion, which are not included in our ideal MHD calculations.

Another effect that we neglected was the H ii region driven by the massive star's UV radiation. The expansion of H ii regions provides a way to remove the clump gas, and perhaps terminate the cluster formation. It needs to be included in future simulations that aim to model the entire history of cluster formation (see, e.g., Gritschneder et al. 2009). Such studies may also need to include massive star winds, which are observed to have dramatic effects in some regions (e.g., the Carina Nebula, see Smith & Brooks 2008 for a review). They are the main alternative to the H ii regions as the means for terminating the cluster formation.

A further limitation is the periodic boundary condition used in our simulation. It precludes any communication between the cluster-forming clump and its surrounding environment. If there is energy injection into the dense clump from the ambient medium, the externally supplied energy may aid the outflows in regulating the cluster and massive star formation. Large-scale external compression may, on the other hand, lead to rapid clump collapse and massive star formation. In this case, the feeding of massive stars may extend beyond the parsec-scale clump, and the massive star formation may simply be part of the rapid clump formation.