Dust growth by accretion of molecules in supersonic interstellar turbulence

We show that the growth rate of dust grains in cold molecular clouds is enhanced by the high degree of compressibility of a turbulent, dilute gas. By means of high resolution (10243) numerical simulations, we confirm the theory that the spatial mean growth rate is proportional to the gas-density variance. This also results in broadening of the grain-size distribution (GSD) due to turbulence-induced variation of the grain-growth rate. We show, for the first time in a detailed numerical simulation of hydrodynamic turbulence, that the GSD evolves towards a shape which is a reflection of the gas-density distribution, regardless of the initial distribution. That is, in case of isothermal, rotationally forced turbulence, the GSD tends to be a lognormal distribution. We also show that in hypersonic turbulence, decoupling of gas and dust becomes important and that this leads to an even further accelerated grain growth.


INTRODUCTION
Dust plays an important role in the formation of stars and planets, but also the evolution of galaxies throughout the whole history of the Universe. Moreover, Bernstein et al. (2002) estimated that dust re-radiates about 30% of starlight in the Universe. Therefore, dust regulates emission and absorption of stellar light.
Even at redshifts as high as z ≈ 4 − 8 galaxies display high dust fractions, which requires a very fast dust-formation process and a sudden appearance of dust as soon as there is enough metals to form it (Rowlands et al. 2014;Mattsson 2015Mattsson , 2016Watson et al. 2015). Dust grains are formed mainly by condensation of metals into solid-state material, e.g., silicates, graphite and amorphous carbon (see Mathis 1990;Draine 2003, and references therein). But accretion of volatiles, forming "ice mantles", also accounts for some of the dust mass (Jones et al. 1994).
Type II Supernovae (SNe) and asymptotic giant branch (AGB) stars not only supply the metals in the cosmic matter cycle, but also a significant amount of dust (Valiante et al. 2009). However, the destruction of dust in the interstellar medium (ISM) is significant and therefore a replenishment mechanism is needed (Mattsson 2011;Valiante et al. 2011;Rowlands et al. 2014). Growth of dust grains by depleting metals (molecules) in the ISM is the most efficient in cold Corresponding author: Xiang-Yu Li xiangyu.li@pnnl.gov molecular clouds (MCs), which is indicated by absorptionline observations (see, e.g., Savage & Sembach 1996;Jenkins 2009;De Cia et al. 2016)  showed that the evolutionary history of the GSD is important for understanding the dust abundance in galaxies. A key question in this context is what fraction of the dust mass is in large (micron sized) grains. Dust grains grow by accretion of molecules in MCs and by coagulation of interacting grains (Li & Mattsson 2020). Thus, growth of dust grains in the ISM is an essential process for the initial phase of planet formation, which is not yet fully understood (Johansen et al. 2007).
MCs in the ISM are so dilute that the typical size of dust grains is much smaller than the mean free path of the gas molecules. Baines & Williams (1965) derived an expression for the accretion, concluding that the mass growth rate of dust grains is proportional to, and primarily regulated by, the surface area of the grains. This was done under the assumptions that dust grains are spherical, a fixed proportion of the colliding gas molecules adhered to the grain and that the gas density is homogeneous. However, they also showed that the relative velocity of grains moving in the gas plays an important role, which is interesting in a context of high Mach number turbulence.
Interstellar turbulence displays high Reynolds numbers and is highly compressible with typical Mach numbers larger than 10 (Nolan et al. 2015). This results in vigorous variations of the gas density, which could potentially change the mass growth rate. Mattsson (2020) included gas-density variation in the accretion process and found that the spatial mean growth velocity of the dust grain radius scales with the square of the Mach number.
Growth and destruction of dust grains shape the GSD, which in turn determines the extinction curve of galaxies and is an important parameter in planet and stars formation (Relaño et al. 2020). Mathis et al. (1977) obtained a power law distribution of interstellar dust grains with an exponent of about -3.5 by fitting observed interstellar extinction in the diffuse ISM. This power law distribution has been widely used and has become the canonical GSD, known as the "MRN" distribution. However, the typical shape of the GSD in MCs is not known.
In this study, we will validate the theory of Mattsson (2020, henceforth "M20") with state-of-the-art high resolution numerical simulations of hydrodynamic turbulence including the accretion process of molecules onto dust grains. We will show below that the M20 theory agrees perfectly well with the simulation results. Based on M20 theory, we will demonstrate that this grain growth facilitated by gas density variations leads to a lognormal GSD that is independent of the initial shape of GSD.

Momentum equation of interstellar turbulence
Interstellar turbulence is governed by the Navier-Stokes equation: where f is a forcing function (see Brandenburg 2001), p is the gas pressure, and ρ is the gas density that follows the continuity equation, The viscosity term F ν visc is given by where S = 1 2 ∇u + (∇u) T − 1 3 (∇ · u) I is the rate-of-strain tensor with I the unit tensor and ν is the kinetic viscosity of the gas. However, for high Mach-number simulations, an artificial shock viscosity is needed to smooth the shock such that the shocks can be captured without causing a singularity in the velocity field. With this shock viscosity included, F visc becomes a superposition of F ν visc and F shock visc , where the shock viscosity ζ shock is given by ζ shock = c shock max[(−∇ · u) + ] (min(δx, δy, δz)) 2 .
Here c shock is a constant defining the amplitude of the shock viscosity (Haugen et al. 2004), δx, δy, and δz, are the lengths of the sides of a mesh cell.
The stochastic solenoidal forcing f is given by where k(t) is the wave vector, x is position, and φ(t) (|φ| < π) is a random phase. The normalisation factor is given by N = f 0 c s (kc s /∆t) 1/2 , where f 0 is a dimensionless factor, k = |k|, and ∆t is the integration time step (Brandenburg & Dobler 2002). In the present study, we choose a completely nonhelical forcing, i.e., where e is the unit vector.
The two dimensionless parameters that characterizes compressible turbulence are the Reynolds number Re and the root-mean-square (rms) Mach number M rms . Re is defined where u rms is the rms turbulent velocity and k f is the forcing wave number. M rms is defined as where c s is the sound speed. We assume that the gas is isothermal, so that c 2 s = γp/ρ = constant, where γ = c P /c v = 1 with the specific heats being c P and c V at constant pressure and constant volume, respectively.

Particle dynamics
The momentum equation for dust grains is given by and The stopping time τ i due to the kinetic drag is given by where a is the radius of dust grains and ρ gr is the material density of dust grains. Eq. (12) is derived under the assumption that the particle radius is much smaller than the mean-freepath λ, which is the case in most astrophysical contexts (large Knudsen number, Kn = λ/a 1 Armitage 2010; Mattsson et al. 2019). In the limit of low relative Mach numbers (W = |u − v i |/c s 1), Eq. (12) reduces to When W 1, the correction term in the parenthesis of eq. (12) is dominating and the momentum equation for dust becomes nonlinear (see Schaaf 1963;Kwok 1975;Draine & Salpeter 1979).

Accretion process
The rate of change of the grain mass due to accretion of molecules is (Mattsson 2020) where X i ∼ 10 −3 is the mass fraction of the relevant growthspecies molecules i in the gas, S is the sticking probability for a molecule hitting the grain and is the thermal mean speed of the molecules. Here, k B is the Boltzman constant and m is the molar mass of the molecular gas. Althoughū t is dependent on the composition of molecular species, we here focus on how supersonic turbulence affects the accretion process for a given generic composition. That is, we assume a constantū t for convenience. As Eq. (15) can be written asū t = √ 8/πc s , we shall assume c s is constant as well, i.e., isothermal conditions. Since m gr = 4 3 πa 3 ρ gr , Eq. (14) can be written in terms of a as We define the growth velocity as ξ = S X iūt ρ ρ gr . Thus, Eq. (16) can be written as where ξ(t) → 0 as t → ∞ due to depletion of growth-species molecules. In the following we will assume ξ = constant, however. There are two reasons for this assumption: (1) a constant ξ saves a lot of computation time; (2) we will have a better chance of reaching a regime where effects on the growth rate due to decoupling of gas and dust (Baines & Williams 1965) can be seen if the growth is not depletion limited. Assuming ξ = constant may not be entirely realistic, but it is reasonable in an early phase of grain growth and could also be justified by the possibility that undepleted gas can be mixed into the modelled region at a rate similar to that of depletion.

Statistical description
According to Eq. (16), with ξ = constant, the change of growth rate is determined by ρ(x, t) only. Using a spatialmean approach, M20 proposed that the mean growth of dust grains due to accretion is given by considering periodic boundary conditions. Furthermore, the standard deviation σ s of the logarithmic density parameter s = ln(ρ/ ρ ) is related to the standard deviation σ ρ of the linear gas density, which in turn is empirically related (via simulations) to the mean Mach number M rms , although the exact relation can be debated (see, e.g., Passot & Vázquez-Semadeni 1998;Lemaster & Stone 2008;Price et al. 2011;Federrath et al. 2010). The hypersonic turbulence of cold MCs is therefore indicative of very rapid grain growth according to the M20 theory.
Depletion is omitted in the present study for mainly two reasons. First, M20 does not consider the effect of dustgas drift as the grains become large as predicted by Baines & Williams (1965). In order to test the validity of neglecting this drift effect, it would be ideal if the growth is not depletion-limited since the drift effect can easily be masked by depletion. Second, simulations with depletion are computationally very demanding. Adding depletion in a physically consistent way is straightforward, but results in a drastic increase of communication between processors. The increase of the computing time is significant and we would not be able to run a whole suite of simulations with our present resources.

Initial parameters
We have designed our simulations such that they can match observed dust grain properties and interstellar turbulence, save for the very high Re of ISM. For the sake of convenience, we define a unit length L 0 = 3.086 × 10 18 cm, unit velocity v 0 = 10 5 cm s −1 , and unit mass density ρ 0 = 6.771 × 10 −23 g cm −3 . This results in an unit time of t 0 = 3.086 × 10 13 s ≈ 10 6 years. Unless anything else is stated we adopt the Gaussian system of units (a.k.a. cgs units) in this study and introduce simulation units defined within this system.
To reproduce the properties of a cold MC we assume an initial mean gas density ρ ini = 49.33ρ 0 and a sound speed is about c s = 0.3v 0 . Since the gas (molecular hydrogen) is assumed to obey the ideal gas law, the corresponding temperature is about 30 K. For the dust we employ a mono-dispersed and a power-law (MRN Mathis et al. 1977) initial GSD with a mean initial radius of dust grains a ini = 10 −23 L 0 . The power index is -3.5 such that f (a, 0) is in the form of f (a, 0) = a −3.5 with cutoff a min = a ini and a max = 3a ini . The material density of dust grains is ρ gr = 3.545×10 22 ρ 0 = 2.4 g cm −3 . We adopt the MRN distribution for the following reasons. First, even though the MRN distribution is derived from diffuse ISM instead of dense MCs in the cold ISM, the initial GSD of an MCs could in fact be rather similar to that of the diffuse ISM because dense MCs are formed from diffuse ISM. Second, the MRN distribution is often used as the "canonical GSD" in ISM and is considered the natural assumption when no actual constraints exist.

Numerical method
We use the Pencil Code to solve equations described above, which is a finite difference code designed, in particular, to solve compressible turbulence. A third order Runge-Kutta time stepping is adopted. There are 1953120 individual dust grains (Lagrangian inertial particles) being tracked in the simulations and the growth of dust grains is monitored at each time step. Tracking such a large number of dust grain is numerically demanding, but ensures sufficient statistics.
We simulate a wide range of M rms (1.38-11) covering transonic to hypersonic turbulence with similar Reynolds numbers. The simulation with M rms ≈ 11 reaches the typical values of M rms observed in large-scale interstellar turbulence (Brunt 2010;Elmegreen & Scalo 2004). We run for 9 simulation time units in the case of M rms ≈ 11, using 1024 CPUs and a wall-clock time of 24 × 4 = 96 hours. There are 61360 time steps integrated with a time step dt = 5.3 × 10 −5 , which corresponds to about 21 turnover times. Accretion is turned on when turbulence is well-developed. A summary of all simulations is presented in Table 1. , and the Eulerian-mapped particle distribution n(x, t) (lower panel) of four simulations with different M rms . As M rms increases, M(x, t) exhibits more intricate structure and a more inhomogeneous distribution. The same is observed for ln[ρ(x, t)], where more pronounced gas filaments are seen at high M rms . Larger M rms also results in stronger spatial clustering of dust grains as shown in the lower panel of Fig. 1. 3.2. Comparison with the M20 theory M20 assumed that dust grains behave like tracers but can grow due to accretion as described by Eq. (18). Hence, we started out with a test simulation (corresponding to Run B) where dust grains are treated as tracers. Fig. 2(a) shows that the simulation result (solid cyan line) agrees more or less perfectly with the theory (red dashed line) in the supersonic regime (M rms = 3.26). The only assumption involved in Eq. (18) is the periodic boundary conditions, which is required to apply the spatial-mean approach on Eq. (16). This agreement confirms that the mean growth rate of dust grains is proportional to the square of mean variance of the gas density in the supersonic regime as expressed in Eq. (18), provided that dust grains are treated as tracers.
Dust grains are inertial particles and their inertia increase as they grow. Therefore, we have focused on simulations where Eq. (10)-Eq. (12) is solved at different M rms . Fig. 2(a) shows the evolution of the mean radius a in these simulations. As expected, the spatial mean growth rate of dust grains increases with increasing M rms . In the transonic regime, the simulation result (solid magenta line) agrees very well with Eq. (18). In the supersonic regime and beyond, the simulation results agree with the theory in the beginning, but deviates at later times. This indicates that Eq. (18) does not hold in the supersonic regime and beyond, which was not the case in the tracer-particle limit. As we can see from Eq. (18), the mean growth rate is determined by the relative variance ρ 2 / ρ 2 , which is shown in Fig. 2(b). Since ρ 2 / ρ 2 becomes stationary, the evolution of a as the time integral of d a /dt, is approximately linear with time as shown in Fig. 2(a), as predicted by Eq. (18) if depletion is negligible.
To understand the deviation from the M20 theory, we shall recall that the theory assumed that dust grains are tracers. However, dust grains have inertia, and the kinetics of such particles will deviate from that of tracers. As shown in Fig. 3, u rms and v p, rms are very similar in the beginning, meaning that both gas and dust must be accelerated to a certain rms velocity before decoupling, at a statistical level, can actually occur. At later times, v p, rms becomes clearly smaller than u rms , indicating that dust grains decouple from the flow even on average. The decoupling effect is stronger for large M rms , which is expected. Therefore, the deviation from the theory seen in Fig. 2 is a result of dynamic decoupling of dust grains from the turbulent flow. Fig. 4 shows GSDs at fixed times (snapshots), where higher M rms results in a broader GSD. This is because higher M rms obviously leads to stronger variations of ρ as shown in Fig. 1 and therefore faster growth of dust grains in highdensity regions. It is noteworthy that the dust grain population was mono-dispersed initially.
The evolution of the width/variance of the GSD in terms of σ a is shown in Fig. 5. First, we note that σ a grows with time, which is due to the variations of ρ. Second, σ a increases with increasing M rms , which demonstrates that turbulenceinduced growth by accretion broadens the GSD. Note, however, that assuming ρ = constant (M rms = 0, no turbulence), leads to σ a = 0 at all times for an initially mono-dispersed Table 1. Re = 100 is kept the same for all the simulations by changing f 0 and ν. c shock = 2. L x = 5 length unit. The number of particles is N p = 1953120. "cgs" unit is adopted. 3.3. Shape of the GSD: a reflection of the gas density PDF?
As we mentioned in 3.2, f (a, t) rapidly develops into a lognormal distribution from a mono-dispersed initial distribution. Mathis et al. (1977) derived a power-law GSD based on observational constraints in diffuse ISM, which is now the canonical, widely used GSD known as the "MRN" distribution. Fig. 6 shows the evolution of f (a, t) starting from a power-law GSD with an MRN slope. A lognormal distribution is the result also in this case, which demonstrates that  the lognormal GSD resulting from grain growth in isothermal, solenoidally forced turbulence is essentially independent of the initial GSD. To understand the mechanism behind this phenomenon, we shall examine the PDF of the gas density. The gas-density PDF in our simulations is a lognormal distribution. Previous studies (e.g., Federrath et al. 2009;Hopkins & Lee 2016;Mattsson et al. 2018) of isothermal hydrodynamic turbulence have also generated lognormal distributions of gas density, which can be expressed as  where s = ± 1 2 σ 2 s depending on whether the PDF is mass or volume weighted (Li et al. 2003). The GSD is determined by summing up the local evolution of da/dt in every mesh cell or, which is equivalent, by following the evolution a large number of fluid elements. The local growth velocity da/dt is proportional to the local density of gas ρ (x, t) as in Eq. (16). Because ρ follows a lognormal distribution, i.e., Eq. (20), it is very likely that f (a, t) develops into a lognormal distribution. From a more general point of view, it seems that the resultant GSD is a reflection of the gas-density PDF, whatever form this PDF may have. The mathematical theory behind this is described by Mattsson (2020, submitted).

DISCUSSION
By applying a spatial mean approach to Eq. (16), M20 predicted that the mean growth rate of dust grains by accretion is proportional to the mean of the second moment of gas density as described by Eq. (18). Simulation results in the present study show that this scaling holds either at early times, when grains are small and their inertia is insignificant, or in the case of small M rms . The analytical theory deviates slightly from the numerical simulations at high M rms , which we attribute to a higher degree of dynamical decoupling between gas and dust at higher kinetic energies. We also see that ρ 2 / ρ 2 becomes stationary and therefore a evolves linearly with time.
Assuming a lognormal gas-density distribution, M20 further proposed that the mean growth rate scales with M rms . However, the semi-empirical nature of the σ s -M rms relation combined with the fact that our simulations have a finite and not very high Re, have led us to think it is generally better to compare the M20 theory with our simulations in terms of ρ 2 (or σ ρ ) directly, rather than M rms . The theory predicts an explicit relation with σ ρ , while the connection with M rms is implicit.
We find/confirm that the high degree of compressibility of ISM turbulence can affect grain growth by accretion of molecules in two ways. On one hand, it increases the mean growth rate of dust grains, i.e., d a /dt ∝ ρ 2 . On the other hand, it broadens the GSD as the strong variations of ρ(x, t) leads to a similar variation of the growth rate. The enhancement is about a factor of 10, going from the transonic to the hypersonic regime. This is contrary to the conventional understanding of grain growth, in which the variance of f (r, t) is unaffected assuming a homogeneous gas density distribution ρ(x, t). In such a case, the radius a of each grain grows by the same amount ∆a in a given time ∆t, which means the absolute variance of the GSD is conserved. It is noteworthy that σ a evolves with time and increases with M rms and σ ρ as a consequence of the accretion process only, in the presence of turbulence.
Previous work on GSD evolution (e.g., Hirashita 2012;Kuo & Hirashita 2012;Asano et al. 2013aAsano et al. ,b, 2014 have, for the most part, not considered the compressibility and density variations of the ISM, which is particularly strong in MCs. However, more recent work has been directed towards grain processing in simulated gas flows (see, e.g. McKinnon et al. 2018;Kirchschlager et al. 2019;Sumpter & Van Loo 2020), but we still have a lot to learn. The compressibility of interstellar gas accelerates the mean growth rate and the width of the GSD increases. Therefore, it seems recommendable that gas-density variance due to turbulence should be taken into account in future studies of dust evolution in galaxies, at least in the form of a parametric prescription.
As discussed in Mattsson (2020), the mean growth rate by accretion can be accelerated by an order of magnitude or more if M rms ∼ 10, which is quantitatively confirmed by numerical simulations in the present study. Such an enhancement due to the compressibility of interstellar gas makes it possible to deplete essentially all growth-species molecules in an MC onto dust grains within the expected lifetimes of MCs. This indicates that a turbulence-facilitated accretion process can produce a sufficient amount of dust mass to match the observed dust fractions at high redshifts (z = 7−8), without SNe as extreme dust producers and assuming low dust destruction rates as in, e.g., (Dwek et al. 2007). In addition, M rms ∼ 10 also leads to roughly a factor of 10 broadening of the GSD. This can potentially accelerate other dust growth processes such as coagulation (Li & Mattsson 2020), leading to the formation of very large grains (a > 10 µm).
Depletion is excluded in this study, but we must emphasize that the depletion time scale is of the same order (or less) as the expected lifetime of a MC. In the framework of the M20 theory, assuming dust sublimation is negligible, the rate of mass increase of an element i locked up in dust can be expressed where Υ i is the mass fraction of i in the generic dust species. The absolute depletion δρ i = ρ i − ρ d, i of an element i takes place on timescale which can be estimated as τ dep = δρ i (0) |dδρ i /dt| −1 . Assuming a constant ρ i , we also have that δρ i (0) = ρ i . Thus, τ dep ∝ (1 + σ 2 s ) −1 , which means that the depletion time decreases quickly as the density variance increases. In case of a homogeneous gas medium, the depletion time is to first order just inversely proportional to the gas density, corresponding to 10 6 -10 7 yr for a typical MC (see, e.g., Boland & De Jong 1982). This is roughly of the same order as the expected lifetime of an MC (see M20). From this we may conclude that many elements can be fully depleted well within the lifetime of a MC and that it may be important to implement depletion in future work.
The GSD is a crucial parameter, as it determines the slope and shape of the ISM extinction curve and also affects the flux-to-mass ratio and temperature distribution of a dust population (see Mattsson et al. 2015, and references therein). We find that the effect of gas-density variance produce a lognormal GSD regardless of the shape of the initial distribution if the grain-size evolution is dominated by growth from molecule accretion. This is due to the lognormal gasdensity distribution generated by rotationally forced isothermal turbulence (Mattsson 2020, submitted). To our knowledge, there is no directly observational evidence of a lognormal GSD. However, Weingartner & Draine (2001) found that the observed galactic emmision is best reproduced if the size distribution of small dust grains is lognormal.

SUMMARY AND CONCLUSIONS
Numerical simulations covering a wide range of M rms (1.38-11.0) have been performed to study the growth of dust grains by accretion of molecules in ISM turbulence. The size evolution of the dust grains was tracked in a Lagrangian manner. Consistent with the theory developed in M20, we have confirmed that the mean growth rate of dust grains is proportional to the second moment of the gas density PDF, i.e., d a /dt ∝ ρ 2 . In the M rms range explored here, ρ 2 / ρ 2 reaches a more or less steady state, which results in a linear mean growth rate of the radius a .
In case of high Mach number turbulence, the simulation results deviate slightly from the theory of Mattsson (2020) due to the decoupling of dust grains from the gas, which, on the other hand, is consistent with the theoretical predictions by Baines & Williams (1965). Furthermore, we have found that strong variations of gas density due to turbulence in the ISM leads to significant broadening of the GSD. This result is indeed noteworthy and, in fact, contrary to the conventional understanding that the absolute width of the GSD is essentially unaffected by growth due to molecule accretion, assuming a homogeneous gas density field. Larger M rms obviously results in stronger density variations and therefore a wider GSD.
It is indeed remarkable that the simple M20 theory is able to capture the effects of gas-density variance on the accretion process so accurately. This provides a possibility to prescribe turbulence accelerated dust growth in models of galactic dust evolution without much loss of generality. Moreover, we found also that when gas density variations were taken into account, the growth process leads to a lognormal GSD, regardless of the shape of the initial distribution. This is a novel result and likely a reflection of the gas-density PDF, as shown by Mattsson (2020, submitted) in a recent attempt to construct a mathematical theory of the phenomenon. As mentioned above, a follow-up study based on this result in the present study is made by Mattsson (2020, submitted), which shows that the lognormal GSD can be derived by modelling turbulent gas-density variation as a Markov process.
The thermal velocityū t of molecular clouds is assumed to be constant in the present study. However, Eq. (15) shows thatū t ∝ m −1/2 . Considering a more realistic chemical composition of MCs could reduce the growth rate. But the overall enhancement of the growth rate due to turbulence should not be affected as it is independent of the chemical composition. Depletion of growth-species molecules is not included in this study (to save computing time), but according to Mattsson (2020) the predicted net effect is merely to slow down the growth rate. However, including depletion and other types of grain processing is a currently ongoing project, which we will return to as soon as possible.