Explaining the stellar initial mass function with the theory of spatial networks
Abstract
The distributions of stars and prestellar cores by mass (initial and dense core mass functions, IMF/DCMF) are among the key factors regulating star formation and are the subject of detailed theoretical and observational studies. Results from numerical simulations of star formation qualitatively resemble an observed mass function, a scalefree power law with a sharp decline at low masses. However, most analytic IMF theories critically depend on the empirically chosen input spectrum of mass fluctuations which evolve into dense cores and, subsequently, stars, and on the scaling relation between the amplitude and mass of a fluctuation. Here we propose a new approach exploiting the techniques from the field of network science. We represent a system of dense cores accreting gas from the surrounding diffuse interstellar medium (ISM) as a spatial network growing by preferential attachment and assume that the ISM density has a selfsimilar fractal distribution following the Kolmogorov turbulence theory. We effectively combine gravoturbulent and competitive accretion approaches and predict the accretion rate to be proportional to the dense core mass: . Then we describe the dense core growth and demonstrate that the powerlaw core mass function emerges independently of the initial distribution of density fluctuations by mass. Our model yields a power law solely defined by the fractal dimensionalities of the ISM and accreting gas. With a proper choice of the lowmass cutoff, it reproduces observations over three decades in mass. We also rule out a lowmass star dominated “bottomheavy” IMF in a single starforming region.
Subject headings:
stars: formation, stars: luminosity function, mass function, ISM: clouds, ISM: structure1. Introduction and Motivation
Six decades ago the stellar initial mass function (IMF) was derived from star counts (Salpeter, 1955) as a scalefree power law (; ) with more frequent lowmass stars than highmass stars. Since then, it has attracted attention as one of the principal star formation characteristics that controls stellar feedback and, therefore, governs galaxy evolution. Furthermore, explaining the IMF will help us to understand the star formation physics. From observations, the IMF shape appears to be universal across different starforming regions (Kroupa, 2002). Resembling a unimodal (Salpeter, 1955) or bimodal (Kroupa, 2001) power law or a lognormal distribution with a powerlaw tail (Chabrier, 2003), it sharply declines at the low end at masses solar mass (). The dense core mass function (DCMF) derived from observations of giant molecular clouds (Alves et al., 2007; André et al., 2010) is an IMF precursor: First, dense cores grow from density fluctuations by attracting surrounding material, cool down, then protostars form inside them and evolve into stars. The DCMF shape also looks like a power law at high masses and declines below 1/3 , offset by a factor of to higher masses compared to the stellar IMF, illustrating a 25% gastostars transformation efficiency. Thus, if robust arguments were provided to explain the DCMF shape, the IMF shape would follow through that heuristic conversion rule.
All existing analytic and numerical IMF theories (see the review by Hennebelle & Chabrier, 2011) consider either the accretion of material on protostars (Zinnecker, 1982; Bonnell & Bate, 2006) or the gravitational fragmentation of the interstellar medium (ISM) (Padoan & Nordlund, 2002). In a simple model, a nonlinear stage of the molecular cloud fragmentation yields a lowmass IMF decline (Silk & Takahashi, 1979) but does not reproduce a powerlaw highmass tail. Diverse physical mechanisms of the molecular cloud cooling that affect the collapse and fragmentation are often hidden in a complex equation of state for molecular clouds where the polytropic index depends on density, temperature, and chemical composition (Spaans & Silk, 2000). One of the most complete analytic IMF theories up to date (Hennebelle & Chabrier, 2008, 2013) explains the observed overall IMF shape by analyzing the evolution of density fluctuations in a selfgravitating turbulent ISM in the presence of a magnetic field. Similarly to other gravoturbulent theories, in order to reproduce a powerlaw part of the mass function consistent with the Salpeter slope, it relies on a specific choice of the scaling relation that connects the density fluctuation amplitude with the mass contained within that fluctuation, and a lognormal initial ISM density probability density function (PDF) required for the analytic computation of the mass function shape. However, the density PDF scale dependence chosen in Hennebelle & Chabrier (2008) relies on results from numerical simulations. Also, density PDFs observed in molecular clouds deviate substantially from the lognormal shape and vary across different starforming regions (Lombardi et al., 2015). Moreover, gravoturbulent theories do not consider any external accretion on dense core progenitors, so there is no guarantee that a powerlaw DCMF holds as the system evolves.
Recent developments by Hopkins (2013) introduce a variety of modifications to the gravoturbulent fragmentation conditions and propose different grounds for the density PDF scaling relation. The author claims that the input ISM density distribution does not have to be lognormal. However, this might be the result of the mathematical simplification applied (the tophat filtering in the Fourier space and the consequent Taylor expansion) that transforms an arbitrary function into a lognormal like shape. Nevertheless, this theory yields a fluctuation mass function consistent with the Salpeter slope after some fine tuning of the model parameters. At the same time, no argument is provided for the fundamental reasons of the power law and a particular value of the exponent.
The generalization of gravoturbulent theories to nonlognormal or nonanalytic PDF shapes becomes cumbersome and not very straightforward. Schmidt et al. (2010) demonstrate analytically and numerically that gravoturbulent theories can reproduce the powerlaw IMF for arbitrary initial PDF shapes, however, in their analytic computation they rely on specific simplifications (e.g. excluding the largest cores from consideration in the case of the Hennebelle–Chabrier theory), which de facto restricts the input density PDFs to certain functional families.
Whereas analytic theories that deal with the competitive accretion scenario reproduce the powerlaw IMF tail, they fail to match its observed exponent and cannot be applied to nonclustered star formation. A recent IMF theory from that family (Basu et al., 2015) generates a lognormal distribution with a powerlaw tail via the process of quenched accretion with the exponential distribution of accretion timescales. However, that theory neither provides a justification for the exponential distribution nor does it proposes a quantitative argument for the slope of the powerlaw tail.
The theory by Maschberger (2013) considers both linear and nonlinear accretion of mass onto dense cores and takes into account the stochasticity of the process using the Stratonovich stochastic calculus. While it successfully reproduces highmass powerlaw tails, the accretion laws considered there are not connected to the ISM properties.
Hence, the following questions still remain unanswered by existing IMF theories: (i) an a scalefree powerlaw distribution of dense core masses become established by some physical processes independent of the initial density distribution; (ii) does it hold as the system evolves by accreting external material; and (iii) what does the powerlaw exponent depend on?
Here we present a different approach to the analytic DCMF theory. We describe the mass accretion onto prestellar dense cores in the fractal ISM as preferential attachment, a key phenomenon studied in the field of network science (Merton, 1968; Newman, 2010; Barthélemy, 2011). We use probabilistic accounting of small parcels that join dense cores, subject to gravitational attraction and stochastic noise. We limit our model to the early stages of the dense core growth by mass accretion without complex physics following the protostar formation. Therefore, it applies to starless dense cores. The principal result of our theory is an analytic expression for the powerlaw tail that develops for any initial distribution of dense core progenitors by mass. This is the first example of a theory that consistently connects the accretion rate to global properties of the turbulent ISM and, thus, effectively combines the gravoturbulent and competitive accretion approaches.
2. Mass distribution and a system of units
We describe a system of dense cores as an array of masses distributed according to some timedependent function . At a given moment, the system has a total mass and a total number of dense cores . The mass distribution shape is governed by the two processes: accretion of mass parcels onto existing cores and generation of new core progenitors. We first treat events of parcel accretion and, thus, masses of dense cores as discrete and further take the continuum limit.
Let us consider an event where a gas parcel of a small mass per one unit of time attaches to one of the cores. Regardless of which core it joins, the total mass of the system grows by , hence . The accretion summed over the entire system (e.g. the global mass growth of the system) happens at some characteristic rate . Over a time step every mass bin can only be affected by bins within from it. Hence, the timescale choice uniquely defines the mass bin size. A volume density of new dense core progenitors created per unit time in each bin is described by some function . This function is the only mechanism that increases , therefore . We stress that the and values are, in general, timedependent. As the accretion exhausts available material in the surrounding ISM, both rates should slowly decay to zero. However, as we show in Section 5, neither their exact time dependence nor their absolute values matter for the final DCMF shape.
Then, the mass distribution of cores is expressed as with , discrete time , and normalization . We are interested in its long term behavior, such that .
3. The fractal matter distribution in the ISM
Observations suggest that density and velocity distributions in the ISM are predominantly defined by turbulent motions on scales from hundredths of a parsec to hundreds of parsecs (Elmegreen & Scalo, 2004). A consequence of the Kolmogorov theory (Kolmogorov, 1941) is a selfsimilar or fractal density distribution in a turbulent flow (Sreenivasan et al., 1989) with the predicted fractal dimensionality . It stays in agreement with measurements obtained from observations of giant molecular clouds (Falgarone et al., 1991; Elmegreen & Falgarone, 1996) () and in laboratory studies of turbulence (Sreenivasan et al., 1989) (). Numerical simulations of the ISM evolution with an input fractal density field (Elmegreen, 1997, 2002) yield powerlaw mass distributions of overdensities that correspond to dense cores in starforming regions.
We describe ISM as a twophase medium where the two phases may have different fractal dimensions. One phase corresponds to dense cores. The other one corresponds to parcels, small gas/dust fragments with individual masses substantially below the turbulent and gravitational Jeans masses (Hennebelle & Chabrier, 2008; Hennebelle & Falgarone, 2012) that do not have to obey the same law. Dense cores arise from the initial turbulent medium and, therefore, their positions trace the initial overdensities in the turbulent flow. Parcels correspond to all remaining material of the ISM that did not enter dense core overdensities initially, but can be accreted by them. The dense core phase defines the gravitational field profile in the ISM, and the parcel phase moves in that field and accretes on dense cores. The spatial distributions of the two phases are governed by different mechanisms and thus are neither positively nor negatively correlated.
Our mathematical description of the twophase ISM follows a fractal model by Tarasov (2005). Both phases have a characteristic “microscopic” lengthscale (pore size in Tarasov (2005)) on which the discreteness of the medium is visible. For the dense core phase this characteristic lengthscale corresponds to the average distance between adjacent dense cores. It can be inferred from the Jeans mass and the total system mass. For the parcel phase, the lengthscale corresponds to the average separation between particles of gas or dust. Presumably , which corresponds to dense cores being much larger and less frequent than parcels. However, for the purpose of this calculation we are interested in the statistics on much larger scales.
As Tarasov (2005) suggests, for a fractal medium observed on scales or , uniform fractal scaling is observed, i.e. the mass confined in a sphere of radius grows as , where is the fractal dimensionality of the density distribution. This relation has two important properties, fractality and homogeneity. For the uniform distribution of matter in 3D space , but for the fractal distribution . Consequently, the average density in a sphere is not constant, but is scaledependent. The homogeneity means that the fractal powerlaw scaling of mass confined in a sphere is independent of the position of that sphere. Features and structures of the fractal distribution are associated with observing it on different scales rather than at different positions. Selfsimilar (or fractal) scaling appears only on scales when discrete features of characteristic scales and blur out.
Since the dense core distribution traces the initial Kolmogorovlike supersonic turbulent distribution, we take it to have . Confining our system to a box of linear size and normalizing its mass to , the mass confined in a thin concentric sphere becomes . The fractal dimension of parcels is left as a free parameter for now and is discussed below in more detail. Therefore, the average number of parcels confined within a thin concentric sphere, properly normalized, is .
Even if a substantial mass fraction is contained in the diffuse phase (e.g. an order of 50%), our calculations will remain valid, because the gravitational field gradients and, correspondingly, the basins of attraction (see Fig. 1 right) will still be defined by “pointlike” dense cores.
4. Dense core growth by preferential attachment
Preferential attachment is a stochastic process in which a set of objects possessing some property acquire discrete units of this property in a partly random fashion such that the probability of a given unit to be attached to a given object increases with the increase of the amount of that property already contained in this object. It is also referred as a “Yule process” in speciation (Yule, 1925), a “Matthew effect” in science organizations (Merton, 1968), a “cumulative advantage” in bibliometrics (Simon, 1955; Price, 1976), and as a “capital gain” in economics (Yakovenko & Rosser, 2009). Preferential attachment in random networks naturally explains powerlaw distributions (Barabási & Albert, 1999) of node sizes defined by the number of links. This approach explained power laws emerging across different fields of science, e.g. in the World Wide Web structure (Albert et al., 1999), protein interactions (Jeong et al., 2001), metabolics (Ravasz et al., 2002), transportation, and social networks, and scientific collaborations (Barabási et al., 2002; Newman, 2004). Here we describe a system of dense cores growing in a molecular cloud by preferential attachment. Gravitational forces representing “links” between dense cores are distancedependent, hence we exploit the spatial network formalism (Barthélemy, 2011).
When a new parcel emerges in the system, it becomes subject to multiple competing attractive gravitational forces from existing dense core progenitors, and at the same time to drag forces as it moves through the ISM. We assume that drag forces dominate over the inertia, so that the exact dense core which will acquire a given parcel is determined only by the competition of forces at the parcel’s starting position (see Fig. 1, right panel). We set the parcel accretion probability by a given core proportional to the initial gravitation acceleration toward it (Fig. 2, left panel). In close vicinities of dense cores, where the gravitational field is totally dominated by one mass, our description becomes equivalent to the deterministic accretion onto that particular core. However, at the border separating areas of dominant attraction (Fig. 1) from two cores, a parcel can be tipped over it by stochastic pushes from other particles of the ISM. The probabilistic approach allows us to model that situation. Hence, the probability of a newly emergent parcel to join an existing dense core is:
(1) 
where is the mass of the th core, is the Euclidean distance between the points, is the gravitational law exponent ( in the 3dimensional space). A possible alternative choice is which corresponds to the decay of gravitational potential rather than acceleration. This, however, is not very important for the preferential attachment description. We choose a normalization so that the probability of a parcel joining some dense core is unity by using the continuous random fractal approximation to sum over all cores:
(2) 
where is a distance from a new parcel to the nearest dense core. Assuming and , we can neglect so that the statistics is dominated by distant dense cores, owing to the long distance nature of gravitation. Thus, the normalization factor is the same for all parcels, regardless of where they emerge.
From examining the integrals in Eqs. 2 and 3 below we can see that distant spherical layers of the ISM contribute very little to the dense core growth or the parcel accretion, since their contribution is proportional to or and vanishes in the limit of large .
The average dense core mass increase per time step is a probability weighted sum over all possible positions where parcels emerge.
(3) 
We notice the difference between our calculated accretion rate for an individual core and the classical Bondi–Hoyle–Lyttleton model (Hoyle & Lyttleton, 1939; Bondi, 1952) of spherical accretion (). This inconsistency is trivially explained by the two facts: (a) we assume a fractal distribution of the infalling matter so that a thin spherical layer no longer contains the mass used in the Hoyle–Lyttleton calculations; and (b) our model considers motions of fractally distributed gas parcels in space to be overdamped as opposed to ballistic motions, therefore the whole orbital computation including the impact parameter and the escape velocity from Hoyle & Lyttleton (1939) is not applicable to our case.
We denote the growth exponent . It characterizes the growth rate of individual masses. We can illustrate its contribution to the dense core growth by using a simple, but manifestly unrealistic assumption of a constant accretion rate in the system. If at some moment , a dense core has mass and the global mass growth rate is quasiconstant ( and ), then it grows in time according to a sublinear power law:
(4) 
In reality, this law does not hold because the global mass growth rate might not be constant. However, one cannot directly observe the growth of a single dense core because it lasts tens of thousands of years. As we show in the next section, the directly observable quantity is a snapshot of the DCMF in the limit, which in turn is not affected by the specific time dependence of and .
5. The powerlaw distribution from the master equation in networks
The growth law for an individual dense core is not sufficient to derive the mass distribution shape. Therefore we use the master equation (Dorogovtsev & Mendes, 2002; Newman, 2010) for the distribution evolution (Schnakenberg, 1976) that describes probability flows between different states of a system, in our case, different masses of dense cores described by the DCMF .
The DCMF declines at low masses because dense cores cannot form below the Jeans mass () where the gravitational contraction cannot overcome the gas thermal pressure (Jeans, 1901) or the turbulent support (Hennebelle & Chabrier, 2008). In our model, dense core progenitors are generated across a finite range of masses according to the initial probability distribution called source function .
The three processes change the number of dense cores in a cell over one time step : some cores of mass grow and enter the cell, some cores of mass grow and leave the cell, and new cores are created in this cell. The accretion rate given by the growth equation (Eq. 3) is the same for all dense cores in a given cell. Putting these contributions together:
(5) 
As the evolution runs for a long time, converges to a constant shape even if the number of dense cores and the total mass of the system keep growing. We take the dynamic equilibrium limit, so that . We also now go to the continuousmass and continuoustime description, such that , while . In that case, we replace the difference between the two accretion terms in Eq. 5 above with a differential.
(6) 
We can substitute and take the steadystate limit where . This is actually a very weak assumption: we do not presuppose any specific functional law for either or , we only need to assume that the ratio of those two rates is asymptotically constant. As both accretion and generation of new dense core progenitors are governed by the same physical processes, we expect them to slow down at the same rate. With this simplification, we obtain the growth equation:
(7) 
The Eq. 7 acts as a filter (a linear functional map) that converts an initial density fluctuation spectrum for dense core progenitors into a DCMF . Note that all timedependent quantities, such as or independently standing (now it only normalizes the differential source function) have canceled out, thus the DCMF shape does not depend on how the system slows down in time. This equation preserves the normalization because . An exact analytical solution is only possible for some simple functional shape of , but a number of numerical solutions are presented for illustrative purposes in Fig. 3. At high masses, for any choice of the DCMF develops the same powerlaw tail with an exponent defined only by the fractal mass distribution properties, while at low masses it essentially preserves the source function shape, with a smooth transition in between. In order to match observations, we take a lognormal source function of a form (Chabrier, 2003; Hennebelle & Chabrier, 2008; Hennebelle & Falgarone, 2012). Its maximum lies at . In Fig. 4 we pickthe and that best resemble the observed distribution.
While the equation 7 allows us to accurately match the observed DCMF, its relevance and generality stretches beyond that. To calculate the highmass tail of the distribution analytically, in that limit we can neglect a rapidly decaying (e.g. a decaying exponent, Gaussian, or lognormal). Then, Eq. 7 becomes homogeneous and has an analytic solution of a form , regardless of the input source function shape:
(8) 
6. Nonlinear accretion
Maschberger (2013) considers the dense core growth through accretion that is in general both nonlinear and stochastic. The accretion rate derived above is linear (), although in general, nonlinear cases are also possible with and . Accretion can be either sublinear (, e.g. Bonnell et al., 2001) or superlinear (, e.g. in Bondi, 1952).
Maschberger (2013) describes the process of competitive accretion using the Stratonovich stochastic calculus formalism in order to predict the mass distribution of dense cores when the accretion rate is partially random and fluctuating. The possible fluctuations need to be restricted to be exclusively nonnegative to rule out the mass loss. In our theory, we account for the stochasticity of accretion using the master equation (Eq. 7). Since in our case, the only possible transition from each bin in mass is to the next bin, our calculation is also restricted to nonnegative fluctuations. Therefore, by plugging an alternative accretion rate , the tail part of the DCMF becomes:
(9) 
Here the constant is now dimensionful for . For sublinear accretion, the DCMF decays at high mass as stretched exponential, while for the superlinear growth it results in a shallow power law .
Here the nonlinear accretion is directly analogous to the nonlinear preferential attachment in network science (Krapivsky et al., 2000; Newman, 2010). The sublinear preferential attachment similarly results in a stretched exponential type distribution. The superlinear preferential attachment results in a situation where a few network nodes accumulate a macroscopic fraction of all edges in the network. This issue is recognized in Maschberger (2013) as an “explosion” of dense core masses in the absence of noise. Since there is no observational evidence of starforming regions, where the entire mass is dominated by a few very massive stars, the explosive growth scenario seems unrealistic.
Because the phenomenology of accretion in the fractal media is not clear, we restrict the further analysis to the simple case of linear accretion following from our model and given by Eq. 7.
7. Discussion and Summary
7.1. Dependence of the slope on input parameters
Having obtained an analytic expression (Eq. 8) for the powerlaw exponent , we can now explore how it behaves as we vary the three possible parameters of the system , and . An important thing to stress is that none of the three parameters bears any dimensional units. On one hand, this is due to the term “scalefree distribution”: in the highmass tail there is no characteristic mass or scale that defines the shape of the distribution (as opposed to other functional forms, such as normal, lognormal, or exponential). On the other hand, these parameters are directly related to fundamental scaling laws of statistical physics relevant on a broad spectrum of length, time, and mass scales.
By substituting the observed value (Falgarone et al., 1991; Elmegreen & Falgarone, 1996), we obtain . The uniform density distribution of gas parcels corresponds to and . The Salpeter value corresponds to . This fractal dimensionality is predicted and observed in a number of physical systems governed by Brownian processes such as the diffusionlimited aggregation (Meakin, 1983) known to take place for dust in the ISM (Planck Collaboration et al., 2011). Whether or not an actual Brownian process stays behind the value is beyond the scope of our work, however, finding how observed physical properties of the ISM may affect its fractal dimensionality at the lowmass end might provide a clue to our understanding of IMF variations.
We also notice that if , the second term in the expression turns into unity and yields . For generic values of the parameters, the numerator in the last fraction of Eq. 8 represents the variety of choice of parcels for a given dense core (Fig. 2 (right)), whereas the denominator represents the variety of choices of dense cores for a parcel (Fig. 2 (left)). If both are distributed in space following the same fractal dimension, then the patterns of parcels attaching to dense cores no longer depend on the spatial coordinates. Effectively, for the spatiality of the problem cancels out and it reduces to the “regular” preferential attachment process (avoiding the doublecounting of network edges as in Barabási & Albert, 1999). The expression also becomes independent of , thus removing the necessity of our model choice to weigh the accretion probabilities by gravitational accelerations or gravitational potentials.
A specific possible value of at is described by Hopkins (2013) as “a generic scalefree distribution, allotting equal mass to each equal logarithmic interval in mass.” The actual value obtained by Hopkins is equal to plus a small addition coming from various effects related to the properties of the turbulent ISM and magnetic fields. In our theory, that addition appears naturally from considering a twophase medium, i.e. different spatial distributions of dense cores and parcels.
Our result favors the unimodal IMF shape over the bimodal. The broken power law (Kroupa, 2001) is acceptable as a fitting approximation for the smooth transition between the lowmass decline and the highmass powerlaw tail. This agrees with the conclusions drawn from numerical simulations (Elmegreen, 1997, 2002) of the fractal ISM evolution. Clauset et al. (2009) specifically discuss the difficulties of fitting power laws and other fattail distributions to empirical data and point out that it is often hard to distinguish which model represents the data better.
Then, given no evidence that the turbulence induced ISM fractal dimensionality varies across different starforming regions, the variation of remains the only channel to explain possible IMF nonuniversality. The two hard limits are () and ().
For a system with a finite number of cores we estimate the mass ratio of the largest to smallest cores in the powerlaw regime as . This explains the observed correlation between the most massive star mass and the total star cluster mass that we can calculate for any specific solution of Eq. 7 (Kroupa et al., 2013).
Larson (1992) attempted to relate the ISM fractal dimensionality to the IMF shape by assuming that the entire mass from some fragment of the molecular cloud surrounding a core accretes onto it. Then, the IMF powerlaw exponent becomes equal to the fractal dimensionality. In our model, however, we do not make the assumption that every core grows by accreting matter from a distinct region of the cloud but rather consider the competitive accretion (or preferential attachment) in order to account for overlapping basins of gravitational attraction.
A filamentary distribution of parcels will correspond to . This will change the convergence of integrals in Eq. 2–3 but will still result in a powerlaw mass function. In principle, it is possible to introduce a scaledependent fractal structure where and/or change at some characteristic scale . This will, however, make the calculations bulky and will also introduce additional free parameters so that the solution behavior will be more difficult to investigate and explain analytically.
7.2. Bottomheavy mass functions
The lowmass star dominated bottomheavy IMF suggested by recent observations (van Dokkum & Conroy, 2010; Cappellari et al., 2012) has a slope at certain masses steeper than the asymptotic value (Fig. 5). The logarithmic slope is given by . We derive it directly from the master equation (Eq. 7) in a selfreferential form, without solving it for any specific :
(10) 
is always nonnegative because dense cores in our model are never destroyed, and is nonnegative as a probability distribution. Asymptotically, their ratio because has exponential or faster decay and is a “slower” power law. Thus, the logarithmic slope from above, and it can never become steeper than unless the nonnegativity condition is violated. We solved Eq. 7 for a fiducial source function that is negative for some masses (Fig. 5) in order to illustrate how a bottomheavy DCMF can be established. Because only depends on the fundamental scaling exponents , and , and because it serves as a hard lower bound on the DCMF slope, bottomheavy mass functions are ruled out by our theory for the linear accretion regime.
This conclusion comes into tension with the results that suggest a bottomheavy IMF shape in elliptical galaxies (van Dokkum & Conroy, 2010; Cappellari et al., 2012) with . It is worth mentioning, that those conclusions have recently been challenged by statistical data analysis (Smith, 2014; Clauwens et al., 2015), observations of extragalactic Xray binaries (Peacock et al., 2014), and strong gravitational lensing (Smith et al., 2015). However, one has to keep in mind that it is impossible to observationally measure the IMF slope at masses () in old stellar populations, because stars at that mass range have already evolved into remnants. Therefore, the unimodal steep IMF slope cannot be excluded as a solution satisfying both our theory and observations, if future studies explain how the parcel fractal density dimensionality depends on a galaxy mass or the ISM metal content.
Also, if we admit variations of the unimodal IMF slope across different starforming regions in the same galaxy, the observational appearance of an IMF to be bottomheavy becomes plausible for composite stellar populations (e.g. galaxies formed by major dry mergers). That can happen if, for example, a combined IMF shape is determined for a stellar system that consists of several building blocks having different intrinsic IMF slopes and comparable masses. Then, the combined stellar distribution will be dominated by lowmass stars preferably from a bottomheavy building block, while its highmass end will be defined by a shallow (topheavy) IMF stellar component. This explains why until now, no standalone star cluster or a starforming complex with a bottomheavy IMF has been found with the same integrated light spectral diagnostics as those used to derive the bottomheavy IMF shape in giant earlytype galaxies (van Dokkum & Conroy, 2010, 2011). Stellar systems that can be reasonably well represented by simple stellar populations, such as ultracompact dwarf galaxies and massive globular clusters, exhibit stellar masses corresponding to the lowmass IMF slopes between Kroupa and Salpeter (Chilingarian et al., 2011; van Dokkum & Conroy, 2011; Podorvanyuk et al., 2013).
7.3. Summary
We presented a simple analytic approach that addresses the following major points formulated in the introduction and left unexplained by existing IMF theories:

The scalefree distribution of dense cores by mass is established by the process of preferential attachment (competitive accretion) of mass onto dense cores.

When the system mass grows, the distribution shape asymptotically stabilizes.

The powerlaw exponent depends only on two parameters, and , fractal dimensions of the turbulent ISM and accreting gas, directly connected to their fundamental physical properties.
Our theory relies on the qualitative description of the supersonic turbulence that follows from the basic Kolmogorov theory. The real structure of the supersonic turbulent flow in the ISM might be different and will potentially affect our results. However, if the density distribution can still be described as fractal, it will only affect the powerlaw slope as suggested by Eq. 8. We explain the bimodality of the Kroupa IMF as a result of a twocomponent fitting of the intrinsically unimodal distribution in the transition region () between a power law at high masses and a declining part at low masses (Fig. 4). By our calculation of the lower bound on the logarithmic slope, in the transition region it should never be steeper than that at higher masses, therefore we rule out a bottomheavy IMF shape for any single standalone starforming region.
References
 Albert et al. (1999) Albert, R., Jeong, H., & Barabási, A.L. 1999, Nature, 401, 130
 Alves et al. (2007) Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17
 André et al. (2010) André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102
 Barabási & Albert (1999) Barabási, A.L., & Albert, R. 1999, Science, 286, 509
 Barabási et al. (2002) Barabási, A. L., Jeong, H., Néda, Z., et al. 2002, Physica A Statistical Mechanics and its Applications, 311, 590
 Barthélemy (2011) Barthélemy, M. 2011, Phys. Rep., 499, 1
 Basu et al. (2015) Basu, S., Gil, M., & Auddy, S. 2015, MNRAS, 449, 2413
 Bondi (1952) Bondi, H. 1952, MNRAS, 112, 195
 Bonnell & Bate (2006) Bonnell, I. A., & Bate, M. R. 2006, MNRAS, 370, 488
 Bonnell et al. (2001) Bonnell, I. A., Clarke, C. J., Bate, M. R., & Pringle, J. E. 2001, MNRAS, 324, 573
 Cappellari et al. (2012) Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2012, Nature, 484, 485
 Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
 Chilingarian et al. (2011) Chilingarian, I. V., Mieske, S., Hilker, M., & Infante, L. 2011, MNRAS, 412, 1627
 Clauset et al. (2009) Clauset, A., Shalizi, C. R., & Newman, M. E. J. 2009, SIAM Review, 51, 661
 Clauwens et al. (2015) Clauwens, B., Schaye, J., & Franx, M. 2015, MNRAS, 449, 4091
 Dorogovtsev & Mendes (2002) Dorogovtsev, S. N., & Mendes, J. F. F. 2002, Advances in Physics, 51, 1079
 Elmegreen (1997) Elmegreen, B. G. 1997, ApJ, 486, 944
 Elmegreen (2002) —. 2002, ApJ, 564, 773
 Elmegreen & Falgarone (1996) Elmegreen, B. G., & Falgarone, E. 1996, ApJ, 471, 816
 Elmegreen & Scalo (2004) Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211
 Falgarone et al. (1991) Falgarone, E., Phillips, T. G., & Walker, C. K. 1991, ApJ, 378, 186
 Hennebelle & Chabrier (2008) Hennebelle, P., & Chabrier, G. 2008, ApJ, 684, 395
 Hennebelle & Chabrier (2011) Hennebelle, P., & Chabrier, G. 2011, in IAU Symposium, Vol. 270, Computational Star Formation, ed. J. Alves, B. G. Elmegreen, J. M. Girart, & V. Trimble, 159–168
 Hennebelle & Chabrier (2013) —. 2013, ApJ, 770, 150
 Hennebelle & Falgarone (2012) Hennebelle, P., & Falgarone, E. 2012, A&A Rev., 20, 55
 Hopkins (2013) Hopkins, P. F. 2013, MNRAS, 430, 1653
 Hoyle & Lyttleton (1939) Hoyle, F., & Lyttleton, R. A. 1939, Proceedings of the Cambridge Philosophical Society, 35, 592
 Jeans (1901) Jeans, J. H. 1901, Royal Society of London Proceedings Series I, 68, 454
 Jeong et al. (2001) Jeong, H., Mason, S. P., Barabási, A.L., & Oltvai, Z. N. 2001, Nature, 411, 41
 Kolmogorov (1941) Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301
 Krapivsky et al. (2000) Krapivsky, P. L., Redner, S., & Leyvraz, F. 2000, Physical Review Letters, 85, 4629
 Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
 Kroupa (2002) —. 2002, Science, 295, 82
 Kroupa et al. (2013) Kroupa, P., Weidner, C., PflammAltenburg, J., et al. 2013, The Stellar and SubStellar Initial Mass Function of Simple and Composite Populations, ed. T. D. Oswalt & G. Gilmore, 115
 Larson (1992) Larson, R. B. 1992, MNRAS, 256, 641
 Lombardi et al. (2015) Lombardi, M., Alves, J., & Lada, C. J. 2015, A&A, 576, L1
 Maschberger (2013) Maschberger, T. 2013, MNRAS, 436, 1381
 Meakin (1983) Meakin, P. 1983, Physical Review Letters, 51, 1119
 Merton (1968) Merton, R. K. 1968, Science, 159, 56
 Newman (2010) Newman, M. 2010, Networks: An Introduction (OUP Oxford)
 Newman (2004) Newman, M. E. J. 2004, PNAS, 101, 5200
 Padoan & Nordlund (2002) Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870
 Peacock et al. (2014) Peacock, M. B., Zepf, S. E., Maccarone, T. J., et al. 2014, ApJ, 784, 162
 Planck Collaboration et al. (2011) Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2011, A&A, 536, A24
 Podorvanyuk et al. (2013) Podorvanyuk, N. Y., Chilingarian, I. V., & Katkov, I. Y. 2013, MNRAS, 432, 2632
 Price (1976) Price, D. D. S. 1976, Journal of the American Society for Information Science, 292
 Ravasz et al. (2002) Ravasz, E., Somera, A. L., Mongru, D. A., Oltvai, Z. N., & Barabási, A.L. 2002, Science, 297, 1551
 Sadavoy et al. (2010) Sadavoy, S. I., Di Francesco, J., Bontemps, S., et al. 2010, ApJ, 710, 1247
 Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
 Schmidt et al. (2010) Schmidt, W., Kern, S. A. W., Federrath, C., & Klessen, R. S. 2010, A&A, 516, A25
 Schnakenberg (1976) Schnakenberg, J. 1976, Reviews of Modern Physics, 48, 571
 Silk & Takahashi (1979) Silk, J., & Takahashi, T. 1979, ApJ, 229, 242
 Simon (1955) Simon, H. A. 1955, Biometrika, 42, 425
 Smith (2014) Smith, R. J. 2014, MNRAS, 443, L69
 Smith et al. (2015) Smith, R. J., Lucey, J. R., & Conroy, C. 2015, MNRAS, 449, 3441
 Spaans & Silk (2000) Spaans, M., & Silk, J. 2000, ApJ, 538, 115
 Sreenivasan et al. (1989) Sreenivasan, K. R., Ramshankar, R., & Meneveau, C. 1989, Royal Society of London Proceedings Series A, 421, 79
 Tarasov (2005) Tarasov, V. E. 2005, Physics Letters A, 336, 167
 van Dokkum & Conroy (2010) van Dokkum, P. G., & Conroy, C. 2010, Nature, 468, 940
 van Dokkum & Conroy (2011) —. 2011, ApJ, 735, L13
 Yakovenko & Rosser (2009) Yakovenko, V. M., & Rosser, J. B. 2009, Rev. Mod. Phys., 81, 1703
 Yule (1925) Yule, G. U. 1925, Philosophical Transactions of the Royal Society B, 213, 21
 Zinnecker (1982) Zinnecker, H. 1982, Annals of the New York Academy of Sciences, 395, 226