LIMITING ACCRETION ONTO MASSIVE STARS BY FRAGMENTATION-INDUCED STARVATION

Understanding the formation of massive stars is of pivotal importance in modern astrophysics. Massive stars are rare and short lived. However, they are also very bright and allow us to reach out to the far ends of the universe. For example, the most distant galaxies in the Hubble Ultra Deep Field (Beckwith et al. 2006) are all characterized by vigorous high-mass star formation. Understanding the origin of massive stars, at present and at early times, therefore is a prerequisite to understanding cosmic history. In our own Milky Way, high-mass stars contribute the bulk of the UV radiation field. This radiation can ionize hydrogen and dissolve molecules. Expanding H ii regions, bubbles of ionized gas surrounding massive stars, have been identified as an important source of interstellar turbulence. The complex interplay between gas dynamics and radiation thus is a key element of the evolution of the interstellar medium (ISM) and the Galactic matter cycle. Furthermore, stars are the primary source of chemical elements heavier than the hydrogen, helium, and lithium that were produced in the big bang. Again, massive stars have a disproportionally large share in the enrichment history of the Galaxy and the universe as a whole. It is the violent supernova explosions associated with the death of high-mass stars that contribute most of the metals. Because they insert large amounts of energy and momentum input, supernovae at the same time strongly stir up the ISM and provide for the effective mixing of the newly bread elements.

Despite their importance, the physical processes that initiate and control the build-up of massive stars are still not well understood and subject of intense debate (Mac Low & Klessen 2004; Zinnecker & Yorke 2007; McKee & Ostriker 2007). Because their formation time is short, of order of 10⁵ yr, and because they form deeply embedded in massive cloud cores, very little is known about the initial and environmental conditions of high-mass stellar birth. In general, high-mass star-forming regions are characterized by more extreme physical conditions than their low-mass counterparts, containing cores of size, mass, and velocity dispersion roughly an order of magnitude larger than those of cores in low-mass regions (e.g., Jijina et al. 1999; Garay & Lizano 1999; Kurtz et al. 2000; Beuther et al. 2007; Motte et al. 2008). Typical sizes of cluster-forming clumps are ∼1 pc, they have mean densities of n ∼ 10⁵ cm⁻³, masses of ∼10³ M_☉ and above, and velocity dispersions ranging between 1.5 and 4 km s⁻¹. Whenever observed with high resolution, these clumps break up in even denser cores, that are believed to be the immediate precursors of single or gravitationally bound multiple massive protostars.

Massive stars usually form as members of multiple stellar systems (Ho & Haschick 1981; Lada 2006; Zinnecker & Yorke 2007) which themselves are parts of larger clusters (Lada & Lada 2003; de Wit et al. 2004; Testi et al. 1997). This fact adds additional challenges to the interpretation of observational data from high-mass star-forming regions as it is difficult to disentangle mutual dynamical interactions from the influence of individual stars (e.g., Goto et al. 2006; Linz et al. 2005). Furthermore, high-mass stars reach the main sequence while still accreting. Their Kelvin–Helmholtz pre-main-sequence contraction time is considerably shorter than their accretion time. Once a star has reached a mass of about 10 M_☉ its spectrum becomes UV dominated and it begins to ionize its environment. This means that accretion as well as ionizing and non-ionizing radiation needs to be considered in concert (Keto 2002, 2003, 2007; Peters et al. 2010a, 2010b). It has been realized decades ago that in simple one-dimensional collapse models the outward radiation force on the accreting material should be significantly stronger than the inward pull of gravity (Larson & Starrfield 1971; Kahn 1974; Wolfire & Cassinelli 1987), in particular when taking dust opacities into account. Since stars as massive as 100–150 M_☉ have been observed (Bonanos et al. 2004; Figer 2005; Rauw et al. 2005) a simple spherically symmetric approach to high-mass star formation is doomed to fail.

Consequently, two different models for massive star formation have been proposed. The first one takes advantage of the fact that high-mass stars always form as members of stellar clusters. If the central density in the cluster is high enough, there is a chance that low-mass protostars collide and so successively build up more massive objects (Bonnell et al. 1998). As the radii of protostars usually are considerably larger than the radii of main-sequence stars in the same mass range, this could be a viable option. However, the stellar densities required to produce very massive stars are still extremely high and seem inconsistent with the observed values of Galactic star clusters (e.g., Portegies Zwart et al. 2010, and references therein). An alternative approach is to argue that high-mass stars form just like low-mass stars by accretion of ambient gas that goes through a rotationally supported disk caused by angular momentum conservation (McKee & Tan 2003). Indeed such disk structures are observed around a number of high-mass protostars (Chini et al. 2004, 2006; Jiang et al. 2008; Davies et al. 2010). Their presence breaks any spherical symmetry that might have been present in the initial cloud and thus solves the opacity problem. Radiation tends to break out along the polar axis, while matter is transported inward through parts of the equatorial plane shielded by the disk. Hydrodynamic simulations in two (Yorke & Sonnhalter 2002; Kuiper et al. 2010a) and three (Krumholz et al. 2009; R. Kuiper et al. 2010, in preparation) dimensions using a flux-limited diffusion approach to the transport of non-ionizing radiation strongly support this picture. Krumholz et al. (2009) find non-axisymmetric Rayleigh–Taylor instabilities using gray radiation transfer, whereas Kuiper et al. (2010b), using frequency-dependent radiation transfer, do not find such instabilities, but nevertheless find strong disk accretion. Strong accretion continues when ionizing radiation is included using ray-tracing methods in three dimensions (Dale et al. 2005, 2007; Peters et al. 2010a, 2010b). In a clustered environment or in individual collapsing cores where the disk becomes gravitationally unstable, the three-dimensional models show that material flows along dense, opaque filaments whereas the radiation escapes through optically thin channels of low-density material. It has been demonstrated that even ionized material can be accreted, if the accretion flow is strong enough. H ii regions are gravitationally trapped at that stage, but soon begin to rapidly fluctuate between trapped and extended states, in agreement with observations (Peters et al. 2010a; R. Galván-Madrid et al. 2010, in preparation). Over time, the same ultracompact H ii region can expand anisotropically, contract again, and take on any of the observed morphological classes (Wood & Churchwell 1989; Kurtz et al. 1994; Peters et al. 2010a, 2010b). In their extended phases, expanding H ii regions drive bipolar neutral outflows characteristic of high-mass star formation (Peters et al. 2010a).

Another key fact that any theory of massive star formation must account for is the apparent presence of an upper mass limit. No star more massive than 100–150 M_☉ has been observed (Massey 2003). This holds for the Galactic field, however, it is also true for young star clusters that are massive enough so that purely random sampling of the initial mass function (IMF; Kroupa 2002; Chabrier 2003) without upper mass limit should have yielded stars above 150 M_☉ (Weidner & Kroupa 2004; Figer 2005; Oey & Clarke 2005; Weidner et al. 2010, see however, Selman & Melnick 2008). This immediately raises the question of what is the physical origin of this apparent mass limit. It has been speculated before that radiative stellar feedback might be responsible for this limit (see, e.g., Zinnecker & Yorke 2007) or alternatively that the internal stability limit of stars with non-zero metallicity lies in this mass regime (Appenzeller 1970a, 1970b, 1987; Baraffe et al. 2001). However, fragmentation could also limit protostellar mass growth. Indeed, this is what we see in the simulations discussed here. The likelihood of fragmentation to occur and the number of fragments to form depends sensitively on the physical conditions in the star-forming cloud and its initial and environmental parameters (see, e.g., Krumholz et al. 2010; Girichidis et al. 2010). Understanding the buildup of massive stars, therefore, requires detailed knowledge about the physical processes that initiate and regulate the formation and dynamical evolution of the molecular clouds these stars form in (Vázquez-Semadeni et al. 2009).

We argue that ionizing radiation, just like its non-ionizing, lower-energy counterpart, cannot shut off the accretion flow onto massive stars. Instead it is the dynamical processes in the gravitationally unstable accretion flow that inevitably occurs during the collapse of high-mass cloud cores that control the mass growth of individual protostars. Accretion onto the central star is shut off by the fragmentation of the disk and the formation of lower-mass companions that intercept inward moving material. We call this process fragmentation-induced starvation and argue that it occurs unavoidably in regions of high-mass star formation where the mass flow onto the disk exceeds the inward transport of matter due to viscosity only and thus renders the disk unstable to fragmentation. We speculate that fragmentation-induced starvation is important not only for present-day star formation but also in the primordial universe during the formation of metal-free Population III stars. Consequently, we expect these stars to be in binary or small number multiple systems and to be of lower mass than usually inferred (Abel et al. 2002; Bromm et al. 2009). Indeed, current numerical simulations provide the first hints that this might be the case (Clark et al. 2008; Turk et al. 2009; Stacy et al. 2010).

In the current study, we analyze the simulations by Peters et al. (2010a) with special focus on the mass growth history of the individual stars forming and the physical processes that influence their accretion rate. We briefly review the numerical method we use and the assumptions and approximations behind it in Section 2. Then we describe our findings in Section 3 and discuss them in the context of present-day and primordial star formation in Section 4. We summarize and conclude in Section 5.

We present three-dimensional, radiation-hydrodynamical simulations of massive star formation that include heating by ionizing and non-ionizing radiation using the adaptive-mesh code FLASH (Fryxell et al. 2000). We use our improved version of the hybrid-characteristics ray-tracing method (Rijkhorst et al. 2006; Peters et al. 2010a) to propagate the radiation on the grid and couple sink particles (Federrath et al. 2010), which we use as models of protostars, to the radiation module via a prestellar model (Peters et al. 2010a).

The simulations start with a 1000 M_☉ molecular cloud. The cloud has a constant density core of ρ = 1.27 × 10⁻²⁰  g cm⁻³ within a radius of r = 0.5 pc and then falls off as r^−3/2 until r = 1.6 pc. The initial temperature of the cloud is T = 30 K. The whole cloud is set up in solid body rotation with an angular velocity ω = 1.5 × 10⁻¹⁴ s⁻¹ corresponding to a ratio of rotational to gravitational energy β = 0.05 and a mean specific angular momentum of j = 1.27 × 10²³ cm² s⁻¹.

We follow the gravitational collapse of the molecular cloud with the adaptive mesh until we reach a cell size of 98 AU. We create sink particles at a cutoff density of ρ_crit = 7 × 10⁻¹⁶ g cm⁻³. All gas within the accretion radius of r_sink = 590 AU above ρ_crit is accreted to the sink particle if it is gravitationally bound to it. The Jeans mass on the highest refinement level is M_jeans = 0.13 M_☉.

The adaptive-mesh technique allows us to resolve the gravitational collapse of the gas from the parsec scale down to a few hundred AU. At higher spatial resolution of only several ten astronomical units, the gas becomes optically thick to non-ionizing radiation, and scattering effects must be taken into account. Since our cutoff density is more than two orders of magnitude smaller than the onset of the optically thick regime at ∼10⁻¹³ g cm⁻³ (e.g., Larson 1969) and we are focusing in our analysis on large-scale effects on the stellar cluster scale, we expect the ray-tracing approximation to be valid. Feedback by radiation pressure, which plays a role on the very small scales and is neglected in our model, is dynamically unimportant on these large scales (Krumholz & Matzner 2009).

We discuss three simulations (see Table 1). In the first simulation (run A), we only allow for the formation of a single sink particle and suppress the formation of secondary sink particles artificially by introducing the dynamical temperature floor

$$\begin{equation} T_\mathrm{min}= \frac{G \mu }{\pi k_\mathrm{B}} \rho (n \Delta x)^2 \end{equation} \tag{ 1 }$$

with Newton's constant G, mean molecular weight μ, Boltzmann's constant k_B, local gas density ρ, and cell size Δx. The temperature floor ensures the sufficient resolution of the Jeans length with n ⩾ 4 cells, which avoids artificial fragmentation (Truelove et al. 1997). The suppression of secondary sink formation guarantees that the accretion flow around the massive star is not weakened by the fragmentation of the disk, which would otherwise inevitably lead to the formation of companion stars that limit the growth of the massive stars in the cluster (see Section 3). In the second simulation (run B), this dynamical temperature floor is not applied, and a whole stellar cluster forms during the course of the simulation. The third simulation (run D) is a control run that allows us to study the influence of radiation feedback on the stellar cluster evolution. As in run B, a small stellar cluster forms, but the stars emit neither ionizing nor non-ionizing radiation.

To demonstrate convergence of our study, we have also run an additional simulation at twice the numerical resolution. This simulation, run D+, which we ran for 0.66 Myr, is identical to run D, except that the cell size at highest grid refinement is 49 AU, the sink particle radius is r_sink = 312 AU and the sink particle cutoff density is ρ_crit = 2.5 × 10⁻¹⁵ g cm⁻³. Since we ran it for only a fraction of the time covered by the other simulations, we will not discuss it at length, but the model ran for long enough to show that the relation between the total stellar cluster mass and the mass of the largest star is the same as in the lower resolution simulations (see Section 4.2.2 for further discussion).

The numerical method along with its inherent physical limitations is discussed in detail in Peters et al. (2010a).

3.1. Accretion History

In this section, we compare the protostellar mass growth rates from our three runs, with a single sink particle (run A), multiple sinks and radiative heating (run B), and multiple sinks with no radiative heating (run D). As already discussed by Peters et al. (2010a), when only the central sink particle is allowed to form (run A), nothing stops the accretion flow to the center. Figure 1 shows that the central protostar grows at a rate $\dot{M} \approx 5.9 \times 10^{-4}\,$M_☉ yr⁻¹ until we stop the calculation when the star has reached 72 M_☉. The growing star ionizes the surrounding gas, raising it to high pressure. However, this hot bubble soon breaks out above and below the disk plane, without affecting the gas flow in the disk⁵ midplane much. In particular, it cannot halt the accretion onto the central star. Similar findings have also been reported from simulations focusing on the effects of non-ionizing radiation acting on somewhat smaller scales (Yorke & Sonnhalter 2002; Krumholz et al. 2007b, 2009; Sigalotti et al. 2009). Radiation pressure is not able to stop accretion onto massive stars and is dynamically unimportant, except maybe in the centers of dense star clusters near the Galactic center (Krumholz & Matzner 2009).

Zoom In Zoom Out Reset image size

The situation is different when the disk can fragment and form multiple sink particles. Initially the mass growth of the central protostar in runs B and D is comparable to the one in run A. As soon as further protostars form in the gravitationally unstable disk, they begin to compete with the central object for accretion of disk material. However, unlike in the classical competitive accretion picture (Bonnell et al. 2001a, 2004), it is not the most massive object that dominates and grows disproportionately fast. Figure 1 shows that, although the accretion rates of the most massive stars (M ⩾ 10 M_☉) steadily decrease, the low-mass stars (M < 10 M_☉) keep accreting at the same rate. Although the detailed mass distribution of the low-mass stars depends on numerical resolution (Federrath et al. 2010), the net effect of their accretion should not, so long as they can form at all. The successive formation of low-mass objects in the disk at increasing radii limits subsequent growth of the more massive objects in the inner disk. Material that moves inward through the disk accretes preferentially onto the sinks at larger radii.

The key to the process, as already discussed by Bate (2000), is that fragments that form later at larger radii will preferentially accrete larger-angular momentum material. This will happen on a much shorter time than the timescales required to redistribute this angular momentum to other parts of the disk at larger radii by viscous or gravitational torques. Similar effects are discussed in the disk fragmentation studies by Kratter et al. (2010).

This behavior is found in models of low-mass protobinary disks, where again the secondary accretes at a higher rate than the primary. Its orbit around the common center of gravity scans larger radii and hence it encounters material that moves inward through the disk before the primary star. This drives the system toward equal masses and circular orbits (Bate & Bonnell 1997; Bate 2000). In our simulations, after a certain transition period hardly any gas makes it all the way to the center and the accretion rate of the first sink particle drops to almost zero. This is the essence of the fragmentation-induced starvation process. In run B, it prevents any star from reaching a mass larger than 25 M_☉. The Jeans mass in run D is smaller than in run B because of the lack of accretion heating, and consequently the highest mass star in run D grows to less than 15 M_☉. In comparison, Krumholz et al. (2007b) described simulations starting with a 10 times less massive core than ours, using both an isothermal equation of state and including accretion heating. In both cases, objects with roughly half the mass of our most massive object formed. The isothermal case formed an only slightly less massive object than the one including accretion heating, just as in our models. The smaller final masses of their objects compared with ours seem to be mainly just a result of their smaller initial core masses.

Inspection of Figure 1 reveals additional aspects of the process. We see that the total mass of the sink particle system increases at a faster rate in the multiple sink simulations, runs B and D, than in the single sink case, run A. This is understandable, because as more and more gas falls onto the disk it becomes more and more unstable to fragmentation, so as time goes by additional sink particles form at larger and larger radii. Star formation occurs in a larger volume of the disk, and mass growth is not limited by the disk's ability to transport matter to its center by gravitational or viscous torques (compare Section 3.2). As a result the overall star formation rate is larger than in run A.

Since the accretion heating raises the Jeans mass and length in run B, the total number of sink particles is higher in run D than in run B, and the stars in run D generally reach a lower mass than in run B (compare Section 4.2.2). These two effects cancel out to lead to the same overall star formation rate for some time. At one point in the evolution, however, also the total accretion rate of run B drops below that of run D. At time t ≈ 0.68 Myr the accretion flow around the most massive star has attenuated below the value required to trap the H ii region. It is able to break out and affect a significant fraction of the disk area. A comparison with the mass growth of run D clearly shows that there is still enough gas available to continue constant cluster growth for another 50 kyr or longer, but the gas can no longer collapse in run B. Instead, it is swept up in a shell surrounding the expanding H ii region.

It is also notable that the expanding ionization front around the most massive stars does not trigger any secondary star formation, which suggests that triggered star formation (Elmegreen & Lada 1977) may not be as efficient as expected, at least on the scales considered here.

Figure 2 shows the individual accretion histories of each of the sink particles in runs B and D. Radiative heating cannot prevent disk fragmentation but raises the local Jeans mass. Hence, as discussed above, much fewer stars form in run B than in run D, and the mass of the most massive stars in run B is higher. It is also evident from the figure that star formation is much more intermittent in the case with radiation feedback (run B). The reason for this behavior is that the star formation process is controlled by the local Jeans mass, which depends to a large degree on how the filaments in the disk shield the radiation (compare Section 3.2), as previously noted by Dale et al. (2005, 2007) and Krumholz et al. (2009). Shielding can lower the Jeans mass temporarily and thereby allow gravitational collapse that would not have occurred otherwise.

Zoom In Zoom Out Reset image size

The accretion histories also reveal large differences between sink particles within the same simulation. For example, in run D a sink particle that forms around t ≈ 0.688 Myr accretes very rapidly and becomes one of the most massive stars in the cluster while other stars only accrete sporadically, interrupted with long intervals of no accretion at all, and some others even accrete only shortly after their formation. There are several reasons for these differences. First, accretion can stop when stars are ejected from the cluster by N-body interactions. This happens mostly to low-mass stars with masses M ≲ 2M_☉. Second, intermediate-mass stars can be starved of material by the formation of companions. This is the essential element of the fragmentation-induced starvation scenario. A star that can grow to intermediate masses is embedded in an accretion flow from a larger gas reservoir. If this reservoir is gravitationally unstable, it can collapse to form companions. This may explain why the intermediate-mass stars in run B seem to have shorter time intervals of no accretion. Since the Jeans mass is higher in this case, the formation of companions is suppressed, and this allows the star to accrete for longer times. Third, the stars are not fixed at a position in the disk plane but move quickly within the disk by gravitational interaction with their surrounding dense gas and with nearby stars. Hence, they are not tied to the filament they form in but can move away from it. If the star moves into a low-density void in the disk, accretion will stop until the star is embedded in higher density gas again. The sink particle in run D mentioned above accretes so vigorously because it moves along with its parental filament over a long time and no companions form around it. Fourth, accretion can stop when the star resides within a low-density H ii region. This applies to the first two massive stars in run B. Their accretion flow had become so weak that the ionizing radiation was able to isolate the stars from the high-density gas in the disk and stopped accretion.

3.2. Disk Growth and Disk Instability

To quantify the accretion of mass onto the disk and the subsequent instability, we calculate the amount of mass in a control volume that encloses the disk at all times. This control volume is the same for all runs and has dimensions 0.24 pc × 0.24 pc × 0.015 pc. We show in Figure 3 the mass M_disk of non-accreted gas contained in the control volume, the mass M_sinks contained in all sink particles as well as the total mass M_tot = M_disk + M_sinks for runs A, B, and D. The disk mass in run A is much larger than in the multiple sink runs B and D since the mass cannot be absorbed by secondary sink particles. Instead, the mass accumulates in the disk plane. Consequently, the star in run A is embedded in a strong accretion flow at all times, which facilitates accretion at an approximately constant rate despite radiation feedback. The evolution in runs B and D differs from the mass growth in run A as soon as secondary sink particles form. The continuous formation of new sink particles in an expanding region around the central star in these runs keeps M_disk almost constant. Figure 3 shows that at t ≈ 0.64 Myr, outflow of material from the central H ii region (Peters et al. 2010a) reduces M_disk even further in run B, so that it falls slightly below the case without feedback (run D) at late times. However, accretion onto secondary sink particles had already cleared the central region, so that this is a relatively minor effect. The accretion onto sink particles is not affected by the central ionization and stays at a constant rate until the H ii region can break free (compare Figure 1).

Zoom In Zoom Out Reset image size

We investigate the initial instability of the disk with an analysis of the Toomre Q-parameter

$$\begin{equation*} Q = \left| \frac{2 c_\mathrm{s} \Omega }{\pi \Sigma G} \right| \end{equation*}$$

with sound speed c_s, angular velocity Ω, surface density Σ, and Newton's constant G. The disk is linearly stable for Q > 1 and linearly unstable for Q < 1 (Toomre 1964; Goldreich & Lynden-Bell 1965). The initial phases of disk instability in run B are displayed in Figure 4. The figure shows slices of density, temperature, and Q in the disk plane for four different times. One can clearly see that the most unstable parts of the disk are the filamentary structures that form as the disk grows in mass. The heating by stellar radiation can suppress instability locally, but shielding by the dense filaments prevents the whole disk from becoming stable and restricts the heating to small regions near the center of the disk that are surrounded by filaments. This shielding makes it possible for star formation to progress radially outward despite accretion heating by the stars. The disk remains sufficiently cool at the inner edge for gravitational instability to set in and star formation proceeds inside–out in the disk plane. Hence, the effect of the filamentary structures in the disk is twofold: they are so dense that they render the disk unstable locally; and, because of their high density, they can effectively shield the radiation from the stellar cluster near the center of the disk, so that radiative heating does not stabilize the outer parts of the disk. Indeed, high-resolution observations of high-mass accretion disk candidates (Beuther et al. 2009) provide some evidence for fragmentation and the presence of substructure on ∼1000 AU scales (as also proposed by Krumholz et al. 2007a).

Zoom In Zoom Out Reset image size

Since the filaments that shield the outer parts of the disk from radiation are optically thick, with an optical depth for non-ionizing radiation of several tens, it is important to estimate the degree to which our simulations, based on a ray-tracing technique to propagate the radiation on the grid, are affected by the lack of heating by diffuse radiation. To test this, we use the adaptive-mesh radiative transfer program RADMC-3D.⁶ RADMC-3D is based on the standard Monte Carlo method of Bjorkman & Wood (2001) in combination with Lucy's method of treating optically thin regions (Lucy 1999). It is the successor of the RADMC code (Dullemond & Dominik 2004) and has been used previously to generate maps of dust emission from these simulations (Peters et al. 2010b). We have calculated self-consistent dust temperatures of the simulation snapshot shown in the last row of Figure 4. Assuming that dust and gas temperatures in the disk plane are equal, we can compare the Monte Carlo dust temperature with the simulation gas temperature, as shown in Figure 5. The comparison demonstrates that direct heating dominates over diffuse radiation in regions more than ∼500 AU away from the stars. The diffusely heated regions lie completely within the region heated in any case by direct radiation, so our ray-tracing method accurately describes the shielding by the filaments of the cold disk region where secondary fragmentation proceeds.

Zoom In Zoom Out Reset image size

To illustrate the tendency of star formation to occur at increasingly larger disk radii, we show the disk radius at which new sink particles form as a function of time for runs B and D in Figure 6. Because the accretion heating sets in already with the very first stars that form in run B, sink formation is suppressed at small disk radii initially. The massive stars slowly spiral outward with time, so that at t ≈ 0.67 Myr their radiation can be shielded by filaments in the disk. Within these filaments, the gas then cools until local collapse sets in and two sink particles form near the center of the disk. Once the filament has dissolved, the gas heats up again and no further sink particles form in the inner disk region.

Zoom In Zoom Out Reset image size

4.1. Initial Conditions

4.2. Relation to Other Theoretical Models

4.2.1. Monolithic Collapse Models

The monolithic collapse model rests on the similarity between the shape of the observed core mass distribution (Motte et al. 1998; Johnstone et al. 2000, 2006; Lada et al. 2008) and the stellar initial mass function, IMF (Kroupa 2002; Chabrier 2003). It assumes a one-to-one relation between the distributions with only a constant efficiency factor separating the two functions. Each molecular cloud core collapses to form a single star or at most a close binary system, with protostellar feedback processes determining the efficiency of the process (Matzner & McKee 2000). This model reduces the problem of the origin of the IMF to the problem of determining the mass spectrum of bound cores, although strictly speaking the idea that the IMF is set by the mass spectrum of cores is independent of any particular model for the origin of that mass spectrum. Arguments to explain the core mass distribution generally rely on the statistical properties of turbulence (Klessen 2001; Padoan & Nordlund 2002; Hennebelle & Chabrier 2008, 2009), which generate structures with a pure power-law mass spectrum. The thermal Jeans mass in the cloud then imposes the flattening and turn-down in the observed mass spectrum.

However, there are a number of caveats. Many of the prestellar cores found in observational surveys appear to be stable entities and thus are unlikely to be in a state of active star formation. In addition, the simple interpretation that one core forms on average one star, and that all cores contain the same number of thermal Jeans masses, leads to a timescale problem (Clark et al. 2007) that requires differences in the core mass function and the IMF. Other concern about this model comes from hydrodynamic simulations, which seem to indicate that massive cores should fragment into many stars rather than collapsing monolithically (Dobbs et al. 2005; Clark & Bonnell 2006; Bonnell & Bate 2006; Federrath et al. 2010). This objection is, however, weakened by the fact that magnetic fields are able to reduce the level of fragmentation on scales of molecular clouds as a whole (Heitsch et al. 2001) as well as of collapsing cloud cores (Hennebelle & Fromang 2008; Hennebelle & Teyssier 2008). In addition, radiative feedback also is able to reduce the number of fragments that form during the collapse of high-mass cores (Krumholz et al. 2007b) as well as low-mass cores (Offner et al. 2009). But again, all current simulations indicate that neither accretion heating nor ionizing radiation can prevent the fragmentation of the massive dense disk that builds up during protostellar collapse, nor can they stop accretion onto the massive star (Yorke & Sonnhalter 2002; Krumholz et al. 2009; Peters et al. 2010a). Instead, the heating merely increases the average mass of the fragments. This is also supported by analytic estimates comparing typical accretion rates onto the disk with its ability to transport matter inward through viscous and gravitational torques (Kratter & Matzner 2006; Kratter et al. 2010). Our current results agree with these studies.

4.2.2. Competitive Accretion

A second model for the origin of the IMF, called competitive accretion, focuses on the interaction between protostars, and between a protostellar population and the gas cloud around it (Bonnell et al. 2001a, 2001b; Bonnell & Bate 2002; Bate & Bonnell 2005). In the competitive accretion picture, the origin of the peak in the IMF is similar to the monolithic collapse model, it is set by the Jeans mass in the prestellar gas cloud. However, rather than fragmentation in the gas phase producing a spectrum of core masses, each of which collapses down to a single star or star system, in the competitive accretion model all gas fragments down to roughly the Jeans mass. Prompt fragmentation therefore creates a mass function that lacks the power-law tail at high masses that we observe in the stellar mass function. This part of the distribution forms via a second phase in which Jeans-mass protostars compete for gas in the center of a dense cluster. The cluster potential channels mass toward the center, so stars that remain in the center grow to large masses, while those that are ejected from the cluster center by N-body interactions remain low mass (Klessen & Burkert 2000; Bonnell et al. 2004). In this model, the apparent similarity between the core and stellar mass functions is an illusion, because the observed cores do not correspond to gravitationally bound structures that will collapse to stars (Clark & Bonnell 2006; Smith et al. 2008).

The competitive accretion picture has been challenged, on the grounds that the kinematic structure observed in star-forming regions is inconsistent with the idea that protostars have time to interact with one another strongly before they completely accrete their parent cores (Krumholz et al. 2005; André et al. 2007).

Taken at face value, competitive accretion models show a correlation between the mass of the most massive star M_max and the total cluster mass M_sinks during the whole cluster evolution that is roughly M_max ∝ M^2/3_sinks (Bonnell et al. 2004). This correlation has been argued to represent a way to observationally confirm competitive accretion (Krumholz & Bonnell 2007) and is in fact in good agreement with observations (Weidner & Kroupa 2006; Weidner et al. 2010). However, we find that our simulations also reproduce the observed relation between M_max and M_sinks.

Figure 7 shows M_sinks as a function of M_max for runs A, B, and D, the higher resolution convergence study run D+ and the relation M_max = 0.39M^2/3_sinks, which was found by Bonnell et al. (2004) as a fit to their simulation data. Over the whole cluster evolution, the curve for run D lies above this fit, while the curve for run B always lies below it. The fit agrees with our simulation data, and with that of Bonnell et al. (2004). The higher resolution simulation run D+ was only carried out for a fraction of the total time of the other simulations, but its data also agree very well with the fit, demonstrating convergence of the results from the lower resolution simulations. This indicates that the scaling is not unique to competitive accretion, but can also be found with the fragmentation-induced starvation scenario and hence cannot be used as an observational confirmation of competitive accretion models.

Zoom In Zoom Out Reset image size

The agreement of both models to this prediction is surprising because the accretion behavior of the most massive stars is totally different. Whereas the massive stars in competitive accretion simulations automatically accrete large fractions of the available gas because they reside in the center of the gravitational potential during the whole cluster evolution, the massive stars in fragmentation-induced starvation models have continuously decreasing accretion rates since they are starved of material by other cluster members and, in the final phase, by feedback from ionizing radiation. It seems that the observed relation between M_sinks and M_max is a very general result of protostellar interaction in a common cluster environment, and not an unambiguous sign of competitive accretion at work.

We can gain further insight from Figure 7. It shows that the accretion heating suppresses low-mass star formation and that for all times the simulation with feedback contains a more massive star relative to the whole cluster mass than the simulation without feedback. When the most massive star in run D reaches 10 M_☉, much more gas is used up to create additional low-mass stars than is funneled toward the most massive one. Thus, the growth of the most massive star in run D is much more ineffective than in run B. Remarkably, all four simulations in Bonnell et al. (2004) show the turn-off toward accretion onto low-mass stars around 10 M_☉, and only one simulation has formed a star more massive than 10 M_☉. A very similar turn-off (indicated by an arrow in Figure 7) at the same mass scale is found in run D, but the radiative heating in run B shifts the turn-off toward higher masses, so that the scaling is followed closely over a larger mass range. In fact, since the turn-off occurs only at the end of the simulation and some of the massive stars are still accreting at this point, it is unknown up to what mass the scaling will continue. However, these results clearly indicate that radiative feedback is necessary to reproduce the scaling for masses beyond 10 M_☉ observed by Weidner et al. (2010).

4.3. Relevance for Primordial Star Formation

We have presented a detailed analysis of the fragmentation-induced starvation scenario, which was first introduced by Peters et al. (2010a). We have studied the accretion history and disk instability in collapse simulations of massive star formation. We have compared the results from a full simulation with radiation feedback and multiple sink particles with two control runs, one with multiple sink particles but without radiation feedback and one with radiation feedback but only a single sink particle. The combined analysis of all three simulations allows us to establish fragmentation-induced starvation as a workable model of massive star formation. We have also compared the new model to the monolithic collapse and competitive accretion models and speculated on the relevance of fragmentation-induced starvation for primordial star formation.

The basic principle of fragmentation-induced starvation is star formation in a rotating gravitationally unstable flow. Gravitational instability in such a flow leads to fragmentation of the accretion flow and the formation of companions around the central massive star that starve the central star of infalling material. In our simulations, star formation occurs within a rotationally flattened, disk-like structure. As more and more material falls onto the disk, star formation proceeds radially outward, keeping the total disk mass roughly constant. The accretion history of individual stars is tightly connected with the position of these stars relative to the dense filaments and low-density voids in the disk. The fact that most if not all high-mass stars are observed in higher order multiple stellar systems that themselves belong to more massive star clusters provides strong evidence for the widespread occurrence of fragmentation that forms the basis of our discussion.

We find that accretion heating does not prevent fragmentation of the disk, but leads to a higher local Jeans mass. As a result, fewer stars form with than without radiation feedback. The accretion heating shifts the masses of the most massive stars up, but leaves the average stellar mass almost unaffected. The accretion heating does not change the overall star formation rate of the whole stellar cluster. Feedback by ionizing radiation is unable to stop protostellar growth if the accretion flow is strong enough. However, if the accretion flow weakens due to starvation, an H ii region can expand and terminate the accretion process. The growing H ii regions reduce the total star formation rate and do not trigger star formation at the ionization front.

Our model is able to explain the observed morphologies (Peters et al. 2010b) and time variability (R. Galván-Madrid et al. 2010, in preparation) of ultracompact H ii regions. We find that we can consistently reproduce the observed relation between the total mass of the star cluster M_sinks and the maximum stellar mass in the cluster M_max, M_max ∝ M^2/3_sinks. This relation seems to be the general outcome of protostellar interaction in a common cluster environment rather than being a signpost of competitive accretion only, as previously claimed. In fact, the dynamical processes discussed here exhibit exactly the opposite behavior of competitive accretion, rather than run-away accretion onto the most massive star together with the suppression of the growth of lower-mass objects, we see that angular momentum conservation and the presence of lower-mass objects limit the mass growth of massive stars.

Our simulations provide evidence for the rejection of proposals that the observed maximum stellar mass of ∼100 M_☉ is set by radiative feedback. When disk fragmentation is artificially suppressed (run A) the central protostar accretes material at very high rate unimpeded by the intense UV radiation it emits without any indications of an upper limit (see also Peters et al. 2010a). When we permit disk fragmentation to occur, it is the process of fragmentation-induced starvation that prevents the stellar mass to become larger than ∼25 M_☉ with our choice of initial conditions. We expect more massive, more centrally condensed, and/or slower rotating cloud cores to lead to more massive protostars. Indeed, extensive parameter studies (Girichidis et al. 2010) show that the initial density profile dominates the accretion behavior, explaining the formation of a 40 M_☉ star from a 100 M_☉ core in Krumholz et al. (2009). The alternative view is to attribute the apparent stellar mass limit to internal stability constraints.

T.P. is a Fellow of the Baden-Württemberg Stiftung funded by their program International Collaboration II (grant P-LS-SPII/18). He also acknowledges support from an Annette Kade Fellowship for his visit to the American Museum of Natural History and a Visiting Scientist Award of the Smithsonian Astrophysical Observatory (SAO). R.S.K. acknowledges financial support from the Baden-Württemberg Stiftung via their program International Collaboration II (grant P-LS-SPII/18) and from the German Bundesministerium für Bildung und Forschung via the ASTRONET project STAR FORMAT (grant 05A09VHA). R.S.K. furthermore gives thanks for subsidies from the Deutsche Forschungsgemeinschaft (DFG) under grant nos. KL 1358/1, KL 1358/4, KL 1359/5, KL 1358/10, and KL 1358/11, as well as from a Frontier grant of Heidelberg University sponsored by the German Excellence Initiative. R.S.K. also thanks the KIPAC at Stanford University and the Department of Astronomy and Astrophysics at the University of California at Santa Cruz for their warm hospitality during a sabbatical stay in spring 2010. M.-M.M.L. was partly supported by NSF grant AST 08-35734. R.B. is funded by the DFG via the Emmy-Noether grant BA 3607/1-1. We acknowledge computing time at the Leibniz-Rechenzentrum in Garching (Germany), the NSF-supported Texas Advanced Computing Center (USA), and at Jülich Supercomputing Centre (Germany). The FLASH code was in part developed by the DOE-supported Alliances Center for Astrophysical Thermonuclear Flashes (ASCI) at the University of Chicago. This work was supported in part by the U.S. Department of Energy contract no. DE-AC-02-76SF00515. We thank the anonymous referee for very useful comments.