On the grain-size distribution of turbulent dust growth

It has recently been shown that turbulence in the interstellar medium (ISM) can significantly accelerate the growth of dust grains by accretion of molecules, but the turbulent gas-density distribution also plays a crucial role in shaping the grain-size distribution. The growth velocity, i.e., the rate of change of the mean grain radius, is proportional to the local gas density if the growth species (molecules) are well-mixed in the gas. As a consequence, grain growth happens at vastly different rates in different locations, since the gas-density distribution of the ISM shows a considerable variance. Here, it is shown that grain-size distribution (GSD) rapidly becomes a reflection of the gas-density distribution, irrespective of the shape of the initial GSD. This result is obtained by modelling ISM turbulence as a Markov process, which in the special case of an Ornstein-Uhlenbeck process leads to a lognormal gas-density distribution, consistent with numerical simulations of isothermal compressible turbulence. This yields an approximately lognormal GSD; the sizes of dust grains in cold ISM clouds may thus not follow the commonly adopted power-law GSD with index -3.5, but corroborates the use of a log-nomral GSD for large grains, suggested by several studies. It is also concluded that the very wide range of gas densities obtained in the high Mach-number turbulence of molecular clouds must allow formation of a tail of very large grains reaching radii of several microns.


INTRODUCTION
Grain growth in cold molecular clouds (MCs) from seed grains present at the formation of MCs is a scenario which has been generally accepted for for a long time (Lindblad 1935;Baines & Williams 1965a,b). This type of dust formation is an important dust-formation channel in many models of various redshifts and galaxy types (see, e.g., Dwek 1998;Calura et al. 2008;Mattsson 2011;Valiante et al. 2011;Asano et al. 2013;Ginolfi et al. 2018) and depletion patterns in ISM gas are indeed consistent with dust depletion due to grain growth in MCs (see, e.g., Jenkins 2009;De Cia et al. 2016;Mattsson et al. 2019b).
Models of grain growth in the ISM usually rely on an assumption that the exact gas-density field can be replaced with the mean density, i.e., a kind of "mean-field approach" (as in the works of Asano et al. 2014;Hirashita & Aoyama 2019;Aoyama et al. 2020). Unfortunately, this "erases" smaller scale variations and other effects of dynamics. However, if the density variations are sufficiently small, this is a reasonable approach. But in case of strong compressible turbulence the gas density can vary by orders of magnitude and a significant fraction of the molecular gas in an MC display densities well above the critical density required for efficient grain growth (see Asano et al. 2013, for more details about this critical density).
In a homogeneous (constant density) environment, grain-growth by accretion is mainly limited by the abundance of the growthspecies molecules, which in turn is limited by the overall metallicity in the ISM. Thus, for a given metallicity, the gas density is decisive for the rate of grain growth. This is indicating that ⋆ E-mail: lars.mattsson@su.se modelling grain growth in the ISM in terms of a locally constant mean density may be an incorrect approach. The cold molecular phase of the ISM is highly inhomogeneous and display strong gas-density variations on sub-parsec scales. Such gas-density variations mean that some regions have number densities of growth species which are high enough to reach very fast grain growth. Moreover, since gas and dust tend to be coupled (at least on average), a majority of dust grains may actually reside in those regions. In a recent paper, Mattsson (2020, henceforth M20) showed that the overall rate of accretion can be increased by as much as two orders of magnitude when turbulent density variations are taken into account. Numerical simulations of interstellar turbulence (e.g., Klessen 2000;Price et al. 2011;Konstandin et al. 2012;Federrath 2013;Nolan et al. 2015) have demonstrated a direct relationship between the gas-density variance and the root-mean.square Mach number M rms , which gives that a relationship between M rms and the effective grain-growth velocity must exist. Hence, a significantly accelerated growth rate is expected in supersonic turbulence (see M20), which may resolve the apparent timescale crisis for dust formation at high redshifts (Mattsson & Höfner 2011;Rowlands et al. 2014;Watson et al. 2015).
In a recent set of simulations by Li & Mattsson (2020, submitted) it was seen that turbulence accelerated growth of dust grains can have also another effect: grains at different locations may experience different overall growth velocities. Consequently, the grain-size distribution may undergo more complex evolution than predicted by simple models (see, e.g., Hirashita & Kuo 2011;Mattsson 2016). The aim of the present paper is to show how density variations due to strong ISM turbulence can be a key factor in shaping the grain-size distribution.

Dust growth by accretion
In case dust growth by accretion of specific molecules takes place in a homogeneous medium, i.e., a gas of constant density, the mathematical description of this process is rather simple. The GSD f used here, can be seen as a probability density function if divided by total number density of grains n d , so that f = n −1 d (dn/da), and must satisfy a "continuity equation" of the form where a is grain radius and t is time. The thermal growth velocity ξ is the rate by which a increases due to thermal collisions (and chemical reactions) with the considered molecular growth species. It is easy to show that ξ is independent of a and given by (see, e.g. Hirashita & Kuo 2011;Mattsson et al. 2014), where X i is the mass fraction of the relevant growth-species molecules i in the gas, S is the sticking probability for a molecule hitting the grain andū t is the thermal mean speed of the molecules (which is assumed to be constant).
A linear transformation of the form where a 0 is the initial mean radius, gives the formal solution to equa- where f 0 is the GSD at t = 0. The choice of f 0 is, from a mathematical point of view, completely arbitrary, but any reasonable GSD must be skewed towards the small-grain end (Mattsson 2016).

Evolution of the GSD
The formal solution can mathematically be classified as a translational-invariant function, a mapping of the form f : R 2 → R, where in this particular case the plane formed by time t and grain radius a can be projected onto a line (in physics often referred to as a "travelling-wave solution"). In case the distribution f is not changing its shape one may say that it has a purely translational evolution. A linear transformation of the form η = a + b will in such a case merely shift f to the left if b > 0 and to the right if b < 0. If only accretion is considered, the translation must be to the right (b < 0). The left panel in Fig. 1 shows how an initially lognormal GSD, f 0 (a) ∝ exp{− 1 2 [ln(a) − ln(a) ] 2 /σ 2 a }, is shifted to the right as the grains grow, while the right panel shows the apparent steepening of the GSD on a log-scale. In the literature it is sometimes said that the GSD becomes steeper and skewed towards the small-grain end (see, e.g., Hirashita & Kuo 2011), which is only true in a relative sense. The relative increase of the radius ∆a/a of a small grain is much larger than ∆a/a for an initially large grain, which becomes apparent when plotted on a log-scale.
However, if ξ is also a function of a, the evolution of the GSD may not be translational on a linear scale any more. In particular, with ξ(a, t) = ξ 0 (t) a/a 0 , the evolution of the GSD can be described by the equation where ϕ(a, t) = a/a 0 f (a, t). This equation is solved by which is a solution that will appear translational-invariant on a logscale. An example of this is shown i Fig. 2, where a lognormal distribution is evolved assuming ξ ∝ a. In such a case, the GSD actually evolves towards being strongly dominated by large grains. This situation could effectively occur if ξ displays significant variance in space and time, in which case there is a certain probability that a grain will grow by a certain amount on an arbitrary time interval t + ∆t. Grains which happen to be in regions of high ξ (presumably with a high gas density ρ) will rapidly become large and may continue to grow faster than other grains which do not become large and continue to grow slowly because they reside in regions of low ξ (ρ, or ξ, in a fluid element can remain below or above the mean for an extended time, as seen in the time-series example in Fig. 3). Thus, the effective (integrated) ξ may be increasing with a.

Statistics of a turbulent gas
Interstellar gas is turbulent and highly compressible. Many numerical simulations as well as observational studies suggest rootmean-square Mach numbers M rms 10 in MCs (e.g., Brunt 2010; Price et al. 2011;Molina et al. 2012;Nolan et al. 2015), which means the turbulence in the cold ISM is strong, highly compressible and displays a wide range of gas densities. Numerical simulations of isothermal hydrodynamic turbulence with rotation-dominated forcing is known to produce roughly lognormal gas-density statistics (see, e.g., Federrath et al. 2010;Mattsson et al. 2019a, and references therein). Magnetohydrodynamic simulations also yield roughly lognormal statistics (see., e.g., Molina et al. 2012, and references therein), but with suppressed density variance for very high-M rms turbulence (Ostriker et al. 2001;Price et al. 2011). However, it is the lowdensity tail that tend to be suppressed (Molina et al. 2012), which means that the effect on processes mainly taking place in highdensity regions (e.g., dust growth by accretion of molecules) is small.
The lognormal distribution is of the form where µ is related to the variance/standard deviation such that mass conservation is obtained (Vazquez-Semadeni 1994; Konstandin et al. 2012). The variance can be estimated from its relation to the root-meansquare Mach number M rms , usually considered to be of the form which has been confirmed by several numerical experiments (e.g., Passot & Vázquez-Semadeni 1998;Federrath et al. 2010). A typical value for the case of for purely solenoidal forcing is b = 1/3. For mixed forcing, a value b ≈ 0.5 is often quoted (Federrath 2013). In the simulations described below, it will be assumed that b = 0.4 and the quoted M rms values are based on this assumption.

Turbulence as a Markov process
Numerical simulations of particle-laden strongly turbulent gas where the particles are interacting with each other or the gas are Figure 1. Evolution of an initially lognormal GSD assuming the growth velocity ξ is independent of grain size. The GSD is simply moving to the right (left panel ), which appear as "steepening" of the GSD when plotted on a log scale (right panel). The reason for this is that the relative growth of small grains is faster than the growth of larger grains.

Figure 2.
Evolution of an initially lognormal GSD assuming ξ ∝ a. In this case the GSD is simply moving to the right when plotted on a log scale, which means its variance is actually increasing.
computationally expensive and it can be difficult to identify the mechanisms behind emergent phenomena, such as the broadening of the GSD seen by Li & Mattsson (2020, submitted). It may thus be useful to try a more analytical approach to test the hypothesis outlined in previous sections. Here, the density variations due to turbulence will therefore be modelled as a stochastic process and not by solving the equations of fluid dynamics. It is well-established that the velocity field of a turbulent fluid can be described as a Markov process (see, e.g., Novikov 1989;Pedrizzetti & Novikov 1994) and the gas-density PDF in highly compressible turbulence has recently been modelled in a similar way (Mocz & Burkhart 2019;Scannapieco & Safarzadeh 2018). In the latter case, the approach is to describe turbulence in terms of the statistical time-evolution of fluid-element densities. The gasdensity PDF can be regarded as made up by volume or mass elements and turbulence as a temporally homogeneous Markov pro-cess in which the future state space only depends on its current values. If the logarithmic density parameter s (see section 2.2) is an Ornstein-Uhlenbeck (OU) process (Uhlenbeck & Ornstein 1930), a special type of Markov process to be described below, the resultant gas-density PDF P(ρ) is a lognormal distribution.
A Markov process R of OU type, essentially a mean-reverting random walk, will obey a Langevin (1908) equation of the form where τ is the relaxation time, D = 2σ 2 s /τ is a diffusion constant and N(µ, σ 2 ) denotes a temporally uncorrelated Gaussian random variable with mean µ and standard deviation σ. The above equation can be equivalently written as a stochastic differential equation with a noise term, where γ(t) is Gaussian white noise with µ = 0 and σ = 1/dt in the present case. The equivalence of the two equations (9) and (10) is a simple consequence of the property N(α + β µ, β 2 σ 2 ) = α + β N(µ, σ 2 ) of a Gaussian random variable. The Markov process described above has a steady-state solution for the PDF of R. The PDF P(R, t) obeys the Fokker-Planck equation, which for a steady state (∂P/∂t = 0) reduces to a simple ordinary differential equation, to which the only nontrivial solution is a Gaussian distribution. The time it takes for an OU process to reach this steady-state distribution is about twice 1 the relaxation time τ. The timescale (unit) τ is the relaxation time of the process and the density ρ is renormalised such that ρ is the mean also when taken over the time series.

Numerical implementation of an Ornstein-Uhlenbeck process
The first term on the right-hand side of the Langevin equation (9) is a deterministic term (in dt), while the second term is a stochastic term. The infinitesimal step is a Gaussian random variable and the derivative is white noise, as described above. Numerical treatment of an OU process is commonly based on the Euler-Maruyama method, which involves discretising time and adding infinitesimal steps to the process at every time step δt (see Kloeden & Platen 2011). That is, for the Langevin equation above, the scheme is where all parameters are as previously defined. N(0, 1) is a random Gaussian variable with mean µ = 0 and variance σ 2 = 1, which is independent at each time step. The factor (δt) 1/2 arises from the fact that the incremental step for white noise must have σ = (δt) 1/2 . Consequently, the error of the Euler-Maruyama method is of order (δt) 1/2 . The scheme outlined above is a first order scheme and it should be noted that higher-order schemes can also be used (Gillespie 1996).

Modelling a dust-gas system
In order to simulate a dust-gas system a few of simplifying assumptions have to be made. Dust grains are inertial particles, which means they can decouple from the carrier fluid (the interstellar gas) if the frictional drag force acting on the grains is low enough relative to the inertial forces of the grains. As explained in M20, dust-gas coupling in turbulence is a both physically and mathematically complex problem, which requires advanced numerical simulation. To make the problem (grain growth by accretion) analytically tractable it is timescale (Karatzas et al. 1991). After 2.3 e-folding times, R has moved 90% of the distance between its initial and final values, regardless of what those values are, which is here considered to be sufficiently close to the (final) steady state to say the system is in the steady-state phase. more or less necessary to assume that dust and gas are position coupled on the spatial scale of interest. Thus, a grain residing in a given fluid element at t = 0 will reside in the same fluid element at any later time. Another assumption to be made is that the initial mass fraction of the growth-species molecules in the gas is a "universal" constant, i.e., the same for every fluid element. Effects of magnetic fields are here considered unimportant, but it should be noted in passing that the Lorentz force acting upon the grains may lead to more efficient coupling between gas and dust. Moreover one may also assume that the thermal mean speed of the gas particles/molecules is everywhere the same; the modelled system is assumed to be strictly isothermal. This assumption means that the growth velocity ξ will only vary due to variations in the gas density ρ.
In the present study, following the discussion above, it is the growth velocity ξ that is the stochastic process of interest. As mentioned above in section 2.1, one may assume ξ = ξ 0 ρ/ ρ , with ξ 0 a constant, for as long as depletion of the growth species is not significant. Thus, the stochastic variation of s = ln(ρ/ ρ ) and ln(ξ/ξ 0 ) is actually governed by the same OU process, albeit with the addition of an "amplitude correction" due to depletion of the growth species. Discretising time in equation (1) and introducing s as a discrete OU process in that equation and tracking a large number N e "fluid elements", one obtains a system of the form where f n is the discretised/binned GSD after n time steps and ξ 0 is the growth velocity at the initial time step (t = 0). The GSD evolution is obtained by binning the elements at all (or selected) time steps by applying the histogram function HIST[ℓ, ∆ log a](. . .), where ℓ is the number of bins and ∆ log a is the logarithmic bin size. N e must be large to ensure proper statistics and generally N e ≫ N and N ≫ ℓ.
To include the effect of depletion one can add a depletion factor F ⋆ to the system of equations above. More exactly, the growth velocity-equation should then be replaced with where ρ d is the local dust mass density, s d = ln(ρ d / ρ ) and s i = s ln[X i (0)], following the model of depletion in Mattsson (2016).

Numerical simulation of GSD evolution
By numerical solution of the system of equations (14) using the Euler-Maruyama method, the evolution of large number N e of "fluid volume elements" can be followed and the evolution of the GSD can be reconstructed by binning the sizes of the grains associated with each element. In the simulations presented here N e = 10 5 and the initial sizes of the grains associated with each fluid element is either the same for all elements ("delta-distributed" case) or randomly drawn from a power-law distribution with a slope −3.5 ("MRN" case, see Mathis et al. 1977). Each time series is computed with a constant time step δt/τ = 0.01 and the length of the time series is 20 relaxation times, i.e., t max = 20 τ. That is, the density evolution of each one of the 10 5 fluid elements is modelled by a random walk with 2000 steps. Table 1. Parameters of the OU simulations used in the present paper. All parameters are dimensionless except the grain radii, a 0 , a min and a max , where the latter two are the initial minimum and maximum radii of the grains. To explore what role the strength of the turbulence plays, simulations has been performed with different σ s (see Table 1) corresponding to a wide range of M rms values. For convenience, the unit time is equal to the relaxation time of the OU process, i.e., all simulations are made with τ = 1 (arbitrary system of units). The lognormal gas-density PDF is assumed to be mass weighted, which requires that s = 1 2 σ 2 s to ensure mass conservation (Vazquez-Semadeni 1994;Scannapieco & Safarzadeh 2018). The parameters governing the OU process are dimensionless and the initial/average growth velocity ξ 0 is the only parameter that requires physical scaling in order to obtain grain sizes in physical units. In units of τ, the adopted value in all simulations presented here is ξ 0 τ = 10 −3 µm. Following equation (2), an estimate of ξ 0 can be made from where S is the sticking probability, u therm and X i (0) are the mean thermal velocity of the relevant molecules and their fraction of the total gas mass, respectively, and ρ gr is the bulk material density of the grains. Assuming a mean density of ρ ∼ 10 −21 g cm −3 , u therm ∼ 0.1 km s −1 , X i (0) ∼ 0.01, ρ gr ≈ 3 g cm −3 and S ≈ 0.3, the initial growth velocity is of the order ξ 0 ∼ 10 −25 cm s −1 . The assumption ξ 0 τ = 10 −3 µm then implies τ ∼ 1 Myr in physical units. Fig. 4 shows how the GSD evolves in the different cases considered. The general trend appears to be that a power-law tail develops at early times regardless of the initial condition and after some transitional phase of temporary flattening, the evolution progresses towards a lognormal shape and a log-translational phase (see Section 2.1.1 and Fig. 2). It should be noted that larger values of σ s (higher average Mach number M rms ) seem to create a deviation from a lognormal tail for large a in the GSDs even for the case of a monodispersed (delta-distributed) initial GSD. This phenomenon could in part be due to low-number statistics in the tails of the GSD, but a longer run with N e = 10 6 and t max = 40 (Run 13, see Table 1) still shows a power-law tail, which suggests that it could be either a real effect or an effect of the fact that the numerical method is a first-order scheme.
In general, however, there is a consistent irrefutable trend: strong turbulence (high σ s values) clearly leads to a wider GSD and generally larger grains. The average grain radius a at t = 20 is approximately 40 times larger for σ s = 2 compared to σ s = 0.25. To first order log( a ∝ σ 2 s at any given time (see Fig. 7), which can be understood in terms of theoretical predictions in M20, in particular equation (17) in that paper. Note that the simulations deviate from the M20 theory at early times, which may be explained by the fact that the number of fluid elements is finite and not large enough to ensure that the simulated gas-density variation strictly follows a lognormal distribution at any given time (which is the case for the analytical theory of course).
To show the evolution towards a lognormal GSD more clearly, all simulations displayed in Fig. 4 are made without any correction for depletion. In reality, depletion is negligible only at early times or if the initial dust-to-gas ratio is very low. The effect of including depletion (see equation 15) of growth species (Run 11 and 12, see Table 1) is shown in Fig. 5, where one can see that, assuming an initial depletion of 10%, the evolution towards a lognormal GSD is slower and may eventually stop before the GSD has obtained a truly lognormal shape. That is, the GSDs with power-law tails seen at intermediate stages can be the end states of the evolution and the log-translational phase may not occur. It is also noteworthy, that if all grain growth is depletion limited, i.e., if the key growth species is fully depleted within the lifetime of an MC, then a will eventually be the same regardless of what σ s is, provided that all other parameters are unchanged. A central point in the M20 theory is that a state of total depletion will likely not occur within the lifetime of an MC unless the growth is accelerated by turbulence.
As dust extinction as well as emission is strongly correlated with dust mass, it is of interest to consider the dust-mass density ρ d and how it evolves. More precisely, the quantity a dρ d /da is displaying how the total dust mass evolves and how the dust mass is distributed over the various grain sizes a. Fig. 6 shows this grain-mass distribution (GMD) with and without depletion for both a mono-dispersed and a MRN-like initial GSD. From Fig. 6 it is evident that the total amount of dust that can be formed is limited by depletion, which is a rather trivial result, while the overall shape of the GMD is very similar, which is a more interesting result. The GMD is initially peaking at small a in case of a mono-dispersed GSD, and at large a if the GSD is initially MRN-like, but in both cases the GMD evolves towards having a peak at essentially the same a. Obviously, depletion limited growth leads to a peak at somewhat smaller a; roughly a factor of two smaller than without depletion, in the present case.
As indicated above, in section 2.1.1, grains residing in high-density regions will grow fast and become large and as they are likely to remain in a high-density region for a while and, similarly, grains  . Examples of GSD evolution including the effect of depletion, assuming an initial depletion of 10%. The evolution towards a lognormal GSD is slower and may eventually stop before the GSD has obtained a truly lognormal shape. The power-law tail seen at intermediate stages may therefore remain. located in low-density regions will be small. Thus, the growth velocity ξ will effectively depend on the grain radius a. More precisely, ξ would increase with a. For simplicity, one may assume that ξ(a, t) = ξ 0 (t) (a/a 0 ) 1+β , where β 0. Then, by the generalised formal solution in Appendix A, where C = n d /( √ 2 π a 0 ), σ η is the standard deviation of η = a 0 /n (a/a 0 ) −β − A(t) and For the special case β = 0 the solution is a lognormal distribution instead. For β > 0 the solution will have a power-law tail and look very similar to GSDs seen at early times in Fig. 4. Thus, one may say that, effectively, the evolution of the GSD in turbulence corresponds to having a grain-size dependent ξ, where the dependence on the grain radius a is very steep at early times and evolve towards a linear relation and a lognormal GSD. The argument above is indeed very sketchy and phenomenological, but serves to prove that the usual assumption, that growth of spherical grains can be described by a ξ which is independent of a, can be called into question if dynamics is not included in the model. In "zero-dimensional" models, without turbulent gas dynamics and resultant density variation, grain-growth by accretion of molecules will only lead to translational evolution of the GSD (see section 2.1.1). But such models can be modified parametrically to include the effects of turbulence by adding the M rms correction to the overall growth velocity suggested in M20 in combination with some dependence on a, as described above. Developing such a parametric modification in detail goes beyond the scope of the present paper, though.

Scaling with Mach number?
An obvious question to ask in connection to parametric models is whether there exists some kind of scaling with the Mach number M rms . Based on Fig. 4 it is tempting to suggest a relation between the standard deviation of the GSD σ a and M rms similar to the wellestablished relation between σ s and M rms (see equation 8). However, closer scrutiny shows that σ a and σ s may not have a simple functional relation. First, the initial GSD must play a role for the resultant σ a , which can be seen in the upper panels of Fig. 4, showing the two simulations with σ s = 0.25. Second, since the evolution of the GSD is depletion limited, the initial level of depletion is also important as it determines how much grains will be able to grow, which in turn affects σ a . Third, as mentioned above in Section 3.1, grains may grow large enough to decouple from the flow, which can affect the over all rate of growth. In this case the scaling of the problem, e.g., mean gas density ρ and σ s , will be crucial.
Given the dependencies on the initial conditions listed above, it is actually unlikely there exists a universal relation between the width/variance of the GSD and M rms . But the qualitative result, that the GSD becomes broader due to turbulence, is quite clear and the effect must be larger at high M rms , regardless of initial conditions.

Other forms of dust processing in turbulence
The present study has focused on grain growth by accretion of molecules, but there are other forms of dust processing that may become important in turbulent MCs. In particular, growth by coagulation/aggregation is important in regions of very high density; the grain-grain interaction rate in compressible turbulence increases mainly due to over-density effects rather than turbulent velocities (see Li & Mattsson 2020). Grain-grain interactions at sufficiently high energies may also lead to shattering (see, e.g., Slavin et al. 2004;Hirashita & Yan 2009). But if coagulation dominates over shattering, the evolution of the GSD is in fact not so different from the evolution seen here due to accretion of molecules, except that small grains will always remain, creating a stretched-out GSD with a lognormal-type slope at the large-grain end (see Li & Mattsson 2020 , Fig. 2). The reason for this similarity between the large-grain tails is likely that while the grain growth is driven by over-densities of molecules in one case, it is driven by over-densities of small grains accreting onto large grains in the other. Small grains trace the gas quite well and therefore they follow the gas PDF just like molecular growth species. Figure 6. Evolution of the grain-mass distribution dρ d /da as a function of grain radius a. Upper panels show the case without depletion, while the lower panels show simulations with depletion. Left panels show evolution from a mono-dispersed initial GSD, while the right panels show evolution from an MRN-like initial GSD.

SUMMARY AND CONCLUSIONS
In a previous study (M20) it was shown that turbulence can significantly accelerate dust growth by accretion of molecules onto grains, where the growth rate scales with the square of the Mach number. Here, it is has been shown, by simulating isothermal turbulence as an Ornstein-Uhlenbeck process, how turbulence must have a significant impact also on the resultant GSD, as seen in recent hydrodynamic simulations (Li & Mattsson 2020, submitted). In particular, the variance ("width") of the GSD increases with the mean Mach number M rms , although a generic scaling relation may not exist.
The turbulence-induced broadening of the GSD implies that a fraction of very large grains, with radii orders of magnitudes larger than the initial mean radius, can form without any need to assume extreme conditions. For M rms ∼ 10 a significant fraction of micronsized grains could form by growth by accretion only.
The shape of the GSD at later stages of evolution appears to be a reflection of of the gas-density PDF. That is, a lognormal distribution of gas densities tend to eventually produce a lognormal GSD. This corroborate the use of a lognormal GSD for large grains in ISM dust models (Jones et al. 2013). Of course, the initial GSD plays a role, but the "memory" of the initial shape of the GSD gradually fades in all simulations presented here, while the total number density of grains is conserved. In case of highly compressible turbulence (high M rms ), the simulations seem to predict slightly skewed GSDs with a large-grain excess compared to simulations corresponding to lower M rms , but this may be a statistical artefact.
Modelling turbulence as a stochastic process has obvious limitations, which is why the results presented above need to be confirmed with numerical simulations where one is actually solving the equations of fluid dynamics. Furthermore, the theory by Baines et al. (1965) predicts that "drift", i.e., dynamical decoupling of dust and gas, may have an important effect on the growth rate, which is most easily explored with detailed numerical simulations. Such efforts are currently under way (Li & Mattsson 2020, submitted).

ACKNOWLEDGMENTS
The author wishes to thank the anonymous reviewer, whose comments, suggestions and criticism was much appreciated. This work is supported by the Swedish Research Council (Vetenskapsrådet), grant no. 2015-04505.

DATA AVAILABILITY
The code and corresponding output data underlying this article will be shared on reasonable request to the corresponding author.