EXPLAINING THE STELLAR INITIAL MASS FUNCTION WITH THE THEORY OF SPATIAL NETWORKS

Six decades ago the stellar initial mass function (IMF) was derived from star counts (Salpeter 1955) as a scale-free power law (${dN}/{dm}\approx {m}^{\alpha };$ $\alpha =-2.35$) with more frequent low-mass stars than high-mass 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 understand star formation physics. From observations, the IMF shape appears to be universal across different star-forming regions (Kroupa 2002). Resembling a unimodal (Salpeter 1955) or bimodal (Kroupa 2001) power law or a log-normal distribution with a power-law tail (Chabrier 2003), it sharply declines at the low end at masses $\lt 1/12$ solar mass (M_⊙). 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 M_⊙, offset by a factor of ∼4 to higher masses compared to the stellar IMF, illustrating a ∼25% gas-to-stars 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 low-mass IMF decline (Silk & Takahashi 1979) but does not reproduce a power-law high-mass 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 to date (Hennebelle & Chabrier 2008, 2013) explains the observed overall IMF shape by analyzing the evolution of density fluctuations in a self-gravitating turbulent ISM in the presence of a magnetic field. Similarly to other gravoturbulent theories, in order to reproduce a power-law 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 log-normal 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 log-normal shape and vary across different star-forming 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 power-law 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 log-normal. However, this might be the result of the mathematical simplification applied (the top-hat filtering in the Fourier space and the consequent Taylor expansion) that transforms an arbitrary function into a log-normal 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 non-log-normal or non-analytic PDF shapes becomes cumbersome and not very straightforward. Schmidt et al. (2010) demonstrate analytically and numerically that gravoturbulent theories can reproduce the power-law 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 power-law IMF tail, they fail to match its observed exponent $(\alpha =-2.35)$ and cannot be applied to non-clustered star formation. A recent IMF theory from that family (Basu et al. 2015) generates a log-normal distribution with a power-law 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 power-law 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 high-mass power-law 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) can a scale-free power-law 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 power-law 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 protostar formation. Therefore, it applies to starless dense cores. The principal result of our theory is an analytic expression for the power-law 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.

We describe a system of dense cores as an array of masses M_i distributed according to some time-dependent function $p(m,t)$. At a given moment, the system has a total mass ${M}_{\mathrm{tot}}$ and a total number of dense cores N. The mass distribution shape is governed by 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 dm per one unit of time dt attaches to one of the cores. Regardless of which core it joins, the total mass of the system grows by dm, hence ${{dM}}_{\mathrm{tot}}={dm}$. The accretion summed over the entire system (e.g., the global mass growth of the system) happens at some characteristic rate K = dm/dt. Over a time step dt every mass bin m can only be affected by bins within dm = Kdt 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 f(m). This function is the only mechanism that increases N, therefore ${dN}/{dt}={\int }_{0}^{\infty }f(m){dm}=F$. We stress that the K and F values are, in general, time-dependent. 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 $p(m,t)$ with $m\in \{{dm},2{dm},3{dm},\ldots \}$, discrete time $t\in \{{dt},2{dt},3{dt},\ldots \}$, and normalization ${\sum }_{m=1}^{\infty }p(m,t)=1$. We are interested in its long term behavior, such that $p(m,t)\underset{t\to \infty }{\to }p(m)$.

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 self-similar or fractal density distribution in a turbulent flow (Sreenivasan et al. 1989) with the predicted fractal dimensionality $D=7/3=2.33$. It stays in agreement with measurements obtained from observations of giant molecular clouds (Falgarone et al. 1991; Elmegreen & Falgarone 1996) ($D=2.2\mbox{--}2.4$) and in laboratory studies of turbulence (Sreenivasan et al. 1989) (D = 2.35). Numerical simulations of the ISM evolution with an input fractal density field (Elmegreen 1997, 2002) yield power-law mass distributions of overdensities that correspond to dense cores in star-forming regions.

We describe ISM as a two-phase 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 two-phase 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 l 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 d corresponds to the average separation between particles of gas or dust. Presumably $l\gg d$, 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 $r\gg l$ or $r\gg d$, uniform fractal scaling is observed, i.e., the mass confined in a sphere of radius r grows as $M(r)\sim {r}^{{D}_{m}}$, where D_m 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 D_m = 3, but for the fractal distribution ${D}_{m}\lt 3$. Consequently, the average density in a sphere $\rho \sim M(r)/{r}^{3}$ is not constant, but is scale-dependent. The homogeneity means that the fractal power-law 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. Self-similar (or fractal) scaling appears only on scales when discrete features of characteristic scales l and d blur out.

Since the dense core distribution traces the initial Kolmogorov-like supersonic turbulent distribution, we take it to have ${D}_{m}=7/3\simeq 2.33$. Confining our system to a box of linear size L and normalizing its mass to ${M}_{\mathrm{tot}}$, the mass confined in a thin concentric sphere becomes ${dM}(r)={M}_{\mathrm{tot}}{D}_{m}{r}^{{D}_{m}-1}{L}^{-{D}_{m}}{dr}$. The fractal dimension of parcels $2\lt {D}_{p}\leqslant 3$ 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 ${dn}(r)={D}_{p}{r}^{{D}_{p}-1}{L}^{-{D}_{p}}{dr}$.

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 Figure 1 right) will still be defined by "point-like" dense cores.

Zoom In Zoom Out Reset image size

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 power-law 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, 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 distance-dependent, 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 Figure 1, right panel). We set the parcel accretion probability by a given core proportional to the initial gravitation acceleration toward it (Figure 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 (Figure 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 j to join an existing dense core i is:

$$\begin{eqnarray}&&{\pi }_{j\to i}=\displaystyle \frac{1}{{c}_{j}}\displaystyle \frac{{M}_{i}}{| {r}_{i}-{r}_{j}{| }^{\delta }}\end{eqnarray} \tag{ 1 }$$

where M_i is the mass of the ith core, $| {r}_{i}-{r}_{j}|$ is the Euclidean distance between the points, δ is the gravitational law exponent ($\delta =2$ in the three-dimensional space). A possible alternative choice is $\delta =1$, 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:

$$\begin{eqnarray}{c}_{j} & = & \displaystyle \sum _{i}\displaystyle \frac{{M}_{i}}{| {r}_{i}-{r}_{j}{| }^{\delta }}={\displaystyle \int }_{{l}_{i}}^{L}{dM}(r){r}^{-\delta }\\ & = & \displaystyle \frac{{M}_{\mathrm{tot}}{D}_{m}}{{L}^{{D}_{m}}}\displaystyle \frac{{L}^{{D}_{m}-\delta }-{l}_{j}^{{D}_{m}-\delta }}{{D}_{m}-\delta }\end{eqnarray} \tag{ 2 }$$

where l_j is a distance from a new parcel to the nearest dense core. Assuming ${D}_{m}\gt \delta$ and ${l}_{j}\ll L$, we can neglect l_j 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.

Zoom In Zoom Out Reset image size

From examining the integrals in Equations (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 ${r}^{{D}_{m}-1}/{r}^{\delta }$ or ${r}^{{D}_{p}-1}/{r}^{\delta }$ and vanishes in the limit of large r.

The average dense core mass increase per time step is a probability weighted sum over all possible positions where parcels emerge.

$$\begin{eqnarray}{{dM}}_{i} & = & {dm}\displaystyle \sum _{{\boldsymbol{r}}}\quad N(\overrightarrow{{\boldsymbol{r}}}){\pi }_{j\to i}\\ & = & {dm}{\displaystyle \int }_{{l}_{i}}^{L}{D}_{p}{dr}\displaystyle \frac{{r}^{{D}_{p}-1}}{{L}^{{D}_{p}}}\displaystyle \frac{{L}^{\delta }}{{M}_{\mathrm{tot}}}\displaystyle \frac{{D}_{m}-\delta }{{D}_{m}}\displaystyle \frac{{M}_{i}}{{r}^{\delta }}=\beta \displaystyle \frac{{M}_{i}}{{M}_{\mathrm{tot}}}{dm}\end{eqnarray} \tag{ 3 }$$

We notice the difference between our calculated accretion rate for an individual core ${{dM}}_{i}/{dt}\propto {M}_{i}$ and the classical Bondi–Hoyle–Lyttleton model (Hoyle & Lyttleton 1939; Bondi 1952) of spherical accretion (${dM}/{dt}\propto {M}^{2}$). This inconsistency is trivially explained by two facts: (a) we assume a fractal distribution of the infalling matter so that a thin spherical layer no longer contains the mass $\rho {R}^{2}{dR}$ 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 $\beta =(1-\delta /{D}_{m})/(1-\delta /{D}_{p})\lt 1$. 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 t₀, a dense core has mass ${M}_{i,0}$ and the global mass growth rate is quasi-constant (${dm}\sim {dt}$ and ${M}_{\mathrm{tot}}\sim t$), then it grows in time according to a sublinear power law:

$$\begin{eqnarray}&&{M}_{i}={M}_{i,0}{\left(\displaystyle \frac{t}{{t}_{0}}\right)}^{\beta }.\end{eqnarray} \tag{ 4 }$$

In reality, this law does not hold because the global mass growth rate K 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 $t\to \infty$ limit, which in turn is not affected by the specific time dependence of K and F.

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 $p(m,t)$.

The DCMF declines at low masses because dense cores cannot form below the Jeans mass (M_J) 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 f(m).

The three processes change the number of dense cores in a cell m over one time step dt: some cores of mass m − dm grow and enter the cell, some cores of mass m grow and leave the cell, and $f(m){dt}$ new cores are created in this cell. The accretion rate given by the growth equation (Equation (3)) is the same for all dense cores in a given cell. Putting these contributions together:

$$\begin{eqnarray}&&p(m,t+{dt})(N+{dN})=p(m,t)N\\ &&+\beta \displaystyle \frac{(m-{dm})p(m-{dm},t)N}{{M}_{\mathrm{tot}}}-\beta \displaystyle \frac{{mp}(m,t)N}{{M}_{\mathrm{tot}}}+f(m){dt}.\end{eqnarray} \tag{ 5 }$$

As the evolution runs for a long time, p(m) 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 $p(m,t)\to p(m)$. We also now go to the continuous-mass and continuous-time description, such that ${dm},{dt}\to 0$, while ${dm}/{dt}=K$. In that case, we replace the difference between the two accretion terms in Equation (5) above with a differential.

$$\begin{eqnarray}&&p(m){dN}=-\beta {dm}\displaystyle \frac{N}{{M}_{\mathrm{tot}}}\displaystyle \frac{d}{{dm}}({mp}(m))+f(m){dt}.\end{eqnarray} \tag{ 6 }$$

We can substitute ${dN}/{dt}=F$ and take the steady-state limit where ${dm}/{dN}\to {M}_{\mathrm{tot}}/N$. This is actually a very weak assumption: we do not presuppose any specific functional law for either K = dm/dt or F = dN/dt, 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:

$$\begin{eqnarray}&&p(m)+\beta \displaystyle \frac{d}{{dm}}({mp}(m))=\displaystyle \frac{f(m)}{F}.\end{eqnarray} \tag{ 7 }$$

Equation (7) acts as a filter (a linear functional map) that converts an initial density fluctuation spectrum f(m) for dense core progenitors into a DCMF p(m). Note that all time-dependent quantities, such as K or independently standing F (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 ${\int }_{0}^{\infty }p(m){dm}={\int }_{0}^{\infty }f(m){dm}/F=1$. An exact analytical solution is only possible for some simple functional shape of f(m), but a number of numerical solutions are presented for illustrative purposes in Figure 3. At high masses, for any choice of f(m) the DCMF develops the same power-law 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 log-normal source function of a form $f(m)\propto \tfrac{1}{m}\mathrm{exp}\left(-\tfrac{(\mathrm{ln}m-\mu )}{2{\sigma }^{2}}\right)$ (Chabrier 2003; Hennebelle & Chabrier 2008; Hennebelle & Falgarone 2012). Its maximum lies at ${M}_{\mathrm{max}}=O(\mathrm{exp}(\mu ))=O({M}_{{\rm{J}}})$. In Figure 4 we pick the μ and σ that best resemble the observed distribution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While Equation (7) allows us to accurately match the observed DCMF, its relevance and generality stretches beyond that. To calculate the high-mass tail of the distribution analytically, in that limit we can neglect a rapidly decaying f(m) (e.g., a decaying exponent, Gaussian, or log-normal). Then, Equation (7) becomes homogeneous and has an analytic solution of a form $p(m)={{Cm}}^{\alpha }$, regardless of the f(m) input source function shape:

$$\begin{eqnarray}&&p(m)={{Cm}}^{\alpha },\alpha =-\left(1+\displaystyle \frac{1}{\beta }\right)=-\left(1+\displaystyle \frac{1-\delta /{D}_{p}}{1-\delta /{D}_{m}}\right).\end{eqnarray} \tag{ 8 }$$

Maschberger (2013) considers the dense core growth through accretion that is in general both nonlinear and stochastic. The accretion rate derived above is linear (${{dM}}_{i}/{dt}\propto {M}_{i}$), although in general, nonlinear cases are also possible with ${{dM}}_{i}/{dt}\propto {M}_{i}^{a}$ and $a\ne 1$. Accretion can be either sublinear ($a\lt 1$, e.g., Bonnell et al. 2001) or superlinear ($a\gt 1$, e.g., a = 2 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 non-negative, to rule out the mass loss. In our theory, we account for the stochasticity of accretion using the master equation (Equation (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 non-negative fluctuations. Therefore, by plugging in an alternative accretion rate ${{dM}}_{i}/{dt}\propto {M}_{i}^{a}$, the tail part of the DCMF becomes:

$$\begin{eqnarray}&&p(m)={{Cm}}^{-a}\mathrm{exp}\left(-\displaystyle \frac{1}{\beta }\displaystyle \frac{{m}^{1-a}}{1-a}\right)\end{eqnarray} \tag{ 9 }$$

Here the constant β is now dimensionful for $a\ne 1$. For sublinear accretion, the DCMF decays at high mass as a stretched exponential, while for the superlinear growth it results in a shallow power law $\simeq {m}^{-a}$.

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 star-forming 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 fractal media is not clear, we restrict further analysis to the simple case of linear accretion following from our model and given by Equation (7).

7.1. Dependence of the Slope α on Input Parameters

7.2. Bottom-heavy Mass Functions

The low-mass star dominated bottom-heavy IMF suggested by recent observations (van Dokkum & Conroy 2010; Cappellari et al. 2012) has a slope at certain masses steeper than the asymptotic value α (Figure 5). The logarithmic slope is given by $d\mathrm{log}p(m)/d\mathrm{log}m=({dp}(m)/{dm})(m/p(m))$. We derive it directly from the master equation (Equation (7)) in a self-referential form, without solving it for any specific f(m):

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{log}p(m)}{d\mathrm{log}m}=\displaystyle \frac{{dp}(m)}{{dm}}\displaystyle \frac{m}{p(m)}=\alpha +\displaystyle \frac{f(m)}{\beta {Fp}(m)}\end{eqnarray} \tag{ 10 }$$

f(m) is always non-negative because dense cores in our model are never destroyed, and p(m) is non-negative as a probability distribution. Asymptotically, their ratio $f(m)/p(m)\to 0$ because f(m) has exponential or faster decay and p(m) is a "slower" power law. Thus, the logarithmic slope $d\mathrm{log}p(m)/d\mathrm{log}m\to \alpha$ from above, and it can never become steeper than α unless the non-negativity condition is violated. We solved Equation (7) for a fiducial source function that is negative for some masses (Figure 5) in order to illustrate how a bottom-heavy DCMF can be established. Because α only depends on the fundamental scaling exponents δ, D_m, and D_p, and because it serves as a hard lower bound on the DCMF slope, bottom-heavy mass functions are ruled out by our theory for the linear accretion regime.

Zoom In Zoom Out Reset image size

This conclusion comes into tension with the results that suggest a bottom-heavy IMF shape in elliptical galaxies (van Dokkum & Conroy 2010; Cappellari et al. 2012) with $\alpha \approx 3$. It is worth mentioning that those conclusions have recently been challenged by statistical data analysis (Smith 2014; Clauwens et al. 2015), observations of extragalactic X-ray 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 ($M\gt 0.7{M}_{\odot }$) 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 D_p depends on a galaxy mass or the ISM metal content.

Also, if we admit variations of the unimodal IMF slope across different star-forming regions in the same galaxy, the observational appearance of an IMF to be bottom-heavy 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 low-mass stars preferably from a bottom-heavy building block, while its high-mass end will be defined by a shallow (top-heavy) IMF stellar component. This explains why until now, no stand-alone star cluster or a star-forming complex with a bottom-heavy IMF has been found with the same integrated light spectral diagnostics as those used to derive the bottom-heavy IMF shape in giant early-type 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 low-mass IMF slopes between Kroupa and Salpeter (Chilingarian et al. 2011; van Dokkum & Conroy 2011; Podorvanyuk et al. 2013).

7.3. Summary

We are grateful to our anonymous reviewers whose feedback helped us to improve the manuscript. I.C. acknowledges support from the Telescope Data Center, Smithsonian Astrophysical Observatory. His theoretical and observational IMF studies are supported by the Russian Foundation for Basic Research projects 15-52-15050 and 15-32-21062, and the President of Russian Federation grant MD-7355.2015.2. The idea of applying network science formalism to the IMF theory was developed by the authors during the annual Chamonix workshop in 2014 supported by the Russian Science Foundation project 14-22-00041. The authors appreciate fruitful discussions with and useful suggestions from J. Silk, P. Hennebelle, G. Mamon, P. Kroupa, S. Mieske, M. Kurtz, I. Zolotukhin, C. Lada, and M. E. J. Newman.