RADIATION FEEDBACK, FRAGMENTATION, AND THE ENVIRONMENTAL DEPENDENCE OF THE INITIAL MASS FUNCTION

Our numerical method is identical to that of Krumholz et al. (2009), and we give an extensive description in the supporting online material of that paper, so we only summarize our methods briefly here. Our simulations use the parallel adaptive mesh refinement (AMR) radiation-hydrodynamics code ORION. In ORION, the gas plus radiation fluid is represented at every grid point by a vector of conserved quantities (ρ, ρv, ρe, E), where ρ is the density, v is the velocity, e is the specific non-gravitational energy (including kinetic and thermal), and E is the radiation energy density. The computational domain also contains an arbitrary number of point mass "star" particles, each of which is characterized by a position $\textbf {\it x}$, a mass M, a momentum $\textbf {\it p}$, and a luminosity L. The code updates these quantities by solving the equations of radiation hydrodynamics plus gravity in the conservative, mixed-frame form (Mihalas & Klein 1982), retaining terms to order v/c accuracy, and using the flux-limited diffusion approximation to represent the radiation flux (Krumholz et al. 2007b). The equations for the gas are

$$\begin{equation} \frac{\partial }{\partial t}\rho = - \nabla \cdot (\rho v) - \sum _i \dot{M}_i W(\textbf {\it x}-\textbf {\it x}_i) \end{equation} \tag{ 1 }$$

$$\begin{eqnarray} \frac{\partial }{\partial t}(\rho v) & = & -\nabla \cdot (\rho vv) - \nabla P - \rho \nabla \phi - \lambda \nabla E \nonumber \\ & & {} - \sum _i \dot{\textbf {\it p}}_i W(\textbf {\it x}-\textbf {\it x}_i) \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \frac{\partial }{\partial t}(\rho e) & = & -\nabla \cdot [(\rho e+P)v] - \rho v\cdot \nabla \phi - \kappa _{\rm 0P} \rho (4 \pi B - c E) \nonumber \\ & & {} + \lambda \left(2 \frac{\kappa _{\rm 0P}}{\kappa _{\rm 0R}} - 1\right) v\cdot \nabla E - \sum _i \dot{\mathcal {E}}_i W(\textbf {\it x}- \textbf {\it x}_i) \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \frac{\partial }{\partial t}E & = & \nabla \cdot \left(\frac{c\lambda }{\kappa _{\rm 0R} \rho } \nabla E\right) + \kappa _{\rm 0P} \rho (4 \pi B - c E) \nonumber \\ & & {} - \lambda \left(2\frac{\kappa _{\rm 0P}}{\kappa _{\rm 0R}} - 1\right) v\cdot \nabla E - \nabla \cdot \left(\frac{3 - R_2}{2} vE\right) \nonumber \\ & & {} + \sum _i L_i W(\textbf {\it x}- \textbf {\it x}_i), \end{eqnarray} \tag{ 4 }$$

where the summations run over all star particles present; $\dot{M}_i$, $\dot{\textbf {\it p}}_i$, and $\dot{\mathcal {E}}_i$ are the rates at which mass, momentum, and energy are transferred from gas to star particles; and $W(\textbf {\it x}-\textbf {\it x}_i)$ is a weighting function that distributes the transfer over a kernel 4 cells in radius. We calculate these using the Krumholz et al. (2004) sink particle algorithm, which we summarize below. The corresponding evolution equations for the star particles are

$$\begin{equation} \frac{d}{dt} M_i = \dot{M}_i \end{equation} \tag{ 5 }$$

$$\begin{equation} \frac{d}{dt} \textbf {\it x}_i = \frac{\textbf {\it p}_i}{M_i} \end{equation} \tag{ 6 }$$

$$\begin{equation} \frac{d}{dt} \textbf {\it p}_i = -M_i \nabla \phi + \dot{\textbf {\it p}}_i. \end{equation} \tag{ 7 }$$

In these equations, the gravitational potential ϕ is given by

$$\begin{equation} \nabla ^2\phi = -4\pi G \left[ \rho + \sum _i M_i \delta (\textbf {\it x}-\textbf {\it x}_i)\right]. \end{equation} \tag{ 8 }$$

The pressure P is given by

$$\begin{equation} P = \frac{\rho k_B T_g}{\mu m_{\rm H}} = (\gamma -1) \rho \left(e - \frac{v^2}{2}\right), \end{equation} \tag{ 9 }$$

where T_g is the gas temperature, μ = 2.33 is the mean molecular weight for molecular gas of solar composition, and γ is the ratio of specific heats. We adopt γ = 5/3, appropriate for gas too cool for hydrogen to be rotationally excited, but this choice is essentially irrelevant because T_g is set almost purely by radiative effects. The remaining quantities are the comoving frame specific Planck- and Rosseland-mean opacities κ_0R and κ_0P, the Planck function B = ca_RT⁴_g/(4π), and the flux limiter λ and Eddington factor R₂, computed using the Levermore & Pomraning (1981) approximation:

$$\begin{eqnarray} \lambda & = & \frac{1}{R} \left(\mbox{coth} R - \frac{1}{R}\right) \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} R & = & \frac{|\nabla E|}{\kappa _{\rm 0R} \rho E} \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} R_2 & = & \lambda + \lambda ^2 R^2. \end{eqnarray} \tag{ 12 }$$

We obtain the dust opacities κ_0P and κ_0R from a piecewise-linear fit to the models of Pollack et al. (1994); see Krumholz et al. (2009) for the exact functional form. It is worth noting that the Rosseland opacity we use includes absorption but not scattering effects, and as a result is likely something of an underestimate. Using a higher opacity would likely enhance the radiation effect we describe below, as suggested by recent static radiative transfer calculations (Dunham et al. 2010).

We solve these equations in four steps. First, the hydrodynamics module updates Equations (1)–(3) using all the terms on the right-hand side except those involving star particles or radiation. The update is based on a conservative Godunov scheme with an approximate Riemann solver, and is second-order accurate in time and space (Truelove et al. 1998; Klein 1999). Second, the gravity module solves the Poisson Equation (8) to update the gravitational potential using a multigrid iteration scheme (Truelove et al. 1998; Klein 1999; Fisher 2002). Third, the radiation module updates the right-hand sides of Equations (2)–(4) for the terms involving radiation. The module uses the Krumholz et al. (2007b) operator splitting method, in which the dominant terms describing radiation-gas energy exchange and radiation diffusion are updated using a fully implicit method based on pseudo-transient continuation (Shestakov & Offner 2008), and then the sub-dominant work and advection terms are handled explicitly. The fourth step is the star particle module. This portion of the code updates Equations (1)–(7) using the gas-particle exchange terms. The code computes $\dot{M}_i$ by fitting the flow in the vicinity of each particle to a Bondi–Hoyle flow, and then sets $\textbf {\it p}_i$ and $\mathcal {E}_i$ by requiring that, in the frame comoving with the particle, accretion does not alter the radial velocity, angular momentum, or temperature of the gas (Krumholz et al. 2004). The method then updates the luminosity and the internal state of each star based on a simple one-zone protostellar evolution model. Details of the model are given in the Appendices of Offner et al. (2009b).

All of these modules operate with an AMR framework (Berger & Oliger 1984; Berger & Collela 1989; Bell et al. 1994) in which we discretize the computational domain onto a series of levels l = 0, 1, 2, ..., L. The coarsest level is l = 0, which covers the entire computational domain. Subsequent levels cover sub-regions of the computational domain which are described by the union of rectangular grids. Levels are nested such that any grid on level l > 0 must be fully contained within one or more grids of level l − 1. The cell size on level l is Δx_l, and the spacings on levels l > 0 are related to that on level 0 by Δx_l = Δx₀/2^l. The process of advancing the computation through these levels is recursive. We first advance the grids on level 0 by a time Δt₀, and then we advance the grids on level 1 by two time steps of size Δt₁ = Δt₀/2. However, each of these advances is followed by two advances on level 2, and so forth, such that each level 0 advance involves 2^l advances of time step Δt₀/2^l on level l. At the end of each advance of level l > 0 through two steps, we perform a synchronization procedure between levels l and l − 1 to ensure conservation of mass, momentum, and energy across the interface between the two levels. We set the overall time step Δt₀ by computing the Courant condition (including a contribution to the effective signal speed from radiation pressure; Krumholz et al. 2007b) separately on each level at the beginning of each coarse time step. We then set Δt₀ = min(2^lΔt^l), ensuring that each level obeys the Courant condition for its advances.

The ORION AMR framework automatically adds and removes higher-resolution grids throughout a simulation. We determine when higher-resolution grids are required based on the following criteria:

• 1.  
  Any cell with a density greater than half the initial density at the edge of the cloud (see below) must be refined to at least level 1.
• 2.  
  We refine any cell within whose distance d to the nearest star particles is less than 16Δx_l. This ensures that regions around stars are always refined to the maximum level.
• 3.  
  We refine any cell where the density exceeds the Jeans density, given by

  $$\begin{equation} \rho _J = J^2 \frac{\pi c_s^2}{G \Delta x_l^2}, \end{equation} \tag{ 13 }$$

  where c_s is the sound speed and we use J = 1/16. This avoids artificial fragmentation (Truelove et al. 1997).
• 4.  
  We refine any cell where the gradient in the radiation energy density satisfies

  $$\begin{equation} |\nabla E| > 0.15 \frac{E}{\Delta x_l}. \end{equation} \tag{ 14 }$$

  This ensures that we adequately resolve gradients in the radiation energy density.

These conditions are applied recursively to every level up to some pre-specified maximum level L. In practice, the fourth condition is usually the most stringent.

For all our calculations we use symmetry boundary conditions on the hydrodynamics, but we use adaptivity to remove the boundary far enough from the cloud we are simulating so that no part of it ever approaches the boundary. The gravity module uses Dirichlet boundary conditions on the potential, with the potential on the boundary set equal to the value obtained from an octopole expansion of the density distribution inside the computational domain. The radiation module uses Marshak boundary conditions, meaning that radiation energy generated within the domain is free to escape. The incoming radiation flux is set to that appropriate for a blackbody radiation field at a temperature of 20 K.

The initial setup for all of our simulations is similar to those of KKM07. In all cases, the setup consists of an initially spherical cloud of mass M = 100 M_☉, initial mean surface density Σ, and radius $R = \sqrt{M/(\pi \Sigma)}$ in the center of a cubical computational domain. Observed dense star-forming clumps of molecular gas have roughly power-law density structures $\rho \propto r^{-k_{\rho }}$ with k_ρ ≃ 1.5 and considerable scatter (Caselli & Myers 1995; Beuther et al. 2002, 2005, 2006; Mueller et al. 2002; Sridharan et al. 2005), so, following McKee & Tan (2003), we adopt a power-law density structure with k_ρ = 1.5 for our initial conditions. We give our clouds an initial turbulent velocity field chosen to put them in approximate balance between gravity and turbulent ram pressure. The velocity dispersion is

$$\begin{equation} \sigma _v = \sqrt{\frac{GM}{2(k_{\rho }-1) R}}, \end{equation} \tag{ 15 }$$

and the initial velocity field we use is identical to that of run 100A of KKM07, which is a Gaussian-random field with a power spectral density P(k) ∝ k⁻². Finally, we give the gas in the cloud an initial temperature T_g = 20 K, and we set the initial radiation energy density throughout the computational domain to E = 1.2 × 10⁻⁹ erg cm⁻³, the value for a blackbody radiation field at a temperature T_r = 20 K.

Outside the cloud we place a hot, diffuse ambient medium with a density ρ_a = ρ_edge/100, where

$$\begin{equation} \rho _{\rm edge} = \left(\frac{3-k_{\rho }}{4\pi }\right) \frac{M}{R^3} \end{equation} \tag{ 16 }$$

is the density at the cloud edge. The ambient medium has a temperature T_a = 100T_g, ensuring that it is in thermal pressure balance with the cloud. To ensure that the ambient medium does not cool, radiatively heat the cloud, or interfere with radiation escaping from the cloud, we set the opacity of the ambient medium to a numerically small value.

We simulate three different clouds, chosen with values of Σ to be representative of three different types of star-forming environment as discussed in Section 1. The first run has Σ = 0.1 g cm⁻², typical of diffuse star-forming clouds such as Perseus and Ophiuchus. The second has Σ = 1.0 g cm⁻², typical of regions of massive star formation in the Galaxy. The third has Σ = 10 g cm⁻², an extremely high surface density found only in clusters near the Galactic center and in extragalactic super-star clusters. However, this type of star-forming environment is believed to have been more common earlier in cosmological evolution. We call the runs L, M, and H, for low, medium, and high column densities. We summarize the setup of the three runs in Table 1. We run each simulation for a time t = 0.6t_ff, where

$$\begin{equation} t_{\rm ff}= \sqrt{\frac{3\pi }{32G\overline{\rho }}} \end{equation} \tag{ 17 }$$

is the free-fall time computed at the mean density $\overline{\rho } = 3M/(4\pi R^3)$ of the cloud. By this point, the differences between the runs are clearly established, and the fraction of the collapsed mass in the most massive star asymptotes to a constant value (see below).

We emphasize that we have chosen the numerical setup so that the runs are, as much as possible, simply rescaled versions of one another. The initial conditions have identical density structures, virial ratios, and velocity fields, and the refinement criteria and peak resolution are the same in every run. We simulate each cloud for the same number of free-fall times. The homology between the runs is broken only by the influence of radiation. Radiation fixes the gas temperature (and thus causes the initial Mach numbers of the runs to vary slightly) and, much more importantly, the very different optical depths of the different clouds cause them to respond differently once stellar feedback begins. Thus, we expect any difference between the three runs to be almost entirely dictated by their differing response to stellar radiative forcing.

We show the large-scale evolution of runs L, M, and H in Figure 1. As the plot shows, and as expected, the three runs are essentially homologous on large scales. Once we account for the scaling of cloud radius and surface density between the runs, the only noticeable difference is a slight trend toward increasingly filamentary structures as Σ increases. This is easily understood as a radiative effect. As noted above, runs with higher Σ require larger velocity dispersions σ_v to maintain initial virial balance, while the sound speed is fixed by radiative effects. Thus, the Mach number of the initial turbulence increases with Σ, and this produces the observed increase in filamentary structure.

Zoom In Zoom Out Reset image size

In Figure 2, we show the simulations at the same times as in Figure 1, but now zoomed in to a small region centered on the most massive star present, or on the origin before stars form. In the left panel of Figure 2, the images are scaled homologously, so that the regions shown all have a length equal to 10% of the initial cloud radius, and the color scale is in units of surface density divided by initial mean surface density. In the right panel, the scaling is physical, so that all the plots show a region of fixed physical size, using a column density scale in fixed rather than normalized units.

Zoom In Zoom Out Reset image size

In contrast to the large-scale homology seen in Figure 1, on small scales the runs rapidly diverge. Deviations from homology start to appear around 0.2t_ff, and by 0.4t_ff any visual similarity between the runs is completely gone. Progressing from low to high Σ, runs are characterized by increasing surface densities even when normalized to the initial surface density. In run H, the predominant structure is a single large disk concentrated around a single central object. In run M, we have a massive binary with two circumstellar disks and a larger circumbinary disk. Finally, in run L the disks are much smaller and less dense, and they have mostly depleted by the final time shown.

The higher Σ runs also fragment less, producing fewer, more massive stars than the low Σ runs. We show this in Table 2 and Figure 3. The total mass in stars M_*,tot at any given time (normalized to t_ff) is very similar from run to run, varying by less than 10% from run L to run H at times >0.2t_ff. At the end of the simulation, each run has a total of 15 M_☉ of stars. This reflects that the star formation rate in the simulations is driven by large-scale flows that are changed very little by radiation feedback. However, the pattern of fragmentation is quite different. The mass of the most massive object M_*,max in the three runs begins to diverge at 0.2–0.3t_ff, and thereafter it increases much more rapidly in run H than in run L. At the final time, the difference in mass is nearly a factor of 3, with run L having a maximum stellar mass below 6 M_☉ and run H reaching 15 M_☉. In all the runs, the fraction of the total stellar mass in the most massive object f_max asymptotes to a roughly constant value after 0.3t_ff. The asymptotic value ranges from f_max ∼ 0.35 in run L to f_max ∼ 0.9 in run H. In effect, the initial gas cloud in run H is like a single massive protostellar core that forms one massive star plus a few small secondaries, while the cloud in run L instead forms a small cluster that does not include any massive stars. Run M is intermediate.

Zoom In Zoom Out Reset image size

The difference in runs is even more apparent if we focus on a single time. Figure 4 shows the cumulative distribution function of stellar mass at t = 0.6t_ff in each of the runs. In run L, the distribution function rises relatively smoothly between 1 and 5 M_☉, so the system consists of a number of stars of roughly comparable mass. In contrast, in run M 90% of the mass is in just two stars that form a binary system, a result very similar to that in KKM07, which used an initial surface density Σ = 0.7 g cm⁻². In run H, a comparable fraction of the mass is in a single star. Since the system in run L involves a number of stars of comparable mass, these are unlikely to wind up as a bound star system. The system in run M, on the other hand, is stable and is likely to remain a bound binary. Thus in both runs M and H, the result is that most of the mass goes into a single star system, while in run L the mass will end up divided into several star systems.

Zoom In Zoom Out Reset image size

The difference in morphology, fragmentation, and star formation is easy to understand if we examine the temperature structure of the gas. In Figure 5, we show the column-density-weighted temperature over the same small-scale regions as in Figure 2. Clearly, the temperature distribution in the gas in the runs is even less homologous than the density structure. Moreover, the differences in temperature begin to appear at earlier times. At t = 0.1t_ff, the density plots are nearly indistinguishable once the homologous scaling is removed (left panel of Figure 2), while the differences in temperature are already obvious. It is important to point out that at this point there are no stars present that are producing significant amounts of power via nuclear luminosity. The most massive star present in any of the runs is only 0.37 M_☉, and its luminosity is entirely driven by accretion. The difference in temperature is therefore solely due to the two factors pointed out by KM08: compared to run L, run H has both a higher density to produce higher accretion rates and thus higher accretion luminosities, and a higher optical depth to more effectively trap the radiation that is produced.

Zoom In Zoom Out Reset image size

We can understand the difference in the thermal structure of the runs more quantitatively and explore how this difference is likely to influence fragmentation by examining the relationship between temperature and density in each of the runs. Figure 6 shows the locus occupied by the clouds in runs L, M, and H at different times in the simulation in the plane of log density versus log temperature. The color represents the mass density at a given point, and the contours indicate, from lowest to highest, regions in the plane containing 99.9%, 99%, 90%, and 50% of the gas mass in the computational domain.

Zoom In Zoom Out Reset image size

The figure demonstrates how different the thermal structure of the gas is in each of the runs. In run L, the great majority of the gas is near the background temperature of 20 K at all times. Even at the final time only 10% of the gas mass is at noticeably elevated temperatures, and less than 1% of the mass is heated above 100 K. In contrast, in run H the heating is much more extensive. Even at t = 0.1t_ff, when the only source of heating is the accretion luminosity of a 0.37 M_☉ star, the contours containing 50% and 90% of the mass are noticeably elevated above the T = 20 K line. Deuterium burning in the star begins shortly before 0.2t_ff, and by the final time deuterium burning and Kelvin–Helmholtz contraction (the star has not yet reached the main sequence) provide enough luminosity to keep all the mass at temperatures ≳50 K. Run M is intermediate, with small but significant fractions of the cloud mass reaching elevated temperatures at early times, and more mass becoming heated as the run progresses and the stars become more luminous. This is consistent with the analytic predictions of KM08, who find that a column density of 1 g cm⁻² constitutes the rough line between clouds that do and do not experience significantly elevated temperatures over much of their mass as a result of trapped accretion luminosity.

It is important to point out that the temperature does not need to rise to the point where the Jeans mass is above 100 M_☉ in order to inhibit fragmentation. Indeed, such a rise in temperature would be sufficient to halt collapse of the core entirely. Instead, the heating prevents fragmentation by creating an environment where the effective equation of state with γ = 1 + dlog T/dlog ρ>1 throughout the bulk of the cloud mass. Examining Figure 6, we see that the region containing 90% of the cloud mass (the third contour from the outermost one) is almost perfectly horizontal in run L at all times, so γ ≈ 1. In runs M and H, on the other hand, this region has a slope ∼0.2–0.3 in the log ρ–log T plane at all times of 0.2t_ff or more, similar to the result obtained by KKM07, indicating that the effective equation of state is closer to γ = 1.2–1.3. As Larson (2005) points out, and the simulations of Jappsen et al. (2005) confirm, fragmentation is likely as long as the effective equation of state for the gas is γ ⩽ 1, and is unlikely when γ>1. In runs M and H, there is no significant gas mass with γ ⩽ 1, which is why fragmentation is suppressed.

Finally, we emphasize that these phenomena cannot be correctly captured by analytic equations of state that are based on either a barotropic or and optically thin cooling assumption. Examples of such equations of state from Dobbs et al. (2005) and Larson (2005) are shown in Figure 6, and they clearly do not even come close to reproducing the results with radiative transfer, a point also made by Boss et al. (2000), Krumholz (2006), Krumholz et al. (2007a), and Offner et al. (2009b). Any such approximation would give the same temperature–density relation for all three of our simulated clouds, while clearly the results are different at different times and for different initial cloud column densities. Urban et al. (2010) reach the same conclusion for lower-density, larger-scale clouds based on their simulations.

Although we have not tested the Bate (2009) approach of omitting radiation from stars and including only radiative emission by gas on size scales resolved by the computation (which are much larger than stellar scales in both Bate's calculation and ours), it seems unlikely that this approximation could succeed either in the case of clouds with differing initial column densities. It would capture the difference in optical depth between runs, but it would not capture the effect that higher-density runs produce higher accretion rates and thus higher accretion luminosities from the embedded protostars. The analytic models of KM08 suggest that both effects are of comparable importance.

Our results demonstrate that the amount of fragmentation that a cloud undergoes is likely to depend strongly on its surface density, and this has important implications for where we expect massive stars to form. Consider a protostellar core, an object with a mass of a few tenths to a few hundreds of M_☉ that collapses to make one or more stars. Low-mass cores do not have significant internal turbulence (André et al. 2007; Kirk et al. 2007; Rosolowsky et al. 2008), as is expected on theoretical grounds (Offner et al. 2008b). Consequently, while they may fragment into a binary like run M, we expect most of the stellar mass they produce to end up in a single star system. The overall efficiency of turning gas mass into stellar mass is expected to be ≈ 1/3 rather than = 1 as a result of mass ejection by protostellar outflows (Matzner & McKee 2000; Alves et al. 2007; Enoch et al. 2008).

In contrast, as pointed out by McKee & Tan (2003), massive cores such as those we simulate are turbulent, since their masses are many times the thermal Jeans mass. We find that, at low Σ, such cores will fragment so that their stellar mass is divided among many star systems. Figure 7 gives a more detailed picture of how this happens in run L. There are five fragments whose masses appear to be asymptotically approaching fixed fractions of the total stellar mass, ranging from 11% to 37%, plus three more smaller stars whose masses appear to have reached nearly fixed maxima, and that are therefore declining with time in the total fraction of stellar mass that they represent.

Zoom In Zoom Out Reset image size

We can use this result to make a toy model for how the fraction of the stellar mass that is in massive stars is likely to vary with surface density. We begin from the observation that stars form from cores that have a mass distribution with the same functional form as the IMF, so that the IMF is set at the phase when gas fragments into protostellar cores (Motte et al. 1998; Testi & Sargent 1998; Johnstone et al. 2001; Onishi et al. 2002; Beuther et al. 2004; Reid & Wilson 2005, 2006a, 2006b; Alves et al. 2007; Nutter & Ward-Thompson 2007; Simpson et al. 2008; Enoch et al. 2008; Rathborne et al. 2009). We model this core mass function (CMF) using a Chabrier (2005) stellar system IMF shifted to higher mass by a factor of 3 to account for the mass that is ejected by outflows:

$$\begin{equation} \frac{dn_c}{d\ln m_c} = \left\lbrace \begin{array}{@{}ll} A \exp \left[-\frac{(\ln m_c - \ln \overline{m_c})^2}{2\sigma ^2}\right], \quad & m_c < m_{\rm break}, \\ B \left(\frac{m_c}{m_{\rm break}}\right)^{-1.3}, & m_c \ge m_{\rm break}, \end{array} \right. \end{equation} \tag{ 18 }$$

with $\overline{m_c} = 0.75 \,M_{\odot }$, m_break = 3 M_☉, and σ = 0.55. Our values of $\overline{m_c}$ and m_break are chosen so that, when = 1/3, the stellar IMF will have a peak at 0.25 M_☉ and will break from lognormal to power law form at 1.0 M_☉, in agreement with Chabrier's best fit to observations. The normalization factors are related by $B = A \exp [-(\ln m_{\rm break} - \ln \overline{m_c})^2/(2\sigma ^2)]$. Theoretical models are able to explain this distribution of core masses as arising naturally from the properties of supersonic turbulence (Padoan & Nordlund 2002; Hennebelle & Chabrier 2008). In this picture, the final stellar IMF is simply the convolution of the CMF with a simple function that maps core mass to stellar mass and which must be nearly mass independent. Clark et al. (2007) suggest that this correspondence will be disrupted if the core free-fall time is mass dependent. However, observed cores do not have mass-dependent free-fall times (André et al. 2007), and theoretical models predict that, contrary to Clark et al.'s assumption, core free-fall time depends on mass at most very weakly (McKee & Tan 2003; Hennebelle & Chabrier 2009) as a result of turbulent support.

To make a quantitative model of how fragmentation in low Σ regions will affect the IMF, we consider two extreme scenarios that bracket the outcome of run L. The first is that, in regions of low Σ, cores with mass m_c above some minimum fragmentation mass m_frag produce only a single massive star with a mass m_* = m_c/9, i.e., 2/3 of the core mass is ejected by outflows, 1/3 of what remains goes into the largest star, and the remaining 2/3 goes into low-mass stars with m_* < m_frag/3. The alternative possibility is that the cores with mass m_c > m_frag fragment into n_frag stars of equal mass m_* = m_c/(3n_frag), with the factor of 3 to account for ejection by outflows. Based on the outcome in run L, we adopt n_frag = 5 for our toy model, although of course in reality we expect that the number of fragments and their mass distribution will vary stochastically.

In the first scenario, if the total mass of cores with masses between m_c and m_c + dm_c is given by dn_c/dln m_c, the corresponding mass of stars with masses between m_* = m_c/3 and m_* + dm_* is given by dn_*/dln m_* = (1/9)dn_c/dln m_c for any mass m_* > m_frag. In the second scenario, we instead have m_* = m_c/(3n_frag) and dn_*/dln m_* = (1/3)dn_c/dln m_c, since all of the core mass that is not ejected by outflows (i.e., 1/3 of it) goes into stars of mass m_c/(3n_frag). Finally, in high Σ regions where massive cores do not fragment, we also have dn_*/dln m_* = (1/3)dn_c/dln m_c, but now the stellar mass is related to the core mass by m_* = m_c/3 rather than m_* = m_c/(3n_frag).

Given these relations, we can compute the fraction of all stellar mass that is contained in stars with masses greater than m_*, which we denote f(>m_*). We adopt a minimum stellar mass m_*,min = 0.01 M_☉ and a maximum m_*,max = 120 M_☉ (Figer 2005). First consider regions of high Σ where massive cores do not fragment and stellar and core masses are related by m_* = m_c/3. In such a region, we have

$$\begin{equation} f({>}m_*) = \frac{ \int _{3 m_*}^{3 m_{*,\rm max}} \left(\frac{1}{3}\right)\frac{dn_c}{d\ln m_c} \, dm_c }{ \int _{3 m_{*,\rm min}}^{3 m_{*,\rm max}} \left(\frac{1}{3}\right)\frac{dn_c}{d\ln m_c} \, dm_c } \equiv f_H({>}m_*). \end{equation} \tag{ 19 }$$

This equation simply states that the fraction of stellar mass in stars larger than m_* is the same as the fraction of core mass in cores larger than 3m_*. For low Σ regions, in the first scenario we instead have

$$\begin{equation} f({>}m_*) = \frac{ \int _{9 m_*}^{9 m_{*,\rm max}} \left(\frac{1}{9}\right)\frac{dn_c}{d\ln m_c} \, dm_c }{ \int _{3 m_{*,\rm min}}^{9 m_{*,\rm max}} \left(\frac{1}{3}\right)\frac{dn_c}{d\ln m_c} \, dm_c } \equiv f_{L,1}({>}m_*) \end{equation} \tag{ 20 }$$

for stellar masses m_* > m_frag/3. Note the factors of 1/9 in the numerator and 1/3 in the denominator, reflecting that, in this scenario, only 1/9 of the mass in a core of mass m_c is incorporated into a star of mass m_* = m_c/9, but that 1/3 of the core mass goes into stars overall. If we instead adopt the second scenario, where 1/3 of the core mass goes into n_frag stars of mass m_* = m_c/(3n_frag), we obtain

$$\begin{equation} f({>}m_*) = \frac{ \int _{3 n_{\rm frag} m_*}^{3 n_{\rm frag} m_{*,\rm max}} \left(\frac{1}{3}\right)\frac{dn_c}{d\ln m_c} \, dm_c }{ \int _{3 m_{*,\rm min}}^{3 n_{\rm frag} m_{*,\rm max}} \left(\frac{1}{3}\right)\frac{dn_c}{d\ln m_c} \, dm_c } \equiv f_{L,2}({>}m_*). \end{equation} \tag{ 21 }$$

A final complication is that, in both scenarios, we have implicitly assumed that the CMF extends to infinity, or at least to 3n_fragm_*,max = 1800 M_☉. This is not necessarily the case—the origin of the observed cutoff in the IMF is not understood, and one possible explanation for it is that there simply are no protostellar cores whose mass is larger than a few hundred M_☉ that are capable of collapsing to single stars even in regions of high surface density. As a simple example of how this would change the results, suppose that there is a maximum core mass m_c,max = 3m_*,max = 360 M_☉, sufficient to make a 120 M_☉ star if the core does not fragment. In this case, the mass fractions f_L,1 and f_L,2 we have just computed are modified to

$$\begin{eqnarray} f_{L,1a}({>}m_*) & = & \frac{ \int _{9 m_*}^{m_{c,\rm max}} \left(\frac{1}{9}\right)\frac{dn_c}{d\ln m_c} \, dm_c }{ \int _{3 m_{*,\rm min}}^{m_{c,\rm max}} \left(\frac{1}{3}\right)\frac{dn_c}{d\ln m_c} \, dm_c }, \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} f_{L,2a}({>}m_*) & = & \frac{ \int _{3 n_{\rm frag} m_*}^{m_{c,\rm max}} \left(\frac{1}{3}\right)\frac{dn_c}{d\ln m_c} \, dm_c }{ \int _{3 m_{*,\rm min}}^{m_{c,\rm max}} \left(\frac{1}{3}\right)\frac{dn_c}{d\ln m_c} \, dm_c }. \end{eqnarray} \tag{ 23 }$$

Evaluating the functions f_H, f_L,1, f_L,2, f_L,1a, and f_L,2a gives the results shown in Figure 8. As the plot shows, reduced star formation efficiency due to the fragmentation in massive cores of the sort we have found can reduce the fraction of total stellar mass in high-mass stars by a significant amount. The minimum reduction in massive star fraction, by 40%, is for f_2,L, corresponding to the scenario where massive cores fragment into a few equal mass objects and the CMF extends to infinity. In this case, the full mass of cores that fragment is still available to make massive stars, and the massive star fraction declines only because stars of mass m_* must be produced by cores of mass 3n_fragm_* in regions of low Σ rather than by cores of mass 3m_* in regions of high Σ, and the total mass of cores available is lower for higher m_c.

Zoom In Zoom Out Reset image size

Any other scenario gives a much more significant reduction in the massive star fraction in regions of low surface density. In scenario f_L,1, where fragmenting cores produce one massive star of mass m_c/9 and only low-mass stars otherwise, the reduction is by 76%. The reduction is partly for the same reason as for f_L,2, and partly because all the mass in the massive core that does not go into massive stars still goes into low-mass stars. Finally, if the CMF has a cutoff, the reduction in massive star fraction is even more dramatic. In this scenario, fragmentation means that a core of mass m_c = 9m_* (for f_L,1a) or of mass m_c = 3n_fragm_* (for f_L,2a) is required to make a star of mass m_*, and if this exceeds the maximum core mass then no stars of mass m_* can form in regions of low Σ.

Regardless of which scenario is ultimately correct, our results show that in regions of low surface density we expect a significant decline in the stellar mass fraction in massive stars compared to a canonical IMF. The reduction is anywhere from a factor of 1.7 in the most conservative scenario to a very large factor in more liberal scenarios. This provides numerical confirmation of the hypothesis advanced by KM08 that there is an effective threshold for massive star formation.

A variable IMF has numerous implications on scales ranging from the sub-galactic to the cosmological, and a full exploration of them is beyond the scope of this paper. Moreover, a full exploration of this topic would require a larger parameter study than the one we present here, allowing a full exploration of how fragmentation and the stellar mass function vary with environment. Nonetheless, we can identify some important implications.

On the scales of star clusters, preferential formation of massive stars in regions of high surface density suggests that clusters are likely born mass segregated, with more massive stars forming in the center where the surface density is highest. Outer parts of the cluster, where the surface density drops below ∼1 g cm⁻², should form preferentially low-mass stars. There is some evidence for such primordial mass segregation in Orion (Hillenbrand & Hartmann 1998; Huff & Stahler 2006), but debate continues about whether mass segregation in clusters in general is a result of formation (Bonnell & Davies 1998) or dynamical processes during the first few Myr of cluster lifetime (Tan et al. 2006; McMillan et al. 2007). Both may occur together, though recent observational (Fűrész et al. 2008; Tobin et al. 2009) and theoretical (Offner et al. 2008b, 2009a) work pointing out that stars are born sub-virial with respect to their parent clouds (even if the clouds themselves are virialized) suggests that dynamical segregation may be much stronger than earlier estimates found (Allison et al. 2009). In very dense clusters, primordial mass segregation may also strongly influence the dynamical evolution of the cluster (e.g., Chatterjee et al. 2009).

Conversely, we do not expect a correlation between the IMF and cluster mass beyond the trivial one expected from finite sampling of a universal IMF, since there does not appear to be a strong correlation between protocluster gas cloud masses and surface densities (Fall et al. 2010). A correlation between cluster mass and the IMF has been suggested by some models (e.g., Weidner & Kroupa 2006), based on the idea that massive stars form via a process in which all gas clouds fragment down to the small initial Jeans mass, but that some subsequently grow through competitive Bondi–Hoyle accretion (Bonnell et al. 2004). However, the premise on which these models are based—fragmentation of all clouds down to the Jeans mass at the initial, low cloud temperature—clearly fails when radiation is included. While the process that goes on in run L might roughly be described as competition, there is clearly no competition in run M or H, where radiation ensures that most of the proposed competitors never form in the first place.

Observations remain divided on whether there is a correlation between cluster mass and maximum stellar mass. Weidner & Kroupa (2004, 2006) and Weidner et al. (2010) claim to detect one, while Oey et al. (2004), Elmegreen (2006), and Parker & Goodwin (2007) argue that there no correlation is present. de Wit et al. (2004, 2005) find that 4% ± 2% of galactic O stars formed outside of a cluster of significant mass, which is consistent with the models presented here (for example, runs M and H form effectively isolated massive single stars or binaries), but not with the proposed cluster–stellar mass correlation.

No comparable studies have been conducted to search for systematic variation of the IMF with surface density, which we do predict. Such studies are likely to be even more challenging than those searching for a correlation with cluster mass, because cluster surface densities evolve very rapidly once star formation ends and the remaining gas is expelled. For this reason, any search for a correlation between IMF and surface density would have to target clusters that are still embedded in their parent gas clouds, for which the cloud surface density can be measured. Unfortunately, this renders optical observations, the most common method for determining stellar masses, impossible, since a surface density Σ = 1 g cm⁻² corresponds to A_V ≈ 200. Instead, other methods of estimating populations of low- and high-mass stars, such as X-ray observations, would be required (Krumholz & McKee 2008).

On larger scales, our ability to make predictions is limited by our lack of a theoretical model capable of connecting the surface densities measured over the small scales of star-forming regions where fragmentation occurs (∼1 pc) to those averaged over much larger areas of galactic disks (∼1 kpc). Regions of high surface density where massive stars form are clearly found preferentially in regions of high galactic surface density such as spiral arms and near galactic centers, while low surface density clouds such as Taurus are found in regions of lower surface density. However, there is no clear one-to-one mapping from large to small scales. Nonetheless, it seems clear that our results do predict a suppression in the formation of high-mass stars in regions where the galactic surface density is low, such as dwarf galaxies and the outer parts of spiral galactic disks. Such a correlation may have been observed (Boissier et al. 2007; Meurer et al. 2009; Lee et al. 2009), although it too remains controversial (Boselli et al. 2009).