<figure><img src="image_1.png" /></figure>

HAL
open science

# On the mass function of stars growing in a flocculent medium

Th. Maschberger

► To cite this version:

Th. Maschberger. On the mass function of stars growing in a flocculent medium. Monthly Notices of the Royal Astronomical Society, 2013, 436, pp.1381-1389. 10.1093/mnras/stt1656. insu-03614199

HAL Id: insu-03614199

https://insu.hal.science/insu-03614199

Submitted on 20 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# On the mass function of stars growing in a flocculent medium

Th. Maschberger*

UJF-Grenoble 1 / CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG) UMR 5274, Grenoble F-38041, France

Accepted 2013 September 2. Received 2013 September 2; in original form 2013 May 31

## ABSTRACT

Stars form in regions of very inhomogeneous densities and may have chaotic orbital motions. This leads to a time variation of the accretion rate, which will spread the masses over some mass range. We investigate the mass distribution functions that arise from fluctuating accretion rates in non-linear accretion, $\dot{m} \propto m^{\alpha}$. The distribution functions evolve in time and develop a power-law tail attached to a lognormal body, like in numerical simulations of star formation. Small fluctuations may be modelled by a Gaussian and develop a power-law tail $\propto m^{-\alpha}$ at the high-mass side for $\alpha > 1$ and at the low-mass side for $\alpha < 1$. Large fluctuations require that their distribution is strictly positive, for example, lognormal. For positive fluctuations the mass distribution function develops the power-law tail *always* at the high-mass hand side, independent of $\alpha$ larger or smaller than unity. Furthermore, we discuss Bondi–Hoyle accretion in a supersonically turbulent medium, the range of parameters for which non-linear stochastic growth could shape the stellar initial mass function, as well as the effects of a distribution of initial masses and growth times.

**Key words:** accretion, accretion discs – turbulence – stars: formation – stars: luminosity function, mass function.

## 1 INTRODUCTION

Star-forming regions typically show a very inhomogeneous structure with large variations in the gas density due to turbulence and filaments. Thus, accretion rates of forming stars, depending on gas density, will show fluctuations. Sufficiently substantial fluctuations in the accretion rate will spread out the starting masses after accretion occurred, which affects, like many other effects, the stellar initial mass function, the distribution of stellar masses at their ‘birth’. The accretion rate of a point mass in a homogeneous medium follows $\dot{m} \propto \rho m^2$, which is Bondi–Hoyle–Lytleton accretion (Edgar 2004). In a flocculent medium, the density variations change the deterministic Bondi–Hoyle–Lytleton accretion into a stochastic process, which is multiplicative, as $\dot{m}$ depends on $m$, and non-linear, as $\dot{m} \propto m^2$. In this paper, we investigate the distribution function that arises from such a non-linear multiplicative stochastic process.

A linear stochastic process, fragmentation, has been employed for some time to explain the shape of the stellar initial mass function at low masses. The initial mass function has at low masses a lognormal shape, but develops a power-law tail at high masses (Kroupa 2002; Chabrier 2003). Larson (1973), Elmegreen & Mathieu (1983) and Zinnecker (1984) modelled fragmentation as a sequence of discrete fragmentation steps, during each the mass is reduced by a fraction of itself. Applying the central limit theorem leads to a lognormal distribution function after a sufficient number of fragmentation steps.

Because the fragment mass is chosen to depend linearly on the mass of the cloud this is a linear process. Fragmentation alone seems not to embrace the whole star formation process, as the initial mass function has a power-law tail at the massive end. This deviation from lognormality has been explained, amongst other ideas, by competitive accretion, which is non-linear accretion of the fragments without fluctuations in the accretion rate (Larson 1978; Zinnecker 1982; Bonnell et al. 1997, 2001a,b; Bate, Bonnell & Bromm 2003). The power-law tail of the initial mass function then arises from the scatter of the initial masses from which accretion starts.

With (linear) fragmentation and (non-linear) accretion two processes encounter each other that have opposite sign. Fragmentation leads (on average) to a reduction of mass, whereas accretion increases mass. One can therefore speculate, whether it is possible to model the main part of the star formation process as a single stochastic process that is non-linear. This has already been attempted by Marcus (1968), who adopted the random splitting model of Filipov (1961) (who extended the work of Kolmogorov 1941¹). Marcus (1968) assumed a time-discrete non-linear stochastic process. An important aspect of the results by Marcus (1968) is that simultaneously the distribution of the number of fragments *and* the mass distribution of the fragments are derived. We investigate a time-continuous non-linear stochastic process, generalizing the work of Marcus (1968) in that respect, but do not derive a distribution for the number of fragments.

<footnote>* E-mail: thomas.maschberger@obs.ujf-grenoble.fr</footnote>

<footnote n="1">¹ In English: Kolmogorov (1992).</footnote>
In linear or non-linear random fragmentation, it is usually assumed that the fraction that is lost during a step follows a Gaussian distribution. This has a curious side effect: the mass distribution function at the end does not vanish above the initial mass. But nothing can fragment to a mass larger than its initial mass. The reason for this lies in the Gaussian distribution used to describe the fragment distribution; there is a non-vanishing probability that some fragments have negative mass. In stellar growth (which is 'negative fragmentation'), the distribution of accretion rates is due to the fluctuations in the gas density. The gas density cannot be smaller than zero. Therefore, we investigate the case of a non-linear stochastic process with a strictly positive distribution of fluctuations as well. The difference to Gaussian noise can be quite substantial.

In order to model the star formation process, it is, besides growth, also necessary to consider the initial mass from which accretion starts and the time how long accretion lasts. Several models covering all these aspects are present in the literature, for example: Basu & Jones (2005) assume a lognormal distribution of initial masses with growth $\propto m$ or $\propto m^{2/3}$ and an exponential distribution of growth times. Bate & Bonnell (2005) consider growth with a lognormal distribution of accretion rates without a mass dependence, no distribution of initial masses, and also an exponential distribution of growth times. Related to this mass-independent growth with lognormal accretion rates is the discussion of stable distributions in Cartwright & Whitworth (2012). In a series of papers, Myers (2000, 2008, 2009, 2011, 2012) published increasingly elaborate models of the star formation process. Their main components are an exponential distribution of growth times, accretion with a mass-independent and a (non-linearly) mass-dependent contribution, and a constant initial mass (a distribution of initial masses is considered in Myers 2009). The model of Dib et al. (2010) contains a distribution of initial masses originating from gravoturbulent accretion (lognormal with a power-law tail) from which stars grow with a mass-dependent accretion rate that is exponentially dampened in time (cf. also Dib, Kim & Shadmehri 2007). Effects of fluctuating accretion rates are to our knowledge not yet considered in the literature. In Section 5.3, we describe the effects of a distribution of initial masses and a distribution of growth times on the mass distribution arising from fluctuating accretion rates in a stochastic growth process.

Before embarking on the investigation of non-linear growth with a distribution of the accretion rates we would like to mention that our results may have more applications than stellar growth. Non-linear stochastic processes appear in many other contexts, particularly in the context of turbulence. It is, for example, found that the probability distribution function (pdf) of the gas density in a turbulent medium shows power-law tails in case of a polytropic index other than unity (Passot & Vázquez-Semadeni 1998). Given that isothermal turbulence is assumed to be a linear multiplicative process, the occurrence of the power-law tails hints at a non-linear process. Similarly, the pdf of velocity in a turbulent medium show power-law tails (e.g. Krumholz, McKee & Klein 2006), which again may be a sign of non-linearity. Gravity also leads to non-linear stochastic processes, power-law tails appear in self-gravitating turbulent gas if it is isothermal (e.g. Klessen & Burkert 2000). Even stellar dynamics can be seen as a stochastic process (Chandrasekhar 1943).

The outline of this paper is as follows: In the next section, we prepare the necessary mathematical prerequisites with the help of linear stochastic growth. Section 3 contains results for non-linear stochastic growth where the random increments are assumed to be Gaussian. Strictly positive fluctuations are considered in Section 4. In Section 5, we discuss Bondi–Hoyle accretion in a supersonically

turbulent medium, the typical parameter ranges that would affect the stellar initial mass function and the effects of a distribution of initial masses and growth times. The usual summary concludes the paper.

## 2 A TIME-CONTINUOUS STOCHASTIC FORMULATION OF LINEAR GROWTH OR FRAGMENTATION

We start with the simplest form of mass-dependent growth where the accretion rate is linear in mass. Without fluctuations this is described by the differential equation

$$ \frac{dm}{dt} = Am. \qquad (1) $$

The quantity A accounts for all the constants that are involved. Suppose now that the star grows in a flocculent medium. Then, the fluctuating gas density will cause fluctuations in the accretion rate, or, more specifically, A will be varying. The mass dependence remains unchanged.

Mathematically, the fluctuations can be introduced by changing equation (1) from an ordinary differential equation to a stochastic differential equation,

$$ dm = m(a dt + b dW). \qquad (2) $$

Now A is split into two terms, $a dt$, which describes the mean of A, and $b dW$, which describes the fluctuations around the mean. We assume that the fluctuations stem from a normal (or Gaussian) distribution,

$$ \mathcal{N}(x; \mu, \sigma) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(x-\mu)^2}{2\sigma^2}}, \qquad (3) $$

with zero mean and variance $dt$, so that the distribution of $dW$ is

$$ p(dW) = \mathcal{N}(dW; \mu = 0, \sigma = \sqrt{dt}). \qquad (4) $$

For $dt = 1$, the distribution of $dW$ has unit variance. Multiplying $b$ and $dW$ is equivalent to scaling the variance to the desired amount of fluctuations. Alternatively, the distribution of accretion rates can be written as

$$ p(A) = \mathcal{N}(A; \mu = a, \sigma = b). \qquad (5) $$

For integrating equation (2), we need to establish how the integral $\int_0^t dW$ is to be interpreted. The integral can be approximated with the limit of $\sum_{i=1}^n dW_i$, where all the $dW_i$ are independently drawn from a normal distribution. The normal distribution has the property that the sum of two normal variates follows again a normal distribution (in other words, the normal distribution is infinitely visible). The parameters for the sum variate are $\mu = \mu_1 + \mu_2$ and $\sigma^2 = \sigma_1^2 + \sigma_2^2$. Therefore, $\sum_{i=1}^n dW_i$ will again obey a normal distribution with $\sigma^2 = \sum_{i=1}^n \sigma_i^2 = \sum_{i=1}^n (\sqrt{dt})^2 = t$. The limit $dt \to 0$ and $n \to \infty$ gives the integral $W_i = \int_0^t dW$, which is a normally distributed random number with zero mean and variance $t$. A more rigorous derivation this can be found in the literature on stochastic differential equations (e.g. Øksendal 2002).

The stochastic calculus is not unique, there are two ways of defining it, the Itô and the Stratonovich calculus. Depending on the nature of the fluctuations, the one or the other are more appropriate. In our case, the fluctuations in density are not caused by the growth process, but by external effects (e.g. turbulence or a chaotic motion of the growing star). Such external fluctuations require the Stratonovich calculus (van Kampen 2006), which preserves the standard rules of calculus.
With this we can solve equation (2) by integrating it to

$$ \log m - \log m_0 = at + bW_t, \qquad (6) $$

where $m_0$ is the initial mass and growth starts at $t = 0$.

For a population of growing stars, the distribution of $m(t)$ can be found by noting that

$$ at + bW_t \sim \mathcal{N}(\mu = at, \sigma = b\sqrt{t}) \qquad (7) $$

($\sim$ denotes ‘is distributed as’). Thus,

$$ \log m - \log m_0 \sim \mathcal{N}(\mu = at, \sigma = b\sqrt{t}), \qquad (8) $$

or,

$$ p(m) = \frac{1}{m} \frac{1}{\sqrt{2\pi}} \frac{1}{b\sqrt{t}} e^{-\frac{1}{2} \frac{(\log m - \log m_0 - at)^2}{b^2 t}}, \qquad (9) $$

where the factor $1/m$ comes from the transformation from $\log m$ to $m$. The distribution of $m(t)$ is, like in the discrete case, a lognormal distribution.

The distribution function of $m(t)$ (equation 9) is also defined for values $<m_0$, as in the discrete case. $dm$ can have values $<-m$, as $dW$ can reach up to $-\infty$. Thus, the mathematical formulation of the growth/fragmentation process, both continuous and discrete, is not strictly correct. Nevertheless, for a small level of fluctuations the probability for $m < m_0$ can be very small. Note that $m(t)$ is always larger than zero.

Another puzzling property of equation (9) lies in the expectation value of $m(t)$,

$$ E(m_t) = m_0 e^{at + \frac{b^2 t}{2}}. \qquad (10) $$

Here, an additional term $\frac{b^2}{2}t$ appears compared to the solution $m(t) = m_0 e^{At}$ of the deterministic differential equation (1). This ‘spurious drift’ is a consequence of using the Stratonovich stochastic calculus.

## 3 NON-LINEAR GROWTH

The accretion of stars is not linear in mass. Depending on the environment $\dot{m} \propto m^2$ (Bondi–Hoyle accretion) or $\dot{m} \propto m^{2/3}$ (gas-dominated potential; Bonnell et al. 2001b). Small fluctuations in the accretion rate due to a flocculent density or random stellar motions can be described neglecting that $dW$ can take negative values, the corresponding stochastic differential equation is

$$ dm = m^{\alpha} (a dt + b dW). \qquad (11) $$

This can be solved proceeding analogously to the previous section. In this section, we discuss the solution of equation (11) and its properties. For an application to stellar accretion with large fluctuations, we need modify equation (11), see the next section.

The solution for $m(t)$ is

$$ m(t) = \left( (1 - \alpha) \left( \frac{m_0^{1-\alpha}}{1-\alpha} + at + bW_t \right) \right)^{1/(1-\alpha)}, \qquad (12) $$

which reduces for $b = 0$ to the solution of the deterministic growth law.

For exponents $\alpha > 1$, we have to account for the fact that the solutions are exploding in the deterministic case ($b = 0$), infinite masses are reached within a finite time,

$$ t_{\text{ex}} = \frac{m_0^{1-\alpha}}{a(\alpha - 1)}. \qquad (13) $$

When fluctuations are added to the accretion rate, then there is no single explosion time any more. Depending on the particular fluctuations that are encountered by a growing star it may explode earlier or later, at some random time. $m(t)$ of equation (12) becomes undefined for large $W_t$. Thus, we have to require

$$ W_t < -\frac{1}{b} \left( \frac{m_0^{1-\alpha}}{1-\alpha} + at \right) =: u_t \qquad (14) $$

if we consider only the not yet exploded particles. At some time $t$ the fraction of the population that is not yet exploded is

$$ f_{\text{nex,G}}(t) = \Phi(u_t; 0, \sqrt{t}), \qquad (15) $$

($\Phi$ is the cumulative normal distribution), because $W_t \sim \mathcal{N}(0, \sqrt{t})$. In nature, the exploding solutions are suppressed because an infinite reservoir from which material could be accreted does not exist. This should be accounted for in the growth model, for example in the line of logistic growth. Massive stars (which would become an exploding solution) exercise a strong feedback as well, which likewise suppresses exceedingly large accretion rates. Unfortunately, accounting for a finite reservoir and feedback is beyond the scope of this paper.

The mass distribution function for the population is then

$$ p_G(m, t) = \frac{1}{f_{\text{nex,G}}(t)} \frac{1}{m^\alpha} \frac{1}{\sqrt{2\pi}} \frac{1}{b\sqrt{t}} e^{-\frac{1}{2} \frac{(\frac{1}{1-\alpha} (m^{1-\alpha} - m_0^{1-\alpha}) - at)^2}{b^2 t}}. \qquad (16) $$

The factor $1/f_{\text{nex,G}}(t)$ normalizes $p_G$ as probability by accounting for the exploded solutions if $\alpha > 1$. For $\alpha < 1$, it is not required, $f_{\text{nex,G}}(t) = 1$. Fig. 1 shows $p_G$ for various values of $\alpha$. The other parameters are $a = 1$, $b = 1$, $t = 1$ and $m_0 = 1$. For $\alpha = 1$ (top-left panel), the solution corresponds to a lognormal distribution, as discussed in Section 2. For $\alpha \neq 1$ (and $\alpha \neq 0$), a power-law tail appears at one side of $p_G$. If $\alpha > 1$ (top-right panel), the power-law tail $\propto m^{-\alpha}$ develops at the high-mass end of $p_G$. The power law is indicated by the dashed line. At smaller masses, $p_G$ behaves similar to a lognormal distribution which is shown by the dotted curve for comparison (the parameters for the lognormal distribution are chosen to follow $p_G$). For exponents smaller than unity in the growth law, a power-law tail develops at the left-hand side of the

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Comparison of the distribution functions $p_G(m, t)$ for various $\alpha$. The parameters are $a = 1$, $b = 1$, $t = 1$, $t_0 = 0$ and $m_0 = 1$. For $\alpha = 1$ (top left) the solution is a lognormal distribution. In the case of $\alpha > 1$ (top right), a power-law tail at high masses develops (dashed). The dotted line is a lognormal distribution with parameters chosen to fit the lower part of $p_G$. If $0 < \alpha < 1$ (bottom left) a falling power-law tail at low masses develops, whereas if $\alpha < 0$ the power-law increases.</figcaption></figure>
peak at small masses. The other side of the distribution resembles a lognormal distribution. If $0 < \alpha < 1$, the power law is decreasing, which gives the distribution a rather peculiar appearance. For $\alpha < 0$, i.e. for example $\dot{m} \propto 1/m$, the low-mass power law increases, so that $p_G$ appears mirrored to the case of $\alpha > 1$. Although $\alpha < 0$ may not be relevant for stellar accretion, it might be suitable to describe the stochastic process for the distribution of velocities in supersonic turbulence. Krumholz et al. (2006) found in numerical simulations that the distribution of velocities has a power-law tail at the left-hand side.

We do not discuss at this point the time evolution and the behaviour of the mean mass, but postpone it to the next section for a comparison with the case of strictly positive fluctuations.

## 4 STRICTLY POSITIVE NOISE

The fluctuations in the stellar accretion rates originate in the variations of the gas density. If there are only small variations, then their description with a mean density modulated by a Gaussian is sufficiently accurate. However, if the gas density is very flocculent, as for example in supersonic turbulence, then the use of a Gaussian introduces an undesirable side effect: the random variate describing the fluctuations can become so large that it exceeds the mean density. A strongly underdense region would then be assigned a negative density, which is physically impossible. In supersonic turbulence, the gas density pdf has been found to follow a lognormal distribution,

$$p_{i\mathcal{N}}(x) = \frac{1}{x} \frac{1}{\sqrt{2\pi}\sigma_i} e^{-\frac{1}{2} \frac{(\log x - \mu_i)^2}{\sigma_i^2}} \quad (17)$$

(e.g. Vazquez-Semadeni 1994; Passot & Vázquez-Semadeni 1998; Nordlund & Padoan 1999). We could use the lognormal distribution to describe the fluctuations. The lognormal distribution is, like the normal distribution, infinitely divisible (Thorin 1977) and consequently is suitable to describe fluctuations in a stochastic differential equation. However, there are no formulas for the sum distribution of two lognormal variates, so the stochastic integral cannot be solved analytically.

Thus, we use for practical reasons an approximation of the lognormal distribution, the inverse Gaussian distribution (invGauss),

$$p_{\text{invGauss}}(x; v, \lambda) = \left(\frac{\lambda}{2\pi x^3}\right)^{1/2} e^{-\frac{\lambda(x-v)^2}{2v^2x}}. \quad (18)$$

It has expectation value $E(x) = v$ and variance $Var(x) = v^3/\lambda$. A comparison of the invGauss and the lognormal distribution with the same expectation value and variance is shown in Fig. 2. Compared to the invGauss the lognormal has somewhat heavier tails.

As model of strictly positive fluctuations, we choose that an infinitesimal fluctuation

$$\mathrm{d}iG_{a,b} \sim \mathrm{invGauss} \left( v = a \mathrm{d}t, \lambda = \frac{a^3}{b^2} (\mathrm{d}t)^2 \right). \quad (19)$$

It has the mean value $a \mathrm{d}t$ and variance $b^2 \mathrm{d}t$, like fluctuations from a normal distribution. For the approximation of the integral by a sum, we note that the sum of two invGauss random numbers with the same $v$ and $\lambda$ follows again an invGauss distribution with $v' = 2v$ and $\lambda' = 2^2\lambda$. With this follows that

$$iG_{a,b,t} := \sum_{i=1}^{n} \mathrm{d}iG_{a,b} \sim \mathrm{invGauss} \left( v = an \mathrm{d}t, \lambda = \frac{a^3}{b^2} (n \mathrm{d}t)^2 \right), (20)$$

the sum of infinitesimal fluctuations follows again an invGauss distribution. The mean of $iG_{a,b,t}$ is $an \mathrm{d}t = at$ and the variance is

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Comparison of a lognormal distribution (blue dashed curve) and an inverse Gaussian distribution (black solid) as approximation. Both have the same mean (1) and variance (0.6). The lognormal distribution is over a wide range well approximated by the inverse Gaussian distribution.</figcaption></figure>

$b^2(n \mathrm{d}t)^2 = b^2 t$. This also corresponds to the results from the case with a normal distribution. Now we are able to perform the limits $\mathrm{d}t \to 0$ and $n \to \infty$ which gives the integral $\int \mathrm{d}iG_{a,b}$.

With this we can pose the mass-dependent growth equation with strictly positive fluctuations,

$$\mathrm{d}m = m^{\alpha} \mathrm{d}iG_{a,b}. \quad (21)$$

Using standard calculus, we find for $\alpha \neq 1$ that

$$m(t) = \left( (1 - \alpha) \left( \frac{m_0^{1-\alpha}}{1-\alpha} + iG_{a,b,t} \right) \right)^{1/(1-\alpha)} \quad (22)$$

(If $\alpha = 1$, the term $m_0^{1-\alpha}/(1-\alpha)$ has to be replaced by $\log m_0$.) Again, we have to consider the exploded solutions if $\alpha > 1$ and need to require $\frac{m_0^{1-\alpha}}{1-\alpha} + iG_{a,b,t} < 0$ so that $m(t)$ is not infinite. The not exploding fraction is then

$$f_{\text{nex iG}} = P_{\text{invGauss}} \left( -\frac{m_0^{1-\alpha}}{1-\alpha}, v = at, \lambda = \frac{a^3}{b^2} t^2 \right), \quad (23)$$

where $P_{\text{invGauss}}$ is the invGauss cumulative distribution function. For $\alpha \le 1$, nothing explodes and $f_{\text{nex iG}} = 1$. With $m(t)$ and $f_{\text{nex iG}}$, we can write the mass function,

$$p_{iG}(m, t) = \frac{1}{f_{\text{nex, iG}}} \frac{1}{m^{\alpha}} \frac{1}{\sqrt{2\pi}} \frac{t}{b} \left(\frac{X}{a}\right)^{-3/2} e^{-\frac{1}{2} \frac{a}{b^2} \frac{(X - at)^2}{X}}, \quad (24)$$

where

$$X := \begin{cases} \frac{m^{1-\alpha}}{1-\alpha} - \frac{m_0^{1-\alpha}}{1-\alpha} & \alpha > 1 \\ \log m - \log m_0 & \alpha = 1 \end{cases}. \quad (25)$$

Fig. 3 shows the behaviour of $p_{iG}$ for different $\alpha$ with $a = 1$, $b = 1$, $m_0 = 1$, $t = 1$ and $t_0 = 0$. These are the same parameters as in Fig. 1, which shows $p_G$. For $\alpha > 1$, the positive fluctuations lead to the same power law $\propto m^{-\alpha}$ as Gaussian fluctuations. If $\alpha \le 1$, the positive fluctuations change the behaviour compared to Gaussian fluctuations, now appears power-law-like behaviour at high masses as well. This can be characterized with the ‘exponent’ function ($p_{iG} \propto m^{S(m)}$),

$$S(m) = \frac{\mathrm{d} \log p_{iG}}{\mathrm{d} \log m} \quad (26)$$

$$= -\alpha + m^{1-\alpha} \left( -\frac{1}{2} \frac{a}{b^2} - \frac{3}{2} \frac{1}{X} + \frac{1}{2} \frac{a^3}{b^2} t^2 - \frac{1}{X^2} \right). \quad (27)$$
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Distribution functions $p_{iG}(m, t)$ for various $\alpha$ (compare with Fig. 1). The parameters are $a = 1$, $b = 1$, $t = 1$, $t_0 = 0$ and $m_0 = 1$. $p_{iG}$ has a power-law tail for $\alpha \ge 1$ and behaves even for $\alpha < 1$ similar to a power law over some mass range. The dashed lines show power laws whose exponents are discussed in the text.</figcaption></figure>

$\alpha = 1$ leads to a strict power-law behaviour at high masses with exponent $-\alpha - \frac{1}{2}\frac{a}{b^2}$. However, as visible in Fig. 3 this exponent occurs only as a limiting case and can set in not until very large masses. At smaller masses $p_{iG}$ is steeper.

If $\alpha < 1$, then $S(m)$ has $-\infty$ as limit for $m \to \infty$, $p_{iG}$ decays at high masses. Nevertheless, $p_{iG}$ can for some combinations of $a$, $b$ and $t$ be well described by a power law over some mass range above the peak. The terms $-\frac{1}{2}\frac{a}{b^2}m^{1-\alpha}$ and $\frac{1}{2}\frac{a^3}{b^2}\frac{m^{1-\alpha}}{X^2}$ cancel each other around

$$m_{\beta} = ((1-\alpha)at + m_{0}^{1-\alpha})^{1/(1-\alpha)}, \qquad (28)$$

where

$$p_{iG} \propto m^{-\alpha-\beta} \qquad (29)$$

with

$$\beta = \frac{3}{2} \frac{1}{1 - \frac{m_0}{m_{\beta}} 1^{-\alpha}}. \qquad (30)$$

Fig. 4 displays the time-evolution of $p_{iG}$ (solid), as well as $p_G$ (dashed) for $a = 1$ and the ratios $b/a = 0.1$, $b/a = 0.5$ and $b/a = 1.0$. Both distribution develop a power-law tail over time, which happens fastest if the amount of fluctuations is not too large (i.e. $b/a = 0.5$). If the fluctuations are very large, then large accretion rates are occurring only very rarely, as they do if the fluctuations are very small. Thus, the power-law tail is slower populated. Both distributions also shift in mass range over time, which is slower for larger $b/a$. As $p_{iG}$ does not allow for ‘negative’ growth it always moves to higher masses. The Gaussian fluctuations in $p_G$ make the mass shift so slow for $b/a = 0.5$ that the peak effectively does not move. For $b/a = 1$, the shift has reversed, the peak of $p_G$ moves towards masses smaller than $m_0$.

The time development of the expectation value of $p_{iG}$ and $p_G$ (calculated numerically) is given in Fig. 5. The parameters are $\alpha = 2.1$, $a = 1$, $b = 0.5$ and $m_0 = 1$. Also displayed is the deterministic solution as dotted curve which explodes at $t_{ex} = 0.909$. Note that for $1 < \alpha \le 2$ the mean of $p_G$ and $p_{iG}$ is infinity, because the mean for a power law with such an exponent is infinite. For $t < t_{ex}$, the mean mass with a fluctuating accretion rate is larger than expected from deterministic growth. This is similar to the spurious drift introduced by the stochastic formulation discussed in Section 2. Solutions can explode but are not accounted for in the distribution function, so

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Time evolution of the mass functions $p_{iG}$ (solid curve) and $p_G$ (dashed curve) for $\alpha = 2$, $a = 1$, $m_0 = 1$ and $t_0 = 0$. In contrast to the Gaussian noise there is no ‘growth’ below $m_0 = 1$ for strictly positive noise. For $b/a = 0.1$, $p_{iG}$ lies outside the plotting area at $t = 5$ and $t = 10$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Time evolution of the mean of $p_{iG}$ (solid curve). The blue dashed curve shows the expectation for $p_G$ and the red dotted line the expectation of the deterministic case.</figcaption></figure>

that both averages are finite even for $t > t_{ex}$. The mean $p_{iG}$ appears to explode, in contrast to $p_G$, albeit slower than deterministic.

# 5 APPLICATION TO STAR FORMATION

## 5.1 Bondi–Hoyle accretion in a medium with supersonic turbulence

Krumholz et al. (2006) have studied the distribution of accretion rates for Bondi–Hoyle accretion in a supersonically turbulent medium. They performed hydrodynamical calculations of accretion on to stationary point masses in a medium that is not self-gravitating. They found that the normalized accretion rate, $\dot{m}/m_0$, can be fitted
<table><caption>**Table 1.** Parameters for the distribution of accretion rates of Bondi–Hoyle accretion in a turbulent medium (calculated from the results of Krumholz et al. 2006). The first two columns give Mach number *M* and the ratio of Bondi Radius *r*<sub>B</sub> to the extent of the region *l*. The next three columns are for the dimensionless case where *ρ* = 1, *G* = 1, *c*<sub>s</sub> = 1 and *m* = 1. The other columns are for physical units where *c*<sub>s</sub> = 220 m s<sup>-1</sup>.</caption><thead><tr><th rowspan="2">M</th><th rowspan="2">log r<sub>B</sub>/l</th><th colspan="3">ρ = 1</th><th colspan="3">ρ = 1000 M<sub>☉</sub> pc<sup>-3</sup></th><th colspan="2">ρ = 10<sup>4</sup> M<sub>☉</sub> pc<sup>-3</sup></th></tr><tr><th>a<sub>0</sub></th><th>b<sub>0</sub></th><th>b<sub>0</sub>/a<sub>0</sub></th><th>a</th><th>b</th><th>b/a</th><th>a</th><th>b/a</th></tr><tr><td></td><td></td><td colspan="3">dimensionless</td><td colspan="3">(1/M<sub>☉</sub>/Myr)</td><td colspan="2"></td></tr></thead><tbody><tr><td>3</td><td>-5.00</td><td>0.17</td><td>1.8</td><td>11</td><td>1.8</td><td>6</td><td>3.4</td><td>18</td><td>1.1</td></tr><tr><td>3</td><td>-3.00</td><td>0.17</td><td>1.8</td><td>11</td><td>1.8</td><td>6</td><td>3.4</td><td>18</td><td>1.1</td></tr><tr><td>3</td><td>-1.00</td><td>0.071</td><td>0.42</td><td>6</td><td>0.75</td><td>1.3</td><td>1.8</td><td>7.5</td><td>0.57</td></tr><tr><td>3</td><td>1.00</td><td>0.0029</td><td>0.014</td><td>5</td><td>0.031</td><td>0.045</td><td>1.5</td><td>0.31</td><td>0.46</td></tr><tr><td>5</td><td>-5.00</td><td>0.07</td><td>5.4</td><td>77</td><td>0.73</td><td>17</td><td>24</td><td>7.3</td><td>7.5</td></tr><tr><td>5</td><td>-3.00</td><td>0.066</td><td>4.9</td><td>74</td><td>0.7</td><td>16</td><td>23</td><td>7</td><td>7.2</td></tr><tr><td>5</td><td>-1.00</td><td>0.028</td><td>0.93</td><td>33</td><td>0.3</td><td>3</td><td>10</td><td>3</td><td>3.2</td></tr><tr><td>5</td><td>1.00</td><td>0.0017</td><td>0.024</td><td>15</td><td>0.017</td><td>0.078</td><td>4.5</td><td>0.17</td><td>1.4</td></tr><tr><td>10</td><td>-5.00</td><td>0.016</td><td>10</td><td>619</td><td>0.17</td><td>33</td><td>191</td><td>1.7</td><td>60</td></tr><tr><td>10</td><td>-3.00</td><td>0.016</td><td>9.1</td><td>588</td><td>0.16</td><td>30</td><td>181</td><td>1.6</td><td>57</td></tr><tr><td>10</td><td>-1.00</td><td>0.0077</td><td>2.3</td><td>300</td><td>0.081</td><td>7.5</td><td>93</td><td>0.81</td><td>29</td></tr><tr><td>10</td><td>1.00</td><td>0.000 83</td><td>0.072</td><td>87</td><td>0.0087</td><td>0.23</td><td>23</td><td>0.087</td><td>8.5</td></tr></tbody></table>

by a lognormal distribution with parameters depending on the Mach number and the ratio of the Bondi radius, *r*<sub>B</sub>, to the extent of the region, *l*. The normalization constant is

$$ \dot{m}_0 = 4\pi\rho \frac{G^2}{(M c_s)^3} m^2, \quad (31) $$

with the gas density $\rho$, Mach number *M*, sound speed $c_s$ and the gravitational constant *G*. Krumholz et al. (2006) give their results in dimensionless units where $\rho = 1$, $G = 1$ and $c_s = 1$, and used $m = \frac{13}{32}$. In Table 1, we give their results adapted to our formalism, where the average accretion rate *a* and its standard deviation *b* do not depend on mass. Typically, the ratio *b/a* in dimensionless units is much larger than unity and can reach values of a few hundred.

Table 1 gives also *a* and *b* in physical units (solar masses and million years) for a typical star forming region. We assume a temperature of 10 K, which corresponds to a sound speed of 220 m s<sup>-1</sup>. The Bondi radius at that temperature is $r_B = Gmc_s^{-2} = 0.09 m/M_{\odot} pc$, for a typical region size of 1 pc and 1 M<sub>☉</sub> mass log(r<sub>B</sub>/l) = -2.4. Although *b/a* becomes smaller with the scaling, it is still typically larger than unity. This large level of fluctuations will only spread out the initial masses over a range of final masses if accretion lasts for several million years, or the initial masses are ≫1 M<sub>☉</sub>, or the average gas density is very high.

The large fluctuations in the accretion rate do not change the result of Krumholz, McKee & Klein (2005) finding that Bondi–Hoyle accretion should not generate large mass increases in typical star-forming regions. For a low average accretion rate, strong accretion events are too rare to have an effect. However, Bonnell & Bate (2006) argue that in the central regions of star-forming regions the average gas densities are several orders of magnitude higher than the average gas density used by Krumholz et al. (2005), in which case accretion does lead to mass growth and, as shown above, a spreading out of the initial masses.

## 5.2 The required amount of fluctuations for an impact of non-linear stochastic growth

Forming stars are embedded in the flocculent environment of their natal cloud and will accrete from it, but in contrast to Bondi–Hoyle accretion the ambient medium is self-gravitating as well. Bonnell

et al. (2001b) argue that in a gas-dominated potential the accretion rate should follow $m^{2/3}$, whereas in a stellar-dominated potential (uncorrelated velocities of gas and stars) classical Bondi–Hoyle accretion $\propto m^2$ occurs. Fluctuations in the accretion rate of forming stars can be generated either by the turbulence of the gas, or by the self-gravity of the gas cloud which generates filaments, or by the chaotic orbital motion of the forming star, which will bring it into regions of different gas density. Fluctuating accretion can determine the shape of the stellar initial mass function. In this section, we investigate the range of values that parameters can take so that accretion has an impact. A measurement of the fluctuations in numerical studies that account for the self-gravity of the gas is required to answer what are the ramifications of accretion in the star formation process.

Accretion can shape the initial mass function in two ways: a power-law tail appears and the distribution becomes wider than the distribution of the initial masses. An additional distribution of growth times will contribute to both as well. Accretion will only have an impact on mass functions if the initial masses are spread over a sufficiently wide range and populate a power-law tail. This occurs if $p_{iG}$ has a rounded triangular shape.

In Fig. 6, we show in the left plot the probability density, $p_{iG}(m, t = 1)$, and in the right plot the cumulative distribution for a range of average accretion rates *a* and level of fluctuations *b/a*. $\dot{m} \propto m^{2/3}$ is shown as solid curve and $\dot{m} \propto m^2$ as dotted curve. Starting mass is $m_0 = 1$ and growth time is $t = 1$. $\alpha = 2/3$ constrains both *a* and *b/a* to be of order unity, otherwise $p_{iG}$ develops a power-law part too shallow to be consistent with the initial mass function. $\alpha = 2$ allows for a larger range of parameters.

For both $\alpha$ the average accretion rate *a* needs to be larger than unity in order to sufficiently spread out the initial masses. Although a power-law tail appears and becomes stronger with increasing *b/a*, it does not contain many stars (cumulative distribution in the power-law tail already $\approx 1$).

For a mean accretion rate $a \gtrsim 1$, the situation depends on the level of fluctuations. If $b/a \ll 1$, then growth is effectively deterministic, which less affects growth with $\alpha = 2$. If the level of fluctuations is very large ($b/a \gg 1$), $p_{iG}$ is very peaked at $m_0$ and the power-law tail is not strongly populated as well, unless the average accretion
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Probability density $p_{iG}$ (left) and cumulative distribution $P_{iG}$ (right) for various $a$ and $b/a$ with $m_0 = 1$ and $t = 1$. The solid curve is for $\alpha = 2/3$ and the dotted curve for $\alpha = 2$. As reference each panel shows $p_{iG}$ and $P_{iG}$ for $\alpha = 2, a = 1$ and $b/a = 1$ in grey.</figcaption></figure>

<table><caption>**Table 2.** Some values for scaling of the initial mass and average accretion rate or time.</caption><thead><tr><th colspan="10">Scaling with initial mass: k<sub>m</sub> = (m'<sub>0</sub>/m<sub>0</sub>)<sup>1-α</sup>, a' = k<sub>m</sub>a and b' = k<sub>m</sub>b</th></tr><tr><th></th><th></th><th>m'<sub>0</sub>/m<sub>0</sub></th><th>0.01</th><th>0.10</th><th>0.20</th><th>0.50</th><th>1.00</th><th>2.00</th><th>5.00</th><th>10.00</th></tr></thead><tbody><tr><td>α = 2</td><td>k<sub>m</sub> =</td><td>100.00</td><td>10.00</td><td>5.00</td><td>2.00</td><td>1.00</td><td>0.50</td><td>0.20</td><td>0.10</td></tr><tr><td>α = 2/3</td><td>k<sub>m</sub> =</td><td>0.22</td><td>0.46</td><td>0.58</td><td>0.79</td><td>1.00</td><td>1.26</td><td>1.71</td><td>2.15</td></tr></tbody><thead><tr><th colspan="10">Scaling with average accretion rate a' = k<sub>a</sub>a, b' = &radic;<span style="text-decoration: overline">k<sub>a</sub>b</span>, t' = t/k<sub>a</sub></th></tr></thead><tbody><tr><td>a' = k<sub>a</sub>a =</td><td>0.10</td><td>0.20</td><td>0.50</td><td>1.00</td><td>2.00</td><td>5.00</td><td>10.00</td></tr><tr><td>t' = t/k<sub>a</sub></td><td>10.00</td><td>5.00</td><td>2.00</td><td>1.00</td><td>0.50</td><td>0.20</td><td>0.10</td></tr><tr><td>b' = &radic;<span style="text-decoration: overline">k<sub>a</sub>b</span> =</td><td>0.32</td><td>0.45</td><td>0.71</td><td>1.00</td><td>1.41</td><td>2.24</td><td>3.16</td></tr><tr><td>b'/a' =</td><td>3.16</td><td>2.24</td><td>1.41</td><td>1.00</td><td>0.71</td><td>0.45</td><td>0.32</td></tr></tbody></table>

rate is very high. For a lognormal distribution, the ratio between median and mean is $R = \exp(-\sigma_t^2/2) = ((b/a)^2 + 1)^{-1/2}$, depending only on the ratio of standard deviation to mean. When $\sigma_l = 2$ ($b/a = 7.32$), this means that the median accretion rate is only $\approx 1/10$ of the mean accretion rate. For comparison, $\sigma_l = 1$ corresponds to $b/a = 1.31$.

With increasing $a$ and $b/a$ the exploding fraction has to be considered for $\alpha = 2$. For $b/a = 0.1$, it is not relevant at $t = 1$ ($\le 10^{-25}$), as for $a = 0.1$. For $b/a = 1$, the exploding fraction is $3 \times 10^{-5}$ for $a = 1$, and 10 and 15 per cent for $a = 10$ and 100, respectively. For $b/a = 10$, a significant fraction of the seeds explode (34 per cent at $a = 1$, 45 per cent at $a = 10$ and 100). These parameter combinations would require some feedback mechanism to prohibit explosions.

Scaling to a different initial mass leaves the shape of $p_{iG}(m, t)$ unchanged for the same $t$, if $a$ and $b$ are scaled by multiplying both with the scaling factor $k_m = (m'_0/m_0)^{1-\alpha}$. Values for the scaling factor for some typical initial masses are given in Table 2. For $\alpha = 2$, they span a very wide range, whereas for $\alpha = 2/3$ they lie around unity.

Scaling the mean accretion rate does not preserve the shape of $p_{iG}$, unless the ratio $b/a$ and the time-scale are changed. Correspond-

ingly, if the same shape of $p_{iG}$ should be achieved in a different time, both $a$ and $b$ require scaling. Here, the scaling does not preserve the ratio $b/a$. Some values for the scaled parameters are given in Table 2.

The above findings can be summarized in the rule of thumb that both the average accretion rate $a$ and the level of fluctuations $b/a$ have to be of order unity for unit initial mass and unit growth time in order to sufficiently populate the power-law tail. Nevertheless, there are other parameter combinations that may also lead to the desired behaviour of $p_{iG}$. Furthermore, if the initial mass, the average accretion rate and the time are given in physical units an observed ratio $b/a$ may be far from unity.

## 5.3 Effects of a distribution of initial masses and growth times

Fig. 7 shows the effects of a distribution of initial masses and a distribution of growth times. We use a lognormal distribution of initial masses (equation 17) with $\mu_l = -3.15$ and $\sigma_l = 0.55$, corresponding to a mean mass $\bar{m}_0 = 0.05 M_{\odot}$. $p(m_0)$ is shown as dashed curve in Fig. 7. The distribution function for an ensemble
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Ensemble distribution functions for a population growing from a distribution of initial masses (dashed) with different distribution of growth times. $\alpha = 2/3$, $a = 2$ and $b = 1$. In the top panel, all stars accrete for the same time, the middle panel uses an exponential distribution of growth times and the bottom panel a uniform distribution. The dotted line shows a standard system mass function.</figcaption></figure>

growing from a distribution of initial masses is given by

$$p_E(m, t) = \int p(m_0)p_{iG}(m, t; \alpha, a, b, m_0) dm_0. \quad (32)$$

The top panel of Fig. 7 shows $p_E$ for $\alpha = 2/3$, $a = 3 M_{\odot}^{1/3} Myr^{-1}$ and $b = \sqrt{3} M_{\odot}^{1/3} Myr^{-1}$ at $t = 0.3$ Myr, which corresponds to the dimensionless case of $a = 1$, $b = 1$ and $t = 1$. Note that here no explosions occur as $\alpha < 1$. The peak of $p_E$ is significantly shifted compared to the distribution of the $m_0$. Above $m_\beta = 0.45 M_{\odot}$ (calculated using $\bar{m}_0$) a power-law part $\propto m^{-1.8}$ appears ($\beta = 1.11$) which starts to decay for larger masses. At $m = 1 M_{\odot}$ the exponent is $-2.6$ and $-4.4$ at $m = 10 M_{\odot}$. For comparison, we also show the system IMF as dotted (Chabrier 2003; parametrization by Maschberger 2013: $p_{Sys} \propto (m/\mu)^{-\alpha}(1+(m/\mu)^{1-\alpha})^{-\beta}$, $\alpha = 2.3$, $\beta = 2$, $\mu = 0.2$). $p_E$ resembles $p_{Sys}$ for stellar masses on, but underpopulates the brown dwarf region.

Not only the initial masses but also the time how long a star accretes can be distributed. This may have an important impact on the ensemble distribution: Basu & Jones (2004) found that an exponential distribution of growth times leads to a power-law tail instead of a lognormal distribution for linear growth with Gaussian fluctuations. Bate & Bonnell (2005) have similar findings for constant growth with lognormal fluctuations. The ensemble distribution function for both a distribution of initial masses and a distribution of growth

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Cumulative ensemble distribution functions for a population growing with a distribution of initial masses and an exponential distribution of growth times (corresponding to the middle panel of Fig. 7) stopping at $t = 0.1, 0.3, 0.6$ and $1.0$ Myr (from left to right).</figcaption></figure>

times is given by

$$p_E(m) = \int \int p(m_0)p(t)p_{iG}(m, t; \alpha, a, b, m_0) dm_0 dt. \quad (33)$$

The middle panel of Fig. 7 shows $p_E$ for a lognormal distribution of $m_0$, as above, and an exponential distribution of $t$ ($p(t) = \theta^{-1}\exp(-t/\theta)$). We choose $\theta = 0.3$, which is the mean growth time. At large masses the power law is flatter compared to the case without a distribution of growth times, populated by the stars that had more time to grow. $p_E$ is visibly shallower than $p_{Sys}$.

A uniform distribution of growth times (bottom panel, $t$ between 0 and 0.6 Myr) does not have the long tail like an exponential distribution, so that the power-law part of $p_E$ is steeper as it is not populated by the stars growing for a very long time. Here also the brown dwarf regime is more populated.

With a distribution of initial masses and growth times the ensemble mass function is also evolving in time. Fig. 8 shows the cumulative distribution of $p_E$ with an exponential distribution of growth times stopping at $t = 0.1, 0.3, 0.6$ and $1.0$ Myr (left to right curves). Rejection sampling has been used to obtain a sample $p_E$ containing 2000 variates. The mean masses are $\bar{m} = 0.11, 0.24, 0.39$ and $0.61 M_{\odot}$, and their time evolution not dissimilar to the evolution of the mean mass in hydrodynamical simulations of star formation (Bonnell, Vine & Bate 2004; Maschberger et al. 2010; Bate 2012; Krumholz, Klein & McKee 2012). The cumulative distribution functions seem to evolve somewhat more in time than in simulations (Krumholz et al. 2012; Bate 2012).

The parameters used in this section are chosen for illustrating purposes such that the resulting $p_E$ has some resemblance of the observed system mass function. Although the parameters have some reasonable values (for example, a star grows for a few hundred thousand years), it is very necessary to infer them from numerical simulations of star formation in order to constrain the scenario.

# 6 SUMMARY

We investigate the consequences of fluctuations in the accretion rates of non-linearly accreting stars by the means of a non-linear multiplicative stochastic process and find the following:

(i) Non-linear accretion, $\dot{m} \propto m^{\alpha}$, with fluctuations in the accretion rates leads to power-law tails in the distribution function of the final masses.

(ii) The main body of the distribution of final masses resembles a lognormal distribution.
(iii) Gaussian fluctuations produce a power-law tail $\propto m^{-\alpha}$ at high masses for $\alpha > 1$ and at low masses for $\alpha < 1$.

(iv) Lognormal fluctuations, approximated by the inverse Gaussian distribution, always produce a mass distribution function that has a power-law tail at high masses, even if accretion is linear ($\dot{m} \propto m$) or the $\alpha < 1$. Only for $\alpha \ge 1$, the exponent of the power-law tail is not evolving in time. For $\alpha < 1$, the power-law tail decays at very large masses.

(v) The shape of the mass distribution function depends on the initial mass, the average accretion rate $a$, the amount of fluctuations (ratio between standard deviation and average accretion rate $b/a$) and time.

(vi) The power-law tail is more and more populated in time, similar to time-evolution of the mass function obtained in numerical simulations of star formation. If observed at an early time the power-law tail may appear steeper because it is only sparsely sampled.

(vii) The final distribution function can resemble the whole shape of the initial mass function.

## ACKNOWLEDGEMENTS

I acknowledge funding via the ANR 2010 JCJC 0501 1 ‘DESC’ (Dynamical Evolution of Stellar Clusters). I would like to thank Cathie Clarke for very insightful discussions, as well as Estelle Moraux and Christophe Becker for their fruitful comments.

## REFERENCES

Basu S., Jones C. E., 2004, MNRAS, 347, L47

Basu S., Jones C. E., 2005, in Corbelli E., Palla F., Zinnecker H., eds, Astrophysics and Space Science Library, Vol. 327, The Initial Mass Function 50 Years Later. Springer, Dordrecht, p. 437

Bate M. R., 2012, MNRAS, 419, 3115

Bate M. R., Bonnell I. A., 2005, MNRAS, 356, 1201

Bate M. R., Bonnell I. A., Bromm V., 2003, MNRAS, 339, 577

Bonnell I. A., Bate M. R., 2006, MNRAS, 370, 488

Bonnell I. A., Bate M. R., Clarke C. J., Pringle J. E., 1997, MNRAS, 285, 201

Bonnell I. A., Bate M. R., Clarke C. J., Pringle J. E., 2001a, MNRAS, 323, 785

Bonnell I. A., Clarke C. J., Bate M. R., Pringle J. E., 2001b, MNRAS, 324, 573

Bonnell I. A., Vine S. G., Bate M. R., 2004, MNRAS, 349, 735

Cartwright A., Whitworth A. P., 2012, MNRAS, 423, 1018

Chabrier G., 2003, PASP, 115, 763

Chandrasekhar S., 1943, Rev. Mod. Phys., 15, 1

Dib S., Kim J., Shadmehri M., 2007, MNRAS, 381, L40

Dib S., Shadmehri M., Padoan P., Maheswar G., Ojha D. K., Khajenabi F., 2010, MNRAS, 405, 401

Edgar R., 2004, New Astron. Rev., 48, 843

Elmegreen B. G., Mathieu R. D., 1983, MNRAS, 203, 305

Filippov A. F., 1961, Theory Probab. Appl., 6, 275

Klessen R. S., Burkert A., 2000, ApJS, 128, 287

Kolmogorov A. N., 1941, Dokl. Akad. Nauk SSSR, 31, 99

Kolmogorov A. N., 1992, in Shiryayev A. N., ed., Selected Works of A. N. Kolmogorov Vol. 2: Probability Theory and Mathematical Statistics. Kluwer, Dordrecht, p. 281

Kroupa P., 2002, Sci, 295, 82

Krumholz M. R., McKee C. F., Klein R. I., 2005, Nat, 438, 332

Krumholz M. R., McKee C. F., Klein R. I., 2006, ApJ, 638, 369

Krumholz M. R., Klein R. I., McKee C. F., 2012, ApJ, 754, 71

Larson R. B., 1973, MNRAS, 161, 133

Larson R. B., 1978, MNRAS, 184, 69

Marcus A. H., 1968, Technical Report, NASA-CR-95938, Random In- dependent Splitting Model for the Mass Spectrum of Protostars and Interstellar Clouds. Bellcomm, Inc., Washington

Maschberger T., 2013, MNRAS, 429, 1725

Maschberger T., Clarke C. J., Bonnell I. A., Kroupa P., 2010, MNRAS, 404, 1061

Myers P. C., 2000, ApJ, 530, L119

Myers P. C., 2008, ApJ, 687, 340

Myers P. C., 2009, ApJ, 706, 1341

Myers P. C., 2011, ApJ, 743, 98

Myers P. C., 2012, ApJ, 752, 9

Nordlund Å. K., Padoan P., 1999, in Franco J., Carraminana A., eds, Interstellar Turbulence The Density PDFs of Supersonic Random Flows. Cambridge Univ. Press, Cambridge, p. 218

Øksendal B., 2002, Stochastic Differential Equations. Springer-Verlag, Berlin

Passot T., Vázquez-Semadeni E., 1998, Phys. Rev. E, 58, 4501

Thorin O., 1977, Scand. Actuarial J., 1977, 121

van Kampen N., 2006, Stochastic Processes in Physics and Chemistry. Elsevier, Amsterdam

Vazquez-Semadeni E., 1994, ApJ, 423, 681

Zinnecker H., 1982, Ann. New York Acad. Sci., 395, 226

Zinnecker H., 1984, MNRAS, 210, 43

This paper has been typeset from a TeX/Latex file prepared by the author.