Radiation-Hydrodynamic Simulations of Cluster Formation II.

Radiation-Hydrodynamic Simulations of the Formation of Orion-Like Star Clusters II. The Initial Mass Function from Winds, Turbulence, and Radiation


We report a series of simulations of the formation of a star cluster similar to the Orion Nebula Cluster (ONC), including both radiative transfer and protostellar outflows, and starting from both smooth and self-consistently turbulent initial conditions. Each simulation forms stars and brown dwarfs, yielding a stellar mass distribution that ranges from to . We show that a simulation that begins with self-consistently turbulent density and velocity fields embedded in a larger turbulent volume, and that includes protostellar outflows, produces an initial mass function (IMF) that is consistent both with that of the ONC and the Galactic field, at least within the statistical power provided by the number of stars formed in our simulations. This is the first simulation published to date that reproduces the observed IMF in a cluster large enough to contain massive stars, and where the peak of the mass function is determined by a fully self-consistent calculation of gas thermodynamics rather than a hand-imposed equation of state. This simulation also produces a star formation rate that, while still somewhat too high, is much closer to observed values than if we omit either the larger turbulent volume or the outflows. Moreover, we show that the combination of outflows, self-consistently turbulent initial conditions, and turbulence continually fed by motions on scales larger than that of the protocluster yields an IMF that is in agreement with observations and invariant with time, resolving the “overheating” problem in which simulations without these features have an IMF peak that shifts to progressively higher masses over time as more and more of the gas is heated, inconsistent with the observed invariance of the IMF. The simulation that matches the observed IMF also qualitatively reproduces the observed trend of stellar multiplicity strongly increasing with mass. We show that this simulation produces massive stars from distinct massive cores whose properties are consistent with those of observed massive cores. However, the stars formed in these cores also undergo dynamical interactions as they accrete that naturally produce Trapezium-like hierarchical multiple systems of massive stars.

Subject headings:
ISM: clouds — radiative transfer — stars: formation — stars: luminosity function, mass function — turbulence

1. Introduction

The origin of the the stellar initial mass function (IMF) is a classic problem in astrophysics. Since the IMF is most easily measured in young star clusters, and appears to be essentially the same in such clusters and in the field (e.g. Bastian et al., 2010), this problem is closely linked to the problem of how star clusters form. There have been numerous theoretical attacks on these twin problems (see the review by McKee & Ostriker 2007), but a major breakthrough in the past few years has been the realization that the answer is tightly linked to the question of gas thermodynamics. An isothermal gas, even a magnetized one, has no characteristic mass scale (McKee et al., 2010; Krumholz, 2011). This implies that the problem of the origin of the IMF, which is observed to be invariant in both its shape and its characteristic scale, is a separable one.

Models that describe the behavior of a isothermal gas, such as those based on turbulent fragmentation (e.g. Padoan & Nordlund, 2002; Hennebelle & Chabrier, 2008) or competitive accretion (e.g. Bonnell et al., 2001a, b), can predict a shape for the IMF, but in order to determine a characteristic scale must either appeal to additional physics or must define a fiducial “cloud”, whose mean density or other properties (e.g. the normalization of its linewidth-size relation) then determines the location of the IMF peak. In the latter approach, however, it is not clear on what scale one should measure cloud properties: an entire GMC, with a mean density cm obeying the Larson (1981) linewidth-size relation, a massive clump with a mean density cm and a linewidth far above the Larson value (e.g. Shirley et al., 2003), or some other scale? Different choices yield wildly varying characteristic masses. Moreover, cloud properties vary in sufficiently extreme galactic environments, for example showing different linewidth-size relations (e.g. Rosolowsky & Blitz, 2005). Despite this variation, however, there is no evidence for a corresponding variation in the IMF. These problems strongly suggest that the origin of the IMF peak cannot be found in the physics of isothermal gas. Instead, models that seek to explain any characteristic mass scale in the IMF must appeal to departures from isothermality (Rees, 1976; Low & Lynden-Bell, 1976; Spaans & Silk, 2000; Larson, 2005; Krumholz, 2011).

Simulations of star formation mirror these trends depending on the physics they include. Isothermal simulations always produce a characteristic stellar mass that is determined by the initial conditions or the numerical resolution (e.g. Martel et al., 2006), and can always be rescaled to produce an arbitrary stellar mass scale. In contrast, those that include non-isothermality produce characteristic mass scales that are determined by the mechanism that causes them to depart from isothermality, whether it be an imposed equation of state (e.g. Bate & Bonnell, 2005; Jappsen et al., 2005) or the inclusion of radiative transfer, either without (Bate, 2009b, 2012) or with (Krumholz et al., 2007a, 2010, 2011; Offner et al., 2009; Urban et al., 2010) the further step of including stellar radiation. Since comparisons between approximate equations of state and radiative transfer calculations show that the former offer only an extremely poor approximation, progress toward an understanding of the IMF’s characteristic peak therefore requires radiation-hydrodynamic simulations.

The radiation-hydrodynamic simulations that have been conducted thus far have demonstrated several promising features. First, radiation feedback suppresses the formation of brown dwarfs, reproducing the observed turn-down in the IMF at low masses (Bate, 2009b; Offner et al., 2009). Second, simulations including radiation feedback are able to suppress fragmentation in very dense regions, allowing the formation of massive stars when the conditions are approach those seen in real regions of massive star formation (Krumholz et al., 2007a, 2010). Third, radiative simulations produce an IMF that does not vary with the properties of the star-forming cloud in low mass, low-density environments (Bate, 2009b, 2012), nor with the gas metallicity (Myers et al., 2011).

While these results are encouraging, these simulation efforts have for the most part been limited either to single massive cores, or to clouds of low density and/or low mass. For example, Bate (2012) simulates a cloud of mean volume density cm and column density g cm, forming of stars, none larger than . Peters et al. (2010) do form massive stars, but from a cloud with a mean density of cm and a column density of g cm, far below the column density at which radiative effects become important (Krumholz & McKee, 2008; Krumholz et al., 2010) – so low, in fact, as to be optically thin in the near-infrared. In contrast, the mean mass and radius of the star-forming regions studied by Faúndez et al. (2004) implies a volume density cm, a mass of , and a column density of 2 g cm, such that multiple massive stars would be expected, and their radiation would be trapped effectively by the cloud’s high optical depth. Similar Galaxy-wide surveys by Shirley et al. (2003) and Fontani et al. (2005) that target regions of active star formation produce comparable properties. Indeed, the observed cluster mass function is (Lada & Lada, 2003; Fall et al., 2009; Chandar et al., 2010), implying that a majority of stars form in clusters larger than in mass, large enough to possess O stars. The ONC, therefore, is a far more typical star-forming environment than most of the regions explored with radiation-hydrodynamic simulations thus far.

In Krumholz et al. (2011, hereafter Paper I) we reported the first radiation-hydrodynamic simulations to probe this more typical regime of star formation; that calculation followed the collapse of a cloud with a column density of 1 g cm, leading to the production of worth of stars, with an IMF extending from brown dwarfs to O stars. This calculation identified a problem. Radiative suppression of fragmentation, which seems necessary to explain the invariant peak in the IMF and avoid overproduction of brown dwarfs, became too efficient. As the calculation proceeded, the cloud underwent a global collapse, leading to extremely high star formation rates and accretion luminosities. As a result, the gas heated up to the point where further star formation was suppressed. The net result was an IMF that was not invariant, but instead had a peak that moved to systematically higher masses as the calculation proceeded. At early times there were too few massive stars, and at late times too many. Since there is no plausible mechanism to guarantee that all star-forming clouds would stop producing stars at the same point in this evolution, this result was inconsistent with the observed universality of the IMF.

In Paper I, we conjectured that the problem could be resolved by lowering the star formation rate per free-fall time, which would in turn lower the accretion luminosity. Such a change is required by observations even in the absence of the problems rapid star formation creates in the IMF, because observed star formation rates per free-fall time are always a few percent across a very wide range of star-forming environments (Krumholz & Tan, 2007; Evans et al., 2009; Krumholz et al., 2012). In this paper we test that conjecture by performing additional simulations of the formation of ONC-like star clusters, with two extra pieces of physics that should lower the star formation rate per free-fall time. First, rather than simulating an isolated star-forming clump as in Paper I, we embed our initial clump in a larger volume of turbulent gas, and we initialize the simulations such that our clump has self-consistently generated turbulent density and velocity structure. Second, we include protostellar outflows. A number of authors have shown that such outflows can inject significant energy into a star-forming cloud, driving its turbulence and lowering its star formation rate (Li & Nakamura, 2006; Nakamura & Li, 2007; Matzner, 2007; Wang et al., 2010; Cunningham et al., 2011). Ideally, a third piece of physics should also be included: magnetic fields. These both lower the star formation rate by themselves (e.g. Price & Bate, 2009), and also enhance the effectiveness of protostellar outflows (Wang et al., 2010). We plan to do so in future work.

The remainder of this paper is organized as follows. In Section 2 we describe our numerical methods and simulation setup. In Section 3 we report the simulation results, and finally we discuss their implications and draw conclusions in Section 4.

2. Numerical Methods

We simulate star-forming clouds using the orion code, which includes radiative transfer (Howell & Greenough, 2003; Krumholz et al., 2007b; Shestakov & Offner, 2008), hydrodynamics (Klein, 1999), self-gravity (Truelove et al., 1998), accreting sink particles (Krumholz et al., 2004), and a model for protostellar evolution and feedback, including stellar radiation and outflows (Offner et al., 2009; Cunningham et al., 2011). Here we briefly summarize the equations we solve, the code itself, and the initial conditions for the simulations. For the first two of these topics, we refer the reader to Paper I for more details, since the physics included and the numerical methods are identical except where specified below.

2.1. Equations and Algorithms

orion solves the equations of gravito-radiation-hydrodynamics in the two-temperature, mixed-frame flux-limited diffusion approximation. These equations are (Krumholz et al., 2007b)


In these equations, , v, , and are the density, velocity, pressure, and total (thermal plus kinetic) energy density of the gas, is the energy density of the radiation, is the gravitational potential, and are the Planck and Rosseland mean opacities of the dust-plus-gas fluid, is the flux-limiter, is the rate of non-dust cooling (via line and continuum processes in gas at temperatures K where the dust sublimes), and is the proton mass. For more information on the flux-limiter, hot gas cooling rate, and choice of dust opacities, we refer the reader to Paper I.

Terms subscripted by refer to stars; , , and are the position, mass, and momentum of the th star, , , and are the rate at which those stars add or remove mass, momentum, and energy from the gas, is the luminosity of star , and is the weighting kernel that spreads the stellar interaction over some number of computational cells. The equations we solve here differ from those in Paper I in that, in addition to accretion (the terms subscripted with ), we also include protostellar winds (the terms subscripted with ). Stars accrete gas from the computational grid following the sink particle method of Krumholz et al. (2004), and each sink particle is linked to a protostellar evolution code that computes the instantaneous stellar radius and luminosity based on the star’s accretion history, following the method described in the Appendix of Offner et al. (2009). In addition, during each time step, each star returns a portion of the mass it accretes to the grid in the form of a collimated protostellar wind. For details of the numerical implementation, see Cunningham et al. (2011). Our wind parameters are the same as in that paper, i.e. each star ejects a fraction of the gas it accretes (so ), this material is launched with a velocity that of the Keplerian speed at the stellar surface, and the wind gas has a temperature K at launch. It is collimated along the axis defined by the stellar angular momentum vector.

Orion solves Equations (1) – (9) within an overall adaptive mesh refinement (AMR) structure, in which the entire domain is discretized onto a coarse grid of size cells on a side, denoted level 0. Sub-regions within the domain are then covered by progressively finer grids. The grid on level has a resolution a factor of better than that of the coarse grid, and evolves with a time step a factor of smaller. These grids are automatically added and removed on the fly as the calculation proceeds, based on user-specified criteria, up to some pre-specified maximum level .

2.2. Simulation Setup

Name Winds? or
() (pc) (km s) (g cm) (kyr) (pc) (AU)
SmNW No 1000 0.26 2.9 56 1.9 256 5 49
TuNW No 1000 0.46 1.4 23 0.46 256 4 23
TuW Yes 1000 0.46 1.4 23 0.46 256 4 23

Note. – Col. 3: cloud mass. Col. 4: cloud radius (for run SmNW) or box size (for runs TuNW and TuW). Col. 6: mass weighted-mean density at time . Col. 7: free-fall time computed using . Col. 8: size of computational box. Col. 9: number of cells per linear dimension on the coarsest AMR level. Col. 10: finest AMR level. Col. 11: grid resolution on the finest AMR level.

Table 1Simulation Parameters

We compare three different simulations, which we refer to as smooth, no wind (SmNW), turbulent, no wind (TuNW), and turbulent, with winds (TuW); in terms of physics these differ in that the NW simulations have protostellar outflows disabled. We summarize this and other properties of the simulations in Table 1. All simulations consist of a mass of gas with a mean surface density g cm (arranged as described below). Throughout the cloud we set the gas temperature and the radiation energy density to K and erg cm, respectively.

For run SmNW, we use a setup identical to run HR from Paper I (though here we have continued the simulation further in time than we described in that paper), so we only briefly discuss its properties here, and refer readers to Paper I for a fuller description. In run SmNW, the initial gas distribution is a sphere with a radius pc. The density distribution is smooth, and consists of a central core of uniform density that extends to half the cloud’s radius, surrounded by an outer region within which the density falls off with radius as , as suggested by observations of massive clumps (e.g. Sridharan et al., 2005; Beuther et al., 2006). The gas is given an initial turbulent velocity field with a dispersion of km s (one-dimensional), corresponding to an initial virial ratio . The velocity power spectrum is , drawn without imposing any bias in favor of solenoidal or compressive modes following the procedure of Dubinski et al. (1995). Outside the sphere of gas we place a zero-opacity ambient medium with a temperature 100 times larger and a density 100 times smaller than that of the gas at the sphere’s edge. We emphasize that, because the density gradient in the gas only extends to half the initial radius, the overall center to edge density contrast is only a factor of , substantially less than that induced by the turbulent shocks. Thus this initial condition is quite similar to that adopted by other authors who have simulated isolated clouds, e.g. Bonnell et al. (2003) and Bate (2012).7

Figure 1.— Column density distribution in the turbulent initial conditions used for runs TuNW and TuW.

In runs TuNW and TuW we initialize so that, unlike in run SmNW, both the initial density and velocity fields are self-consistently turbulent. We set up a periodic domain of length pc on a side, so that g cm averaged over the box. To initialize the simulation, we impose the same turbulent velocity field as in run SmNW, scaled to a velocity dispersion km s, corresponding to if we use in place of . Although this means the gas is less turbulent initially than in run SmNW, as we see below, damping of the turbulence in run SmNW brings the values closer together as the runs progress. To produce a density field consist with this velocity field, we drive the turbulence and allow the simulation to evolve for two crossing times. During this period we turn off both gravity and radiation, and we hold the gas isothermal at a temperature K by setting the gas ratio of specific heats to ; since, in the absence of stellar sources, molecular cloud gas is close to isothermal, this should be a very good approximation, and ignoring radiation during this setup phase significantly reduces the computational cost. During this setup phase we also fix the computational resolution at cells, with no further refinement. At the end of two crossing times we turn off driving, change the gas ratio of specific heats to , turn on gravity and radiation, and return to our normal refinement criteria (see below). This state represents the initial condition for runs TuNW and TuW. Note that, since the turbulence is driven mostly on large scales, the result of this procedure is essentially a single, dense, turbulent cloud, surrounded by lower density turbulent material; we show this state in Figure 1. This clump is therefore analogous to the isolated one in run SmNW, but is surrounded by a realistic turbulent environment rather than an artificial hot ambient medium.

In all simulations the refinement criteria used to add higher resolution grids are the same. Specifically, we add resolution in any cell that satisfies one of the following three conditions: (1) the density in the cell exceeds the local Jeans density (Truelove et al., 1997), , where is the Jeans number, is the isothermal sound speed, and is the grid spacing on AMR level ; (2) the radiation energy gradient is sharp enough so that (although we sometimes temporarily reduce the coefficient below 0.15 for stability reasons in the TuNW and TuW runs); (3) the cell is within a distance of of any star particle. We refine to a maximum resolution of 49 AU () in run SmNW, and 23 AU () in runs TuNW and TuW. Finally, we note that, while the hydrodynamic and gravitational boundary conditions are necessarily different in the smooth and turbulent runs, the great majority of the star formation in the turbulent runs occurs in subregions much smaller than the entire computational volume, and thus the periodic boundary conditions have minimal impact. We also impose Marshak boundary conditions on the radiation in runs TuNW and TuW in order to let radiation escape the computational volume (c.f. Offner et al. 2009).

3. Results

Before examining the results of our simulations, we first mention two subtleties in the analysis that apply to the remainder of this discussion. First, since we are comparing runs with different initial conditions, it is important to normalize the times so that differences between the runs reflect the underlying physical behavior, and not simply that the dynamical time is different in different cases. Moreover, in the runs with turbulent initial conditions, the strong initial turbulence guarantees that the majority of the mass is compressed into structures that are significantly denser than the volume-averaged density. Given these considerations, the most natural approach, which we adopt, is to measure times in units of the free-fall time evaluated at a density equal to the initial mass-weighted density , since this is the dynamical time appropriate to the bulk of the matter. This approach also has the advantage that it is the most natural basis for observational comparison, since an observation would detect the bulk of the mass, and would be sensitive to the typical density at which this mass resides. For this reason, in what follows whenever we refer to times, we normalize to defined in this manner. We report this quantity in Table 1.8 The second subtlety is that, as in Paper I, we only regard stars as collapsed objects once their mass exceeds , based on one-dimensional calculations of the mass at which second collapse to stellar densities occurs (Masunaga & Inutsuka, 2000). We allow smaller objects to merge with one another and with more massive stars. We therefore restrict our analysis to objects larger than this mass.

3.1. Overall Evolution and Morphology

Figure 2.— Column density (left) and density-weighted mean temperature (right) in run SmNW. The times of each pair of images are indicated in the right column, running from to in steps of 0.25. In the column density plot, white circles indicate the positions of star particles, with the size of the circle indicating the mass of the star. In the right column, the temperature shown is the radiation temperature , defined implicitly by . We show this rather than the gas temperature because the gas and radiation temperatures are nearly equal everywhere in the cloud, except in the hot ambient medium outside the cloud in run SmNW, and in material ejected by protostellar outflows in run TuW. Using the radiation temperature provides a convenient means to filter this contribution.
Figure 3.— Same as Figure 2, but for run TuNW. Note that the color scales are the same, but the size of the region shown is slightly different.
Figure 4.— Same as Figure 3, but for run TuW.

In Figures 2, 3, and 4 we show the large-scale column density and density-weighted temperature distributions in for runs SmNW, TuNW, and TuW, respectively. In all three we observe the same general trend: the turbulence creates an overdense region, which then begins to collapse and form stars. The collapsing structures are filamentary, and the stars are born along the filaments, and particularly at the nodes where the filaments intersect. The temperature is initially small, but as stars form hot spots around individual stars appear, and these gradually spread over time.9

There are a few interesting points to take from these plots. One involves the morphology of the heated regions. In the turbulent runs, even at late times the temperature distribution looks more like a series of islands of heated gas surrounded by a large medium that is either at or quite near to the background temperature. In effect, one can discern something like individual protostellar cores that are heated by the star or star system embedded within them. In contrast, by the end of run SmNW there is simply a single, concentrated region of heating, and one cannot discern individual cores any more. As we show below, this difference proves to be important in determining the evolution of the IMF.

Another interesting point is that the overall morphology is surprisingly similar in runs TuW and TuNW, despite the change in whether we include protostellar winds or not. Partly this is a function of the fact that wind-blown bubbles are fairly low column-density structures, and that we are looking at static slices. In an animation of the column density field, one readily discern outflows driving shells of gas orthogonal to the filaments. However, this clearly has a relatively small effect on the large-scale morphology.

3.2. Star Formation Rate and History

Figure 5.— Total mass in stars (top) and total number of stars (bottom) as a function of time in runs SmNW, TuNW, and TuW.

In Figure 5 we show the star formation history of each of our simulations. The most immediate and striking thing about the Figure is the difference in star formation histories between the smooth and turbulent runs. Run SmNW starts off its star formation more slowly than TuNW or TuW, which is not surprising since its initial density structure is smooth, and possesses no high density peaks that collapse quickly. However, star formation in that run accelerates dramatically as the time approaches . Late in the simulation, the star formation rate approaches . In contrast, in runs TuNW and TuW the star formation rate is roughly constant and fairly low. After a time , only about 10% of the mass has been turned into stars. There is no obvious acceleration with time. Note that, although part of the difference in star formation rates comes from the difference in free-fall times between the turbulent and smooth runs, even if we were to measure in seconds rather than free-fall times run SmNW would have a much larger star formation rate than TuNW or TuW. We summarize the dimensional and dimensionless star formation rates in the simulations in Table 2. Observationally, the dimensionless star formation rate (Krumholz & McKee, 2005)


where the factor of arises because half the cloud mass is above the density used to define , is , with roughly half a dex scatter, over a very wide range of densities and galactic environments (Krumholz & Tan, 2007; Evans et al., 2009; Krumholz et al., 2012).10 Comparing to the table, we see that in runs TuNW and TuW is still roughly an order of magnitude too high compared to observations, but it is roughly an order of magnitude lower than in run SmNW. We discuss the origin of the remaining discrepancy between TuW and the observations further in Section 4.3.

( yr)
SmNW 1.25 0.70 540 16 1.78
TuNW 1.39 0.20 127 7.4 0.33
TuW 1.32 0.15 158 6.2 0.28

Note. – Col. 2: time at which run was stopped. Col. 3: total stellar mass at the end of the run. Col. 4: number of stars present at the end of the run. Col. 5: time-averaged star formation rate in the run, measuring from formation of the first star to the end of the run. Col. 6: dimensionless star formation rate .

Table 2Simulation Outcomes

The difference in star formation rate (SFR) between the runs may be understood readily if we consider what happens to the turbulence, which, in these runs with no magnetic fields, is the main mechanism for regulating the SFR. In run SmNW, the initial turbulence present in the gas decays, and after one crossing time, which is , this decay gives rise to a global collapse and an accelerating star formation rate. In contrast, for runs TuNW and TuW, the box crossing time, and thus the turbulent decay time, is significantly longer than the free-fall time at the mass-weighted mean density. It is comparable to the free-fall time at the volume-weighted mean density, which is much longer. The difference in star formation history between the runs makes a critical point: it matters for their star formation rates that star-forming dense clumps like the one out of which the ONC formed are not isolated objects. They are instead the inner parts of larger turbulent structures, and the energy from those larger scales is able to cascade down to smaller scales and maintain the turbulence for longer than the dynamical times of the small clumps. The turbulent decay timescale in a proto-ONC gas clump is the crossing time of its parent molecular cloud, not the crossing time of the small clump. This point has previously been made by Falceta-Gonçalves & Lazarian (2011) in the context of non-self-gravitating turbulence, and our work strongly confirms their conclusion and extends it to the self-gravitating case.

In contrast, the differences between the two turbulent runs are relatively small. The star formation rate measured by mass (as opposed to number of stars) is lower in TuW than in TuNW. Since our model for protostellar outflows prescribes that 27% of the mass that reaches a star particle (and thus the inner wind launching region) be ejected in an outflow, this reduction in the star formation rate is surprisingly small. This implies that there must be very little entrainment of additional material by the outflows, and even that some of the material that is entrained by outflows must be rapidly stopped and recycled back into the star-forming region. Visual inspection of the morphology of the outflows and accretion flows confirms that this is in fact the case: outflow shocks visible in the animations are generally traveling at right angles to the filaments feeding the stars. This is not an accident. Each star launches its bipolar outflow along the axis specified by its angular momentum vector. If stars are being fed primarily by filaments lying in a plane, as is the case in all our runs, then most of their angular momentum vectors tend to be perpendicular to that plane, producing relatively little entrainment. The minority of outflows that do end up aligning with the filaments possess too little momentum to significantly hinder the accretion flow, and the matter they do eject is stopped by the greater ram pressure of the infalling gas. As a result, it is re-accreted fairly rapidly. Whether this behavior is actually realistic is a separate question, one to which we return in Section 4.

3.3. The IMF

Figure 6.— Evolution of the IMF over time in the three simulations. Thick lines indicate the 50% percentile mass (see main text for formal definition), while the shaded regions indicate the range between the 25th and 75th percentile masses and . Thin lines indicate the mean mass . Colors indicate the run, as described in the legend. Circles long the thick lines indicate the points at which the stellar mass reaches 50 , 100 , 150 , etc. For comparison, thick gray unbroken horizontal lines show , , and for a fully sampled Da Rio et al. (2012) IMF (Equation 13), and the thin gray unbroken horizontal line shows for this IMF. Dashed lines show the equivalent quantities for a Chabrier (2005) IMF.

Another interesting feature in Figure 5 is that the mass in stars in run TuW is smaller than in TuNW at equal times, but that the number of stars is larger in TuW. This indicates an important shift in the stellar IMF between the runs. Although it is less obvious visually from Figure 5, there are also very important differences in the IMF between run SmNW and the other two runs. We now examine these.

For the purposes of quantitative comparison between different simulations, and between simulations and observations, it is helpful to examine percentiles in the cumulative mass distribution function for stars produced in the runs. We define the th percentile mass implicitly via the equation


where is the mass of each individual star, the first sum runs over stars with masses , and the second sum rus over all stars. Thus, for example, is defined by the condition that sum of the masses of all stars smaller than constitutes 50% of the total stellar mass. We also examine the mean stellar mass, defined by


where is the total number of stars. We can measure each of these quantities directly from our simulations at every time. We can also compare the simulation IMFs to observed ones. We select two observational IMFs for comparison. In ONC, Da Rio et al. (2012) find for low mass stars an IMF well-fit by a lognormal function with a width of in , centered on (measured in ; their Table 3). The highest mass bin in Da Rio et al.’s sample is , so to extend this to higher masses we adopt a Chabrier (2003) functional form in which the lognormal at low mass has a powerlaw tail of slope at high mass. Thus the observed ONC IMF to which we compare is


where all masses are in Solar units, over a range from . For this IMF, , , , and . The second comparison IMF is the system IMF of Chabrier (2003, 2005) for the galactic field, which also seems to fit other star clusters reasonably well (Parravano et al., 2011). We use the system rather than the single star IMF because we do not resolve tight binaries. This IMF has the same functional form as Equation (13), but with and . The corresponding percentile and mean values are , , , and . The Chabrier and Da Rio et al. IMFs differ significantly in the number of brown dwarfs and very low mass stars they predict, but converge at masses above a few tenths of . For a discussion of possible origins of the discrepancy between the two IMFs, we refer readers to Da Rio et al. (2012).

In Figure 6 we plot the time evolution of , , , and in each of our simulations, and for the observed Da Rio et al. (2012) and Chabrier (2005) IMFs. The figure immediately reveals some interesting results. First, we see that in run TuW the IMF is in remarkably good agreement with the observed IMFs. At the end of the simulations, the mean mass agrees with the observed Da Rio et al. value to better than 20%, and the 50th percentile mass to less than a factor of 2; we show below that this level of disagreement is consistent with coming simply from statistical sampling variance. Moreover, and perhaps more importantly, the agreement is good at almost all times when there is a significant mass of stars present, because the IMF in run TuW is very stable over time. From , when the total stellar mass reaches 50 , to , when it reaches , we find that , , and stay constant to within ; changes slightly more, almost certainly as a result of under-sampling the high end of the IMF when there are relatively few stars. The change becomes even smaller at later times. From the time when 10% of the mass is in stars () to when 15% is in stars, changes by less than 5% and by less than 10%.

In contrast, for run TuNW the mean mass is relatively stable, but rises systematically with time, increasing by a factor of 2.2 as the stellar mass grows from to , corresponding to times to . This reflects more rapid growth of the more massive stars in the run where winds do not suppress accretion. Moreover, in this run the rate at which new stars form is lower than in run TuW. The agreement with observations in this case is clearly weaker; the run produces an IMF that is too top-heavy.

The changes with time in run TuNW, however, are small compared to those that occur in run SmNW. There and increase by nearly an order of magnitude in a time less than . Each increase in stellar mass of is accompanied by a factor of gain in . This pattern of growth occurs due to the “overheating” problem discussed in Paper I: in run SmNW, star formation is much too rapid and too concentrated, and this produces a rapidly rising accretion luminosity that heats the gas mass to the point where the Bonnor-Ebert mass is too large for stars for new small stars to form. Accretion continues, but it is entirely captured by the existing stellar population, leading to an IMF whose mean and median mass rise with time. Moreover, since all the stars are growing in lockstep the mass distribution in this run is too narrow as well.

Figure 7.— Cumulative mass functions in the simulations, compared to observations. The top set of panels shows the cumulative distribution by mass, and the bottom shows the distribution by number. Within each set of panels, columns show the results from runs SmNW, TuNW, and TuW, as indicated. Rows correspond to times and 1.25, as indicated. In each panel, the colored line indicates the fraction of stellar mass or the fraction of the number of stars in stars with mass less than in the simulation, and the gray band indicates the range from the 10th to 90th percentile resulting from drawing a large number of clusters from the Da Rio et al. (2012) IMF. The hatched band is the 10th to 90th percentile range for a Chabrier (2005) IMF. For details on how this drawing is done, see the Appendix to Paper I. The label in each panel indicates the total mass or number of stars at that time in that simulation.
Figure 8.— Same as Figure 7, but showing differential rather than cumulative mass distributions. The histogram value in each bin shows the total fraction of all stellar mass (for the top panels) or the total fraction of the number of stars (for the bottom panels) falling within that bin. Thick colored lines indicate the simulation result, and gray lines indicate the results of drawing an equal mass stellar population from the Da Rio et al. (2012) IMF. For the gray histogram, the histogram values give the median result, and the vertical lines indicate the range from the 10th to the 90th percentile. We omit the Chabrier (2005) IMF here to reduce clutter.
Figure 9.— Level of statistical agreement between the simulation and observed IMFs as a function of time. At each time, the quantity plotted is the result of a KS test comparing the three simulations (green for SmNW, blue of for TuNW, and red for TuW) to the Da Rio et al. (2012, thick lines) and Chabrier (2005, thin lines) IMFs. The right axis shows the -value returned by the KS test, where is the confidence level at which we can rule out the null hypothesis that the simulation and observed IMFs are drawn from the same underlying distribution. The left axis shows the equivalent confidence level measured in number of standard deviations, which is related by , with erfc the complementary error function. The dashed horizontal black line indicates a confidence level of .

We can also make this comparison more quantitatively. In Figures 7 and 8 we show the cumulative and differential mass distributions produced in our simulations at various times, and compare to observed IMFs. At each time in the simulations, we can quantitatively described the level of consistency or inconsistency between the simulated and observed IMFs using a Kolmogorov-Smirnov (KS) test. We plot the result in Figure 9. Examining the figures, we see that run SmNW is strongly inconsistent with both observed IMFs at most times. At early times the IMF is too bottom-heavy, but as time increases the IMF peak shifts to higher masses. Around run SmNW is fully consistent with the Chabrier IMF, and marginally consistent with the da Rio one, but at later times the IMF peak continues to shift to higher values and becomes inconsistent with both once more. This is the overheating problem described in Paper I. For run TuNW the IMF peak does not shift systematically with time, and so there is no overheating problem. However, the absolute value of the mean mass is systematically too high, as shown in Figure 6. As a result, the overall level of agreement between the simulation and the observed IMFs is poor. On the other hand, run TuW is generally statistically consistent with both the da Rio and Chabrier IMFs at most times.

It is important to add some caveats to this result. First, the KS statistic does not account for the observational and systematic uncertainties in the observed IMFs. Were these uncertainties to be included, it is entirely possible that run TuNW would be consistent within them, and perhaps even that run SmNW would be, at least for a longer period of time. Second, the KS test itself is an imperfect tool. It is most sensitive to differences in distributions near the 50th percentile, and less sensitive to differences on the tail of the distribution. Thus, for example, Figure 8 shows that run SmNW has a slight excess of stars in the range that is clearly visible in a differential mass function on a logarithmic axis. The KS test does not regard this excess as statistically significant, but it is conceivable that a more sensitive statistical test might. Indeed, one can get a sense of the level of statistical power that the KS test provides when applied to our simulations from the fact that, formally, our simulations are consistent with both the da Rio and Chabrier IMFs. This is partly because we are not performing a comparison in the mass range where the two distributions are most different, but it is also partly because, with only 158 stars in run TuW, there is significant sampling noise.

3.4. Gas Thermodynamics

Figure 10.— Phase diagrams of the three runs at different times. The three columns correspond to runs SmNW, TuNW, and TuW, as indicated. The three rows correspond to times , , and 1.25. In each panel, the color indicates the gas mass in a given bin of density and temperature; bins are 0.025 dex wide in both and . The color scale is normalized so that the bin containing the largest amount of mass is 1.0. The long-dashed line indicates the locus in density and temperature at which the code inserts sink particles. The short-dashed lines indicate the locus in density and temperature where the Bonnor-Ebert mass is , , and as indicated. Note that gas in the winds is run TuW is heated to K, well above the temperature range shown here, but there is relatively little mass at these temperatures.

The fragmentation of the gas is driven by its thermodynamics, and we can gain insight into the differences in outcome between the runs by examining the temperature structure of the gas. In Figure 10 we show phase diagrams of the three runs at three different times. Not surprisingly, each of the runs is quite different. First examining run SmNW, we note that, at time , the bulk of the gas in run SmNW is cooler than in the other two runs. This reflects the fact that the total stellar mass in run SmNW is comparable to that in runs TuNW and TuW at this point. Since the free-fall time is longer in run SmNW, this corresponds to a lower total accretion rate and thus a lower accretion luminosity. However, as star formation in run SmNW accelerates, the accretion luminosity rises and the gas heats, while the gas in the other two runs stays relatively cool. Quantitatively, at the final times shown in the bottom row of Figure 10, 42% of all the gas is at temperatures above 50 K in run SmNW (excluding the ambient medium); the equivalent Figures in both runs TuNW and TuW are 7%. It is important to note that this difference is driven by accretion luminosity and not by the intrinsic luminosity of massive stars. If we instead examine run SmNW at time , the most massive star present is 8.8 , smaller than the most massive stars present at time in runs TuNW (13.3 ) and TuW (9.9 ). Nonetheless, we still find that 23% of the mass is at temperatures above 50 K, and 64% is above 30 K, i.e. there is more hot gas in run SmNW even when the individual stars are less massive.

The rapid heating in run SmNW gives rise to the overheating problem identified in Paper I – bulk heating of all the gas makes it impossible for small stars to form, thus shifting the IMF systematically to higher mass as time goes on. Runs TuNW and TuW clearly do not suffer from this problem. Even at late times, the great majority of their gas is at temperatures of no more than K, and there is very little material at temperatures of more than 50 K. Although there clearly is gas being warmed by stars in these runs, there remain pockets of cold, gas at densities g cm and temperatures K that is capable of producing new stars with masses . These are visible in Figure 10, where the phase diagram reveals the presence of material for which the Bonnor-Ebert mass,


is below . In contrast, at late times in run SmNW there no material for which is this small. To be quantitative, at the times shown in the final panel of Figure 10, run SmNW contains only of material in the density and temperature region where , i.e. too little mass to actually create a star. The corresponding figures for runs TuNW and TuW are and , respectively, making it possible for new brown dwarfs to form. It is interesting to note that the amount of cold, high-density gas is generally greater in run TuW than in run TuNW. This is likely an effect of the reduced accretion rate and changed IMF in run TuW compared to TuNW, both of which serve to generally lower the accretion luminosity and thus the heating rate. The spatial distribution of the star formation may also play a role: in run SmNW, because there is no pre-existing density structure at the start of the simulation and because the turbulence decays rapidly, all the stars and gas become concentrated in a single dominant cluster, where stellar heating is very intense. In runs TuNW and TuW, the combination of a pre-existing density structure present in the initial conditions and the non-decay of turbulence throughout the simulation serves to break star formation up into several subclusters, within each of which stellar heating is less intense.

3.5. Massive Cores and Massive Stars

Figure 11.— Images of the initial cores that produced the four most massive stars in simulation TuW. In each column, the upper image shows the column density distribution, centered on the protostars that will grow to be massive stars. The lower image shows the same column density distribution, smeared with a AU Gaussian beam. In the upper panels we indicate the mass of the core (defined as the projected mass within a radius of pc, as indicated by the dashed circles) and the time at which the snapshot is taken. In the lower panels we indicate the core mass that would be inferred from the beam-smeared image and the final mass of the resulting star. Note that the second and fourth columns are nearly identical because two of the final massive stars both form in the same core.
Figure 12.— Same as Figure 11, but for the four stars closest to the median of the final mass distribution.
Figure 13.— Accretion rate versus stellar mass for the four most massive stars present at the end of run TuW. Colored unbroken lines indicate the measured simulation accretion rates, while the dashed black line is the prediction of the McKee & Tan (2003) model (Equation 15). The simulation accretion rates have been smoothed over 500 yr timescales to reduce scatter.

Run TuW is the first published simulation that includes radiative and protostellar outflow feedback, produces an IMF that is in good agreement with the observed IMF over a broad mass range, and forms a large enough cluster for there to be massive stars present. It is therefore important to pay particular attention to the processes by which those massive stars form. We turn now to the properties of the massive stars in run TuW.

There has been considerable discussion in the literature about whether massive stars form from distinct massive protostellar cores (Padoan, 1995; Padoan & Nordlund, 2002; McKee & Tan, 2002, 2003; Krumholz et al., 2007a; Hennebelle & Teyssier, 2008; Hennebelle & Chabrier, 2009), or whether all stars are born from cores with masses , and massive stars subsequently grow from these small seeds by Bondi-Hoyle accretion (Bonnell et al., 1997, 2001a, 2001b, 2004, 2006; Bonnell & Bate, 2002, 2006; Bate & Bonnell, 2005; Smith et al., 2009a, b). A number of authors have also proposed hybrid models, in which massive stars form from gravitationally bound gas structures, but these structures are assembled and fed from larger scales at the same time as they form massive stars (Peretto et al., 2006; Wang et al., 2010). To address this question, we examine the four most massive stars present at the end of run TuW; these have masses of , , , and , respectively, and thus each is large enough that, even if it were to accrete no further, it would be expected to end its life as a supernova. For comparison, we also examine the four stars whose masses are closest to the median mass at the end of the simulation, . For each of these stars, we identify the point in space and time at which that star first appeared in our simulations, and examine the gas density distribution in its vicinity.

We show the results in Figure 11 for the high mass cores and Figure 12 for the low mass cores. To facilitate comparison with observations, in addition to showing the true gas density distribution, we show the distribution smeared with a AU Gaussian beam; we choose this size scale because it is approximately the spatial resolution of the highest published resolution maps of massive cores (e.g. Beuther & Schilke, 2004; Bontemps et al., 2010), though the Atacama Large Millimeter Array (ALMA) will soon produce images at significantly higher resolution. Figure 11 demonstrates that the massive stars in our simulation form in distinct, massive overdensities that can be identified as cores. Their characteristic sizes, determined from visual inspection, are roughly pc. Comparing the gravitational and kinetic energies in this structures shows that they are roughly gravitationally bound and virialized. The flows within them are highly supersonic, producing a filamentary morphology. Nonetheless, these objects are not highly sub-fragmented. There are at most one or two density maxima in each one, not many density maxima. These structures look much like the turbulent cores posited in the McKee & Tan (2003) theory for massive star formation. When smeared on a resolution of 1700 AU, distinct centrally-condensed structures remain visible for three of the four massive stars, indicating that these objects would be detectable as massive cores in an observation.

It is important to understand that our analysis says nothing about the Lagrangian trajectories of the fluid elements that eventually coalesce to form the massive stars in our simulations, a topic that has previously received extensive investigation by Bonnell et al. (2004) and Smith et al. (2009a, b), among others. It may well be that particular fluid elements that are present in the cores at the time shown in Figure 11 do not accrete onto the final star and are instead accreted by other stars or torn off by turbulent motions, while fluid elements not present in the core at the time shown are eventually accreted into the final star. Indeed, McKee & Tan (2003) predicted in their analytic model that turbulent cores should over the course of their lives interact with a surrounding gas mass comparable to that which eventually ends up in their central stars. However, the fact that the Lagrangian elements making up a core change with time is irrelevant to the question of whether, as a massive star forms, it sits at the center of a gravitationally bound Eulerian structure. Figure 11 shows that it does.

We can make the link between the massive cores and the stars they form more quantitative by comparing to the massive core evolution model of McKee & Tan (2002, 2003) and Tan & McKee (2004). This model predicts that the accretion rate onto a star as a function of its mass should be


where is the star’s instantaneous mass, is its final mass, is the surface density of the molecular clump from which it forms, and we have used McKee & Tan’s fiducial parameter choices, with the exception that we have increased the accretion rate by a factor of 2.6 to include subsonic contraction, following Tan & McKee (2004). To evaluate this equation and compare it to our simulations, we take , since this is roughly the mass of our four most massive stars at the end of the simulation. For , we note that, in the simulation, the core is better-defined than the clump, so we adopt McKee & Tan’s result with replaced by . In their fiducial model these agree to within a factor , so this does not significantly affect the accretion rate. As shown in Figure 11, our cores have masses of order in radii of order pc, which corresponds to g cm. With these parameter choices, in Figure 13 we plot the accretion rate as a function of stellar mass for the four most massive stars at the end of the simulation, whose cores are shown in Figure 11, and compare to the McKee & Tan prediction. As the plot shows, the simulation accretion rates agree quite well with the analytic predictions.

In contrast, the cores that give rise to low mass cores (Figure 12) are quite noticeably different from the high mass ones. In three of the four cases (the first, third, and fourth columns in the Figure) they are also centrally-condensed lumps of gas. However, unlike the massive cores they are highly sub-fragmented and show many density maxima. Clearly these objects are not single cores, but instead tightly-packed agglomerations of many smaller cores. For the final low mass core (shown in the second column of Figure 12) the point at which the star forms is a slight overdensity in the middle of a filament, and there is no centrally-concentrated object at all. Thus massive cores and low mass cores have clearly distinct properties. However, we also find that these differences are completely indistinguishable in the smeared images, indicating that it is not possible to distinguish true high mass cores from agglomerations of low mass ones with the resolution available in pre-ALMA telescopes, at least for objects at the kpc distances typical of massive star-forming regions. This conclusion is consistent with that of Offner et al. (2012).

It is important to note that the differences between high and low mass stars is not simply a function of formation time. It is certainly true that the most massive stars at the end of the simulation preferentially began forming early. However, their greater masses are far less a reflection of this than it is of their different formation environments. The four massive stars grow at time-averaged rates of yr, compared to yr for the low mass stars. At the accretion rates typical of the low mass stars, it would require for one of them to grow to the typical of the massive stars. The massive stars are not simply those that form first; they are those that form surrounded by coherent, bound, non-subfragmented structures that provide high accretion rates. This is somewhat similar to the competitive accretion model in that massive stars’ preferred locations at the centers of collapsing regions that provides their high accretion rates. However, it differs from competitive accretion in that these cores are non-subgfragmented and have masses at the same order of magnitude as the final stars, , and therefore intermediate between that of the entire star cluster, and the thermal Jeans mass, . In the competitive accretion model such structures should be absent, because everything fragments down to the thermal Jeans mass (Bonnell et al., 2004; Bate & Bonnell, 2005) and some objects subsequently grow to larger masses by Bondi-Hoyle accretion. There are no objects that do not subfragment in competitive accretion.

Finally, we note that both the high mass and the low mass cores are above the column density threshold g cm for massive star formation posited analytically by Krumholz & McKee (2008) and confirmed numerically by Krumholz et al. (2010). This means that both the high mass and low mass cores have the potential to form massive stars; indeed, in one of the four cases shown in Figure 12, the low mass star is in fact forming in a core that puts most of its mass into a single high mass star. That does not appear to be the case for the other three low mass stars shown in the Figure, however. Thus a high column density is clearly a necessary but not a sufficient condition for massive star formation. A high column density allows radiative heating to suppress the growth of gravitational instabilities that would lead to fragmentation and prevent a massive star from forming. However, if the turbulent density field present before a star begins radiating is already highly non-linearly fragmented, as is the case for several of the low mass cores shown in Figure 12, radiative heating will not undo this fragmentation and prevent the core from forming a small cluster of low mass stars rather than a few massive ones.

3.6. Stellar Multiples

Figure 14.— Multiplicity fraction as a function of system primary mass in run TuW. The thick red line shows as a running average. The light red boxes show computed over discrete bins in . In each case, the width of the box shows the primary mass range for that bin, the asterisk shows the mean multiplicity fraction for stars in that bin, and the vertical extent of the box shows the statistical uncertainty on that value, computed as described in the text. Finally, black crosses indicate observational results, with the horizontal width indicating the mass range for the observations and the vertical range showing the stated uncertainty. The two highest mass observational data points are lower limits, indicated by the upward arrows. The data shown are taken from, from left to right, Basri & Reiners (2006) and Allen (2007) (shown as a single combined point), Fischer & Marcy (1992), Raghavan et al. (2010), Preibisch et al. (1999), and Mason et al. (2009); the data compilation shown here is the same as that in Bate (2012).
Figure 15.— Semi-major axis versus primary star mass for all the binaries in our simulations. For triple and quadruple systems, we plot them only once, showing the properties of the most bound pair of stars. Points are coded by the mass ratio of the system: purple stars for , red squares for , green triangles for , and blue circles for .

It is also illuminating to consider the properties of the stellar multiples that form in run TuW, since producing the correct multiplicity fraction has been proposed as test for star formation models in addition to producing the correct IMF (e.g. Bonnell et al., 2007). We therefore examine the final time slice for this run.11 Extracting the fraction of stars in multiple systems from the simulation requires some care, as pointed out by Bate (2009a). Many of the stars in our simulation form a bound cluster, and thus many stars are bound to many other stars, often in hierarchical structures consisting of dozens of individual stars; for example a binary and a triple system may orbit one another, and these in turn may have additional stars orbiting them. Such agglomerations would be extremely unlikely to survive dynamically even for the lifetime of a massive star, and would break up if we could continue the simulation further.

Thus we follow Bate in defining stellar multiplicity via the following algorithm. We first compute the total energy (gravitational plus kinetic in the center of mass frame) pairwise for each pair of stars in the simulation. We find the most bound system and replace it with a single point mass, with a mass equal to the sum of the two components, a position located at their center of mass, and a momentum equal to the sum of their two momenta. We then continually repeat this process, with the exception that we do not create aggregates consisting of more than four individual stars; should the most bound system contain five our more stars, we proceed to the next most bound pair with fewer than five members instead.12 We terminate the process once there are no more bound pairs consisting of fewer than five individual stars. At the end we are left with a list of star systems, some single and some containing up to four individual stars.

Given this list, we can compute the fraction of multiple systems as a function of primary star mass. For a set of star systems, we define the multiplicity fraction


where , , , and and the numbers of single, binary, triple, and quadruple systems, respectively.13 We choose our sets of star systems in two ways. The first is as a running average; for a primary mass , we compute considering all systems for which the primary mass is within half a dex of . The second is in discrete bins, chosen to roughly match the mass ranges selected in observational surveys. We consider primary mass bins in the range , , , , , , and . In addition to the mean value, we compute the statistical uncertainty in this value for each bin.14

We plot the results in Figure 14. We see that the simulations generally agree well with the observational constraints, with the multiplicity fraction reaching near unity for stars larger than a few , and declining to below for stars smaller than . Our simulations somewhat underproduce binaries at the lowest masses, which is likely a resolution effect, arising because low mass binaries must be very close in order to remain bound, and our simulation resolution makes this difficult. In our simulation we soften particle-particle gravity forces on a scale of 0.25 cells, or AU, and gas-particle gravity forces are necessarily smoothed on the grid scale of 23 AU. Thus we cannot easily make binaries tighter than AU. At this distance the Keplerian speed around a 0.05 object is only 2 km s, comparable to the velocity dispersion in the cluster. Thus low mass binaries formed in our simulation will tend to be broadened and disrupted, and we cannot resolve the tighter binaries that will tend to be hardened. This leads to an artificial reduction in the binary fraction at low masses, a phenomenon also noted by Bate (2012).

In Figure 15 we illustrate some of the properties of our multiple systems. Systems with more massive primaries tend to have the smallest separations, as a result of dynamical hardening and, in some cases, of a companion having been born in the disk of the primary. The companions to the most massive stars also tend to be fairly massive, with mass ratios of ; this is inconsistent with their having been drawn randomly from the IMF. These are often triple or quadruple systems. Thus we see that the massive stars in our simulation tend to form Trapezium-like structures. In contrast, at near-Solar masses, the range of semi-major axes and mass ratios is extremely broad.

3.7. Accretion Variability and Outbursts

Figure 16.— Accretion rate versus time for a sample of 12 randomly-selected stars in run TuW. Thick pale lines show accretion luminosities averaged over 200 yr timescales, while thin darker lines show the accretion luminosity computed over yr timescsles, the finest available given the frequency with which we output simulation data. There is no distinction between red and blue curves; we simply use two different colors to make the two stellar accretion histories shown in each panel more easily distinguishable. Note that the time axes are different for the left and right sides. All times are relative to the instant when a star first appears in the simulations, and plots continue to the end of the simulation.

A final useful datum to be extract from run TuW is the amount of accretion variability for low mass stars, which is of interest for its relevance to the protostellar luminosity problem and the origin of FU Ori outbursts, although due to the duration of our simulation we can only address these issues as they apply to class 0 and I objects, not class II sources. To characterize the degree of luminosity variability we find, we select 12 stars at random at the end of our simulation. The final masses of these stars range from to . For each star we measure its accretion luminosity, which in our code is taken to be


where , , and are the star’s instantaneous mass, accretion rate, and radius, and the factor of accounts for the energy used to drive protostellar outflows (for details see the Appendix of Offner et al., 2009). Our simulation outputs are spaced roughly yr apart, so this represents the minimum timescale on which we can study variability. Outputs are 80 fine grid time steps apart, so this timescale is numerically well resolved.

In Figure 16 we show the accretion history of each of our 12 stars, both at the maximum temporal resolution of the simulations, and smoothed over 200 yr timescales. In the figure, the zero of time is the point at which a given star forms, and we plot the accretion history for the remainder of the simulation. The figure demonstrates several interesting results. First, the majority of stars have relatively smooth luminosity histories when averaged over 200 yr timescales. For only a few examples do order of magnitude variations in the luminosity occur on less than timescales of several kyr. The variability is somewhat larger when measured at the maximum temporal resolution of the simulation, but for most stars this is not a large effect. However, there are three exceptions: the star shown in blue in the upper right panel, the star shown in blue in the lower left panel, and the star shown in red in the lower right panel. All of these stars experience sudden increases in luminosity on timescales below our ability to resolve given the frequency of our output. During these spikes, the luminosity rises by dex compared to the long-term average. These are plausibly FU Ori-type outbursts, although we caution again that these outbursts are occurring in class I sources, not true T Tauri stars.

4. Discussion and Conclusion

4.1. The Role of Protostellar Outflow Feedback

In our simulations, we find that protostellar outflow feedback is not particularly effective. Including outflows reduces the star formation rate by only , comparable to the mass fraction that is ejected from young stars by out subgrid model for protostellar winds. This result at first appears to contradict those of previous studies, including Li & Nakamura (2006), Nakamura & Li (2007), Matzner (2007), Wang et al. (2010), and Cunningham et al. (2011), all of whom find that outflow feedback is important. Our results also contradict those of Hansen et al. (2012), who find that outflow feedback greatly reduces the efficacy of radiative feedback, because it reduces the accretion rate and thus the protostellar luminosity.

Some of our differences from previous results are a function of what effects are included in our simulations, and in previous work. Wang et al. (2010) note that outflow feedback is only effective as a long-term driver of turbulence in the presence of magnetic fields. Fields facilitate transfer of momentum between gas parcels, while in their absence most of the momentum injected into a protocluster by outflows is simply lost, as the outflows break out of the cloud and deposit their momentum and energy outside its boundaries. Since our simulations lack magnetic fields, we likely suffer from a similar underestimate of outflow efficacy.

A second source of difference is likely to be our choice of parameters. We have simulated a fairly massive, high surface density cloud representative of the typical Galactic star-forming region, but with properties quite distinct from those of the star-forming regions closest to the Sun. All of the previous simulations mentioned above have chosen properties typical of these lower density regions. Analytic models suggest that protostellar winds are only able to eject significant mass from clusters with escape velocities below km s (Matzner & McKee, 2000; Matzner, 2007). In run TuW, the escape velocity is km s, within a factor of 2 of the analytic estimate. The comparable figure for the cluster simulated by Wang et al. (2010), which is modeled after NGC 1333, is a factor of two lower: km s. For Hansen et al. (2012), who adopt initial conditions modeled after Ophiuchus, it is km s. Thus it is not surprising that outflows should be much more effective in those simulations than in run TuW. Indeed, placing the cluster masses and surface densities of these simulations on Fall et al. (2010)’s diagnostic diagram for where different sorts of feedback are effective (their Figure 2) immediately predicts this dichotomy.

Given that our simulations likely differ from previous work on the importance of outflow feedback due to both a physical deficiency (lack of magnetic fields) and a choice of parameters that is closer to the typical region than most previous work, it is hard to draw general conclusions about the importance of outflow feedback in regulating star formation. Resolution of this question will have to await future magnetohydrodynamic simulations that probe the higher density regime we have explored in this work.

4.2. Implications for Massive Star Formation

The picture of massive star formation that emerges from our simulations is generally consistent with the turbulent core model proposed by McKee & Tan (2002, 2003). The massive stars form at the centers of well-defined, turbulent, centrally-concentrated structures, and these structures feed mass onto them at a rate that is consistent with the predictions of the McKee & Tan model. In contrast, the regions from which low mass stars form are quite noticeably different. They are either messy regions consisting of many small density peaks and no clear central concentration, or they are small regions of filaments. Thus the basic core to star mapping proposed in the turbulent core model appears to describe our simulation fairly well.

However, we also do see elements of the alternative competitive accretion model (Bonnell et al., 2007, and references therein) operating as well. In particular, our massive stars do all form as part of small sub-clusters and experience significant dynamical interactions. These interactions appear to be important in shaping the multiplicity properties of the resulting stars, and in producing the Trapezium-like systems in which most of our massive stars find themselves at the end of the simulation.

4.3. Implications for the IMF and the Star Formation Rate

Run TuW represents the first simulation published to date that reproduces the observed IMF in a cluster like the ONC that is large enough to contain massive stars, and where the peak of the mass function is determined by a fully self-consistent calculation of gas thermodynamics. Previous simulations that have had success in reproducing the IMF have either examined small, low-density star-forming regions that would not be expected to produce massive stars (Offner et al., 2009; Bate, 2009b, 2012; Urban et al., 2010; Hansen et al., 2012), or have relied on a parameterized, non-self-consistent equation of state to determine the location of the IMF peak (e.g. Bate & Bonnell, 2005; Jappsen et al., 2005).

The success of run TuW, in contrast to the failures of runs SmNW and TuNW, suggests that obtaining the correct IMF from a self-consistent simulation of a typical star-forming environment is not as simple as some authors have posited (e.g. Bonnell et al., 2007). While a lognormal function with a powerlaw tail at high masses appears to be a fairly generic result regardless of what physics is included in the simulations, getting the peak of the lognormal to lie at the correct position in a calculation where it is determined self-consistently rather than through a hand-imposed equation of state requires careful attention to the thermodynamics of the gas, which is in turn determined primarily by the accretion luminosity of protostars.

This requires several ingredients to work correctly. As conjectured in Paper I, the star formation rate cannot be too high, and star formation cannot become too centrally cocentrated; if it is, the resulting accretion luminosity becomes so high that formation of low mass stars is suppressed and the IMF peak marches to ever-increasing mass with time. In addition, outflows appear to make a small but significant contribution by both reducing the masses of individual accreting protostars, reducing the accretion luminosity, and ejecting mass from the warm regions near accreting stars, increasing fragmentation. The combination of these effects shifts the mean mass downward by a factor of . Our results suggest that future simulations of gas collapse and fragmentation, if they are to reproduce the observed IMF while treating the gas thermodynamics self-consistently, must at a minimum include these ingredients.

We obtain a reduction in the rate and degree of concentration of star formation in run TuW mainly because we have starting our simulations with a self-consistent density and velocity field and by embedding the protocluster gas clump in a realistic surrounding turbulent molecular cloud that continually feeds energy into it via a turbulent cascade, rather than simulating a smooth, isolated clump as in most previous work. However, the star formation rate per free-fall time we obtain is still an order of magnitude too high compared to observations, likely as a result of other mechanisms we have not included. For example, Price & Bate (2009) and Padoan & Nordlund (2011) both find that magnetic fields with strengths comparable to those in observed molecular clouds reduce star formation rates by a factor of a few to . Depending on protocluster properties, ionized gas and radiation pressure may also contribute to reducing the star formation rate (e.g. Fall et al., 2010; Lopez et al., 2011). Ultimately, since accretion luminosity plays a critical role in regulating the IMF, the problems of determining the star formation rate and the IMF cannot be fully separated. A truly accurate simulation must reproduce both.

This work was supported by an Alfred P. Sloan Fellowship (MRK); the NSF through grants CAREER-0955300 (MRK) and AST-0908553 (CFM and RIK); NASA through ATFP grant NNX09AK31G (RIK, CFM, and MRK) and a Chandra Space Telescope grant (MRK); and the US Department of Energy at LLNL under contrast DE-AC52-07NA (RIK). Support for computer simulations was provided by an LRAC grant from the NSF through Teragrid resources and NASA through grants from the ATFP and Chandra Programs.


  1. affiliation: Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064; krumholz@ucolick.org
  2. affiliation: Lawrence Livermore National Laboratory, P.O. Box 808, L-23, Livermore, CA 94550
  3. affiliationmark:
  4. affiliation: Department of Astronomy and Astrophysics, University of California, Berkeley, Berkeley, CA 94720
  5. affiliation: Department of Physics, University of California, Berkeley, Berkeley, CA 94720
  6. slugcomment: Accepted for publication in The Astrophysical Journal
  7. An additional difference between our setup and that of Bonnell et al. (2003) and Bate (2012) is that we place an ambient medium outside our cloud that is in thermal pressure balance with the material at the cloud edge, while the smoothed particle hydrodynamics (SPH) simulations of Bonnell et al. and Bate have a vacuum outside their clouds. However, this difference is almost certainly negligible. The thermal pressure of our ambient medium is set equal to the thermal pressure of the cloud, which is smaller than either the ram pressure or the self-gravitational weight of the cloud by a factor of . Thus the extra pressure provided by the external medium will enhance the collapse that would occur due to gravity alone by only . Even this is likely an overestimate of the difference between the two simulation methods, because, while formally the SPH simulations have vacuum outside their clouds, SPH creates an artificial surface tension at density discontinuities (Price, 2008), and this will act very much like a confining external pressure. Our Eulerian simulation method does not suffer from this problem.
  8. Note that in Paper I we instead used the volume weighted mean density to compute the free-fall time for run SmNW; however, because the initial density field is very smooth, the difference between volume- and mass-weighted mean density free-fall times for this run is only .
  9. In runs TuNW and TuW there are sometimes brief increases in the overall background temperature level visible in some snapshots, but these are short-lived and small, generally keeping the temperature K. These flashes are associated with brief increases in the accretion luminosity that are large enough to heat the entire simulation volume above 10 K for short periods.
  10. There is some subtlety in the observational comparison here, because real observations usually have an upper limit on the density to which they are sensitive, for example because the tracer being used depletes or becomes very optically thick at high density. Since we include all the mass above in our computation of , we are not capturing this effect. However, the change in mass it would induce is small, because both real star-forming clouds and our simulated turbulent clouds have density probability distribution functions that are sharply declining at densities above the peak. Thus the amount of mass missed due to the density upper limit in the observations is likely to be very small.
  11. In this section of the paper alone we do not exclude stars smaller than from consideration, but we consider them only as companions to larger stars. We allow them to count in this capacity because to omit them would artificially make it impossible for stars near our cut to have companions.
  12. The results are not particularly sensitive to the choice of four as the maximum size of a system, as long as we stop at some point well short of allowing the entire cluster to be considered a single large star system.
  13. Following Bate (2009a) and Hubber & Whitworth (2005), we measure this quantity rather than either the companion star fraction or the fraction of stars in multiple systems because it is more robustly determined observationally. If a new member of a multiple system is found, for example leading to a binary being reclassified as a triple, does not change, while the companion star fraction and the fraction of stars in multiple systems does.
  14. We determine the statistical uncertainty by assuming that there is a true multiplicity fraction for stars in that bin, and that our sample of systems in that bin represent a series of random drawings that follow a binomial distribution. From these assumptions, we can compute the probability distribution for given the measured multiplicity fraction in the simulations and the number of systems in that mass bin. We compute the uncertainty by finding the range of values for that enclose the central of the probability distribution. Note that this range is not in general symmetric about .


  1. Allen, P. R. 2007, ApJ, 668, 492
  2. Basri, G., & Reiners, A. 2006, AJ, 132, 663
  3. Bastian, N., Covey, K. R., & Meyer, M. R. 2010, ARA&A, 48, 339
  4. Bate, M. R. 2009a, MNRAS, 392, 590
  5. —. 2009b, MNRAS, 392, 1363
  6. —. 2012, MNRAS, 419, 3115
  7. Bate, M. R., & Bonnell, I. A. 2005, MNRAS, 356, 1201
  8. Beuther, H., & Schilke, P. 2004, Science, 303, 1167
  9. Beuther, H., Zhang, Q., Sridharan, T. K., Lee, C.-F., & Zapata, L. A. 2006, A&A, 454, 221
  10. Bonnell, I. A., & Bate, M. R. 2002, MNRAS, 336, 659
  11. —. 2006, MNRAS, 370, 488
  12. Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 1997, MNRAS, 285, 201
  13. —. 2001a, MNRAS, 323, 785
  14. Bonnell, I. A., Bate, M. R., & Vine, S. G. 2003, MNRAS, 343, 413
  15. Bonnell, I. A., Clarke, C. J., & Bate, M. R. 2006, MNRAS, 368, 1296
  16. Bonnell, I. A., Clarke, C. J., Bate, M. R., & Pringle, J. E. 2001b, MNRAS, 324, 573
  17. Bonnell, I. A., Larson, R. B., & Zinnecker, H. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson: University of Arizona Press), 149–164
  18. Bonnell, I. A., Vine, S. G., & Bate, M. R. 2004, MNRAS, 349, 735
  19. Bontemps, S., Motte, F., Csengeri, T., & Schneider, N. 2010, A&A, 524, A18+
  20. Chabrier, G. 2003, PASP, 115, 763
  21. Chabrier, G. 2005, in Astrophysics and Space Science Library, Vol. 327, The Initial Mass Function 50 Years Later, ed. E. Corbelli, F. Palla, & H. Zinnecker, 41–+
  22. Chandar, R., Fall, S. M., & Whitmore, B. C. 2010, ApJ, 711, 1263
  23. Cunningham, A. J., Klein, R. I., Krumholz, M. R., & McKee, C. F. 2011, ApJ, 740, 107
  24. Da Rio, N., Robberto, M., Hillenbrand, L. A., Henning, T., & Stassun, K. G. 2012, ApJ, 748, 14
  25. Dubinski, J., Narayan, R., & Phillips, T. G. 1995, ApJ, 448, 226
  26. Evans, N. J., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321
  27. Falceta-Gonçalves, D., & Lazarian, A. 2011, ApJ, 735, 99
  28. Fall, S. M., Chandar, R., & Whitmore, B. C. 2009, ApJ, 704, 453
  29. Fall, S. M., Krumholz, M. R., & Matzner, C. D. 2010, ApJ, 710, L142
  30. Faúndez, S., Bronfman, L., Garay, G., et al. 2004, A&A, 426, 97
  31. Fischer, D. A., & Marcy, G. W. 1992, ApJ, 396, 178
  32. Fontani, F., Beltrán, M. T., Brand, J., et al. 2005, A&A, 432, 921
  33. Hansen, C. E., Klein, R. I., McKee, C. F., & Fisher, R. T. 2012, ApJ, 747, 22
  34. Hennebelle, P., & Chabrier, G. 2008, ApJ, 684, 395
  35. —. 2009, ApJ, 702, 1428
  36. Hennebelle, P., & Teyssier, R. 2008, A&A, 477, 25
  37. Howell, L. H., & Greenough, J. A. 2003, JCP, 184, 53
  38. Hubber, D. A., & Whitworth, A. P. 2005, A&A, 473, 113
  39. Jappsen, A.-K., Klessen, R. S., Larson, R. B., Li, Y., & Mac Low, M.-M. 2005, A&A, 435, 611
  40. Klein, R. I. 1999, J. Comp. App. Math., 109, 123
  41. Krumholz, M. R. 2011, ApJ, 743, 110
  42. Krumholz, M. R., Cunningham, A. J., Klein, R. I., & McKee, C. F. 2010, ApJ, 713, 1120
  43. Krumholz, M. R., Dekel, A., & McKee, C. F. 2012, ApJ, 745, 69
  44. Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007a, ApJ, 656, 959
  45. —. 2011, ApJ, 740, 74
  46. Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J. 2007b, ApJ, 667, 626
  47. Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250
  48. —. 2008, Nature, 451, 1082
  49. Krumholz, M. R., McKee, C. F., & Klein, R. I. 2004, ApJ, 611, 399
  50. Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304
  51. Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57
  52. Larson, R. B. 1981, MNRAS, 194, 809
  53. —. 2005, MNRAS, 359, 211
  54. Li, Z.-Y., & Nakamura, F. 2006, ApJ, 640, L187
  55. Lopez, L. A., Krumholz, M. R., Bolatto, A. D., Prochaska, J. X., & Ramirez-Ruiz, E. 2011, ApJ, 731, 91
  56. Low, C., & Lynden-Bell, D. 1976, MNRAS, 176, 367
  57. Martel, H., Evans, II, N. J., & Shapiro, P. R. 2006, ApJS, 163, 122
  58. Mason, B. D., Hartkopf, W. I., Gies, D. R., Henry, T. J., & Helsel, J. W. 2009, AJ, 137, 3358
  59. Masunaga, H., & Inutsuka, S. 2000, ApJ, 531, 350
  60. Matzner, C. D. 2007, ApJ, 659, 1394
  61. Matzner, C. D., & McKee, C. F. 2000, ApJ, 545, 364
  62. McKee, C. F., Li, P. S., & Klein, R. I. 2010, ApJ, 720, 1612
  63. McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565
  64. McKee, C. F., & Tan, J. C. 2002, Nature, 416, 59
  65. —. 2003, ApJ, 585, 850
  66. Myers, A. T., Krumholz, M. R., Klein, R. I., & McKee, C. F. 2011, ApJ, 735, 49
  67. Nakamura, F., & Li, Z.-Y. 2007, ApJ, 662, 395
  68. Offner, S. S. R., Capodilupo, J., Schnee, S., & Goodman, A. A. 2012, MNRAS, 420, L53
  69. Offner, S. S. R., Klein, R. I., McKee, C. F., & Krumholz, M. R. 2009, ApJ, 703, 131
  70. Padoan, P. 1995, MNRAS, 277, 377
  71. Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870
  72. —. 2011, ApJ, 730, 40
  73. Parravano, A., McKee, C. F., & Hollenbach, D. J. 2011, ApJ, 726, 27
  74. Peretto, N., André, P., & Belloche, A. 2006, A&A, 445, 979
  75. Peters, T., Banerjee, R., Klessen, R. S., et al. 2010, ApJ, 711, 1017
  76. Preibisch, T., Balega, Y., Hofmann, K.-H., Weigelt, G., & Zinnecker, H. 1999, New A, 4, 531
  77. Price, D. J. 2008, Journal of Computational Physics, 227, 10040
  78. Price, D. J., & Bate, M. R. 2009, MNRAS, 398, 33
  79. Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1
  80. Rees, M. J. 1976, MNRAS, 176, 483
  81. Rosolowsky, E., & Blitz, L. 2005, ApJ, 623, 826
  82. Shestakov, A. I., & Offner, S. S. R. 2008, JCP, 227, 2154
  83. Shirley, Y. L., Evans, II, N. J., Young, K. E., Knez, C., & Jaffe, D. T. 2003, ApJS, 149, 375
  84. Smith, R. J., Clark, P. C., & Bonnell, I. A. 2009a, MNRAS, 396, 830
  85. Smith, R. J., Longmore, S., & Bonnell, I. 2009b, MNRAS, 400, 1775
  86. Spaans, M., & Silk, J. 2000, ApJ, 538, 115
  87. Sridharan, T. K., Beuther, H., Saito, M., Wyrowski, F., & Schilke, P. 2005, ApJ, 634, L57
  88. Tan, J. C., & McKee, C. F. 2004, ApJ, 603, 383
  89. Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179+
  90. —. 1998, ApJ, 495, 821
  91. Urban, A., Martel, H., & Evans, N. J. 2010, ApJ, 710, 1343
  92. Wang, P., Li, Z., Abel, T., & Nakamura, F. 2010, ApJ, 709, 27
This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.