Rapid Elimination of Small Dust Grains in Molecular Clouds

We argue that impact velocities between dust grains with sizes less than $\sim 0.1$ $\mu m$ in molecular cloud cores are dominated by drift arising from ambipolar diffusion. This effect is due to the size dependence of the dust coupling to the magnetic field and the neutral gas. Assuming perfect sticking in collisions up to $\approx 50$ m/s, we show that this effect causes rapid depletion of small grains - consistent with starlight extinction and IR/microwave emission measurements, both in the core center ($n \sim 10^{6}$ cm$^{-3}$) and envelope ($n \sim 10^{4}$ cm$^{-3}$). The upper end of the size distribution does not change significantly if only velocities arising from this effect are considered. We consider the impact of an evolved dust size distribution on the gas temperature, and argue that if the depletion of small dust grains occurs as would be expected from our model, then the cosmic ray ionization rate must be well below $10^{-16}$ s$^{-1}$ at a number density of $10^{5}$ cm$^{-3}$.


Introduction
Knowledge of the dust size distribution is crucial to interpretation of observations in molecular clouds. There is evidence that the size distribution in dense regions contains fewer small grains (Weingartner & Draine 2001) than the dust in the diffuse ISM and some grains of at least micron size (Pagani et al. 2010). There has been theoretical work modeling the coagulation process. Ossenkopf (1993) made estimates for the contribution of several processes, including turbulence, Brownian motion, forces arising from grain asymmetries and gravitational settling, to the relative velocity between dust grains, and determined turbulent motions to be dominant at densities ≤ 10 8 cm −3 . He ran coagulation simulations for regions with densities between 10 5 and 10 9 cm −3 . It was found that the optical properties of the dust distribution change significantly due to the removal of small particles, but the upper edge of the size distribution does not change appreciably within 10 5 years for a density of 10 6 cm −3 . Ormel et al. (2009) used a detailed model for the outcome of grain-grain collisions, taking into account the effect of collisions on both grain size and structure. They assumed collision velocities to be due to gas turbulence, as this was argued to be dominant in Ossenkopf (1993), and they ignored the effect of grain charge on the collisional rates. A great deal of uncertainty exists regarding the lifetimes of molecular cloud cores (Lee & Myers 1999;Ciolek & Basu 2000), although a few million years seems to be a good estimate based on statistical studies by e.g. Könyves et al. (2015). The lifetime of the central regions of a dense core (with densities above 10 5 cm −3 , as in prestellar cores; e.g. Crapsi et al. (2005)), similar to that modeled in this paper, is likely somewhat shorter. Ormel et al. (2009) ran their standard simulation for 5 × 10 7 years with a gas density of 10 5 cm −3 . In this longer time they found that large (tens of µm) aggregates form, and fragmentation becomes important.
Recent non-ideal MHD simulations have pointed out the importance of the removal of very small grains for the formation of protoplanetary disks (Zhao et al. 2016(Zhao et al. , 2018b. It is therefore important to look into this problem quantitatively and for physical conditions relevant for stellar system formation. In particular, we need to determine the time scale for very small grain removal in prestellar cores, which represent the initial conditions in the process of star and planet formation. If this depletion time scale is significantly shorter than the dynamical time scale, then the diffusion of magnetic fields (by ambipolar diffusion and the Hall effect) becomes efficient during the protostellar collapse, which greatly promotes the formation of protoplanetary disks. The depletion of very small grains also affects the chemical evolution of prestellar and protostellar cores. Very small grains both provide a large area for surface chemical processes to proceed and also can become the main carriers of electrons, thus affecting ion-molecule chemistry in the gas phase.
In this work, in addition to turbulence, we study the effect of a different source of relative velocity between grains. Ciolek & Mouschovias (1996) considered the role of ambipolar diffusion in changing the spatial distribution of grains. In a collapsing cloud, the large grains are pulled to the center of the cloud by the infalling neutral material, whereas the small (mostly negatively charged) grains stay coupled to the magnetic field. In this paper, we note that this differential motion leads to a substantial additional source of collision velocity between small and large grains, which operates even in the absence of turbulence. This is particularly noteworthy given evidence (e.g. Fuller & Myers 1992;Pineda et al. 2010) that turbulent velocities in dense regions are substantially subsonic.

Cloud Model
We take the cloud density profile which Tafalla et al. (2002) used to model the prestellar core L1544. This is a spherically symmetric cloud with the H 2 number density given by n(r) = n c 1 + (r/r 0 ) 2.5 . (1) Here n c = 1.4 × 10 6 cm −3 is the central density of the core, and r 0 = 0.014 pc = 2,900 AU is the rough scale of the density peak. In this paper, we focus on two representative environments within the cloud. First we consider the location at r = r 0 , taken to be representative of the environment near the cloud center, but far enough from the center that motions due to ambipolar diffusion and gravitational settling are still present. We denote this our "core center" location. We also consider a location at r = 7.2r 0 where n = 10 4 cm −3 , denoted our "envelope" location.
In what follows, we assume that the primary ion is H + 3 in the core center (Caselli et al. 2003) and HCO +  in the envelope. We assume a composition with a helium to hydrogen atom number ratio of 10% and ignore all other species. From the temperature profile shown in the top panel of Figure  3 in Keto & Caselli (2010) we take the temperature to be 6.5K in the core center and 10K in the envelope, consistent with gas measurements done by Crapsi et al. (2007). We take the ion density from Caselli et al. (2002), (their model 3), but scale it with the square root of the cosmic ray (CR) rate ζ: (2) ζ is rather uncertain, but thought to be a decreasing function of the column density N. We take the model for ζ given in Padovani et al. (2018) (assuming their high proton spectrum), which yields ζ = 10 −16 s −1 in the core center, and ζ = 3 × 10 −16 s −1 in the envelope. This is somewhat higher than previous estimates in the literature. Their model however gives the correct ionization rate at lower densities where ζ can be determined from H + 3 observations (Neufeld & Wolfire 2017). A summary of the important parameters in each environment is provided in Table 1, with references provided at the bottom.

Collision Rates and Velocities
We consider both random and systematic grain motions. Random motions are due to Brownian motion (important for the smallest grains) and motions due to the coupling with gas turbulence, important for larger grains. Systematic grain motions arise due to the fact that grains of different sizes and charge states couple differently to the neutral gas, and to the magnetic and gravitational fields.

Grain Charging
As we will show below, both the collision cross sections, and the grain dynamics are affected by the grain charge. We used the grain charging model from Ivlev et al. (2015). This includes a combination of plasma charging (including the effect of the charge-induced dipole interaction on the plasma charging rates), photoelectric charging from the UV photons produced by the interaction of CRs and H 2 molecules (Prasad & Tarafdar 1983), and additional negative charging from collisions of dust with CR electrons in the low-energy tail of the spectrum. However, instead of calculating the local electron spectra ourselves, we estimated them at the three column densities shown in Figure 2 from Ivlev et al. (2015) and interpolated between them. This will not affect the accuracy of our results, since as demonstrated in Ivlev et al. (2015), direct charging from accumulating low energy CRs has only a very small effect relative to the other charging mechanisms. Table 2 shows the charge states of the grains in the envelope. Each row corresponds to a different grain size. The triples of numbers comprising each entry of the table are the fraction of grains which are negatively charged, neutral, and positively charged. The three columns are for different ionization rates ζ in units of ζ P18 , the value shown in Table 1 for the envelope. The charge distribution is dominated by grains of single charge

Brownian Motion
We assume that each dust grain has a thermal velocity distribution which leads to expected collision velocities where µ = m 1 m 2 /(m 1 + m 2 ) is the reduced mass of the two dust grains.

Motion Due to Coupling with Gas Turbulence
If transonic turbulence is present, motions of dust grains arising in response to this turbulence likely dominate the collision velocities between the larger grains (Ossenkopf 1993). To calculate these motions, we use the results from Ormel & Cuzzi (2007). We assume a Kolmogorov turbulent energy spectrum with no smallest scale, i.e. E(k) ∝ k −5/3 for all k greater than k L . Here k is the spatial frequency, and k L is the k corresponding to the turbulent injection scale. We take the rms collision velocity from Ormel & Cuzzi (2007). We use their equations (5)-(10), in the limit of infinite Reynolds number. This gives the collision velocity as a function of the turbulent velocity v g , and the Stokes numbers of the two grains, defined as S t = τ s /τ L , where τ L is the turbulent crossing time at the injection scale, and τ s is the stopping time of the grain, given by where v * th = 0.92 4k B T/(πm p ) is the thermal velocity scale of the gas (with a factor to account for the helium fraction), ρ g = 2.8m p n is the gas density, and σ g = πa 2 (where a is the grain radius) is the collision cross-section between grains and gas particles.
This treatment likely overestimates the collision velocities between smaller dust grains for two reasons. First there must be a cut-off at some k due to either ion-neutral friction, or neutral viscosity (Xu et al. 2016). Also, even if there isn't a sharp cutoff, there is evidence from simulation (Downes 2012) that the spectrum is steeper than Kolmogorov in realistic ISM conditions. For this reason, we consider models both with and without turbulence, and focus much of our paper on the previously ignored, but actually better constrained, source of relative velocity discussed in the following section.

Systematic Drifts
In this section we consider systematic drifts due to the sizedependence of the grain coupling to the neutral gas, and gravitational and magnetic fields. Ossenkopf (1993) considered the differential drift of grains in a gravitational field and found this effect to be subdominant to the turbulent motions. He was, however, considering turbulence which was sonic at the Jeans length, and there is evidence (e.g. Fuller & Myers 1992;Lada et al. 2008;Pineda et al. 2010), that turbulent velocities in dense cores are much lower. In this regime, gravitational settling can be the dominant source of collision velocity for grains larger than a few tenths of a micron. However, as we show later, this is unlikely to have a significant impact on the size distribution. Ciolek & Mouschovias (1996) considered the drift of small dust grains relative to large ones as a result of their coupling to the magnetic field. To our knowledge this motion has not been considered in studies of coagulation. In this section we show this effect can be the dominant source of relative motion for grains up to ∼ 100 nm, even when significant turbulence is considered.
To estimate the drift rate of the dust grains as a result of coupling to the magnetic and gravitational fields as well as gas drag, we assume that the magnetic field is able to keep the ions stationary in the frame of the cloud center. We then assume that the motion of the neutrals can be determined by balancing the drag force from collisions with ions against the gravitational force. Considering both the molecular hydrogen and the helium, the drift velocity v in between the ions and neutrals is determined by c d gρ g = n i v in µ H 2 n H 2 σv H 2 + µ He n He σv He , where g is the local gravitational field, and µ H 2 and µ He are the reduced masses of H 2 and helium respectively in collisions with the ions. σv is the ion-neutral momentum transfer cross section averaged over a Maxwellian velocity distribution. Following Raizer et al. (2011), σv is given by 4πa 2 0 (α/a 3 0 )I H /µ, where a 0 is the Bohr radius, I H the Rydberg energy, and α is the polarizability of the neutral species. α/a 3 0 = 5.52 for H 2 and 1.39 for He. We have added a fudge factor c d ∈ [0, 1] related to the fact that the field is not purely azimuthal, the ions may be to some extent dragged in along with the neutrals, or the cloud may be supported by gas pressure in addition to the magnetic field.
Solving Equation (5) (with c d = 1) for v in yields v in = 5.2 × 10 3 cm/s in the core center and 4.2 × 10 3 cm/s in the envelope.
In order to determine the grain dynamics, we must understand how the charging timescales compare with the other two relevant dynamical timescales in the problem -the stopping time τ s and the inverse of the grain gyrofrequency Ω g , which is given by Here,q is the time-averaged grain charge, B the magnetic field strength, and c the speed of light. For the sizes where these  Fig. 1. Timescales for processes relevant to grain dynamics. τ s is the stopping time of the grain, given by Equation (4). Ω −1 g is the inverse grain gyrofrequency given by Equation (6). τ ch (0) is the time for a neutral grain to acquire a negative charge, and τ ch(−1) is the time for a grain of charge -1 to become neutral. The top panel is for the core center, and the bottom is for the envelope.
timescales are comparable, nearly all grains are either singly negatively charged, or neutral (see Table 2, and note that there are essentially no multiply negatively charged grains). For this reason, we consider two charging timescales: τ ch(−1) is the time for a grain of charge -1 to become neutral, and τ ch(0) is the time for a neutral grain to acquire a negative charge.
The charging timescales, as well as τ s and Ω −1 g are plotted in Figure 1 as a function of grain size. We calculated τ ch(0) and τ ch(−1) as a function of grain size using Equations (3.1) -(3.5) from Draine & Sutin (1987). These equations ignore the contribution of UV photons from H 2 fluorescence due to CR impacts, but this is not important to the charging of the grain sizes where the charging timescale is comparable to the other dynamical timescales. The top panel is made for the core center, and the bottom is for the envelope. We see that for sizes greater than about 10 nanometers, in both environments, the two charging timescales are shorter than both τ s and Ω −1 g . For that reason, for the calculations in this paper, we use the mean grain charge when calculating the grain dynamics. This approximation should somewhat overestimate the depletion timescale for grains less than 10 nm, as in reality there may be significant collision velocities between grains of the same size, but different charge states. We also performed a simulation in which we considered the opposite limit where we calculate the grain dynamics assuming that an individual grain's charge remains fixed, and found indeed that the smallest grains were depleted even faster, and there was little difference for the larger grains.
To estimate the grain velocities, we assume the following model. Consider a coordinate system in whichr points radially outward, and the magnetic field has no radial component. The ions are assumed to be at rest with respect to the cloud center. Then, the equations of motion for a dust grain with mean charge q are given by where v ⊥ is the velocity in the direction parallel toB ×r. In equilibrium, these equations can be solved to yield where g eff = g + v in /τ s is the effective gravity. We note that for the grain sizes considered in this work, v in /τ s g for c d ≤ 1. The collision velocity due to these drift effects is Small dust grains (those for which Ω −1 g τ s ), are relatively well-coupled to the magnetic field so they move with respect to the gas at velocity v in in ther direction. Large grains (Ω −1 g τ s ) do not feel the magnetic field, and drift in the −r direction at speed gτ s . The transition between these two regimes occurs near size a c , defined as the size where Ω g τ s = 1. Using the expressions for Ω g and τ s , we find We observe that the RHS of Equation (10) depends also on a c throughq. The mass of the grain does not enter Equation (10), so the expression for a c is valid for very porous grains as well.
To estimate B(r), we assume that the core is in hydrostatic equilibrium with the pressure support provided by the magnetic field. Thus B(r) solves the equation where we assume that B outside the cloud is negligible, and g(r) is calculated from the density distribution ρ g (r) determined by Equation (1). This gives B of 395 µG in the core center and 45 µG in the envelope, corresponding to a c of 34 and 88 nm in the two environments respectively. Nakamura et al. (2019) measured a line of sight magnetic field of 117 ±21 µG at an estimated H 2 density of 3 × 10 4 cm −3 in the cloud TMC-1. Multiplying by √ 3 to estimate the total magnetic field gives 203 µG. Our model gives 81 µG at this density.
There are two issues with this expression. First, it does not necessarily accurately represent the mean collision velocities between bodies in which more than one source of relative motion is important. Second, it does not account for the fact that if turbulent motions are important, there is in fact a distribution of collision velocities. Given the large uncertainties on the individual sources of relative motion, both these issues seem rather minor.  (5)], are varied between the panels, as labelled. c s is the isothermal sound speed, given by k B T/(2.33m p ). We have included labelled contours at 1, 10, and 50 meters per second. Figure 2 shows the collision velocities as a function of the sizes of the collision partners. As we are uncertain about the strength of the turbulence and the ambipolar diffusion, the four panels show the four possible cases in which each effect is either present or absent. We see that collision velocities are generally largest for collisions in which one grain is larger than a c and one is smaller, or between two large ( 1µm) grains in the models in which turbulence is present. Figure 3 shows the dominant source of collision velocity for each set of collision partners. Grain collision velocities are dominated by turbulence and ambipolar diffusion. Gravitational settling is dominant for a wide range of particle sizes only if there is no turbulence and we ignore ambipolar diffusion. Brownian motion is typically irrelevant for all except the smallest particles. Turbulence and Brownian motion are able to excite collision velocities between grains of the same size. Systematic drifts are not. For that reason, we see narrow grey and black lines across the panels, corresponding to samesize collision partners where the systematic drifts have no effect.

Collision Rates
In our dust coagulation code, the size distribution is discretized in N bins discrete logarithmically spaced size bins between r min and r max . We can then write the number of collisions per unit time between bodies in mass bin i and j in a given thin spherical 0.01 shell as Here N i and N j are the numbers of particles in mass bin i and j respectively, v i j (q i , q j ) is the collision velocity between a particle in mass bin i with charge q i and a particle in mass bin j with charge q j . Note that for our standard assumption that the charging timescales are rapid compared to the other dynamical timescales, v i j is calculated usingq i andq j , however σ i j is still calculated using the instantaneous charge states q i and q j . f q i is the fraction of particles in mass bin i with charge q i , and f q j is defined analogously for mass bin j. V is the volume of the shell. The collision cross section is given by where a i and a j are the radii of the two dust grains. We find in Section 5.3 that Coulomb attraction/repulsion plays a negligible role in the environments we consider 1 . This is to be expected, given the typical collision velocities. As can be seen from Equation (14), the critical parameter determining the Coulomb effect is 2q i q j / µv 2 i j (a i + a j ) . For 10 nm grains, this is less than unity for v > 2 m s −1 . It scales inversely with the 4 th power of the sizes of the colliding grains, so it will be smaller for larger grains. Figure 2 shows that the typical collision velocities are substantially larger than 2 m s −1 except in the case with no turbulence or ambipolar diffusion.

Characteristic Timescales
In this section we estimate four relevant timescales. We consider τ AD , the time for small grains to collide with large ones due to ambipolar diffusion, τ grav , the growth timescale of the largest grains due to gravitational settling, τ turb , the growth timescale of the largest grains due to turbulent motions, and finally τ cl , the lifetime of the cloud against ambipolar diffusion, which we take as the runtime of our simulations.

Cloud Lifetime
We assume that the lifetime of the cloud at scale r is equal roughly to the ambipolar diffusion time at that scale in the cloud (Tassis & Mouschovias 2004), i.e., the cloud lifetime τ cl is given by This works out to be 2.6 × 10 5 years in the core center, and 2.3 × 10 6 years in the envelope, assuming the standard value for the ionization fraction given in Equation (2).

Elimination of Small Grains Due to Ambipolar Diffusion Drift
We can estimate the timescale on which grains coupled to the magnetic field collide with grains coupled to the neutral gas as follows. First, by setting the mass in dust equal to a fraction f d of the gas mass, and assuming an MRN (Mathis et al. 1977) distribution between sizes a min and a max (dn/da ∝ a −3.5 ), we can write the differential density of dust grains with radius a as f d = 0.01 is the gas to dust mass ratio. Then, let us consider a small grain of size a s < a c . We will approximate its collision velocity with another larger grain of size a l as 0 for a l < a c and v in for a l > a c . We will approximate the collision cross-section as πa 2 l . Then the expected time for the smaller grain to experience a collision with a larger one is given by This gives For the parameters used here, we find τ AD = 1.1 × 10 4 years in the core center, and 2.4 × 10 6 years in the envelope.

Coagulation of Large Grains Due to Gravitational Settling
We can also estimate the time required for the upper end of the mass distribution to grow due to collision velocities arising from gravitational settling. We assume in this case that all grains are large enough that we may ignore the contribution of the magnetic field. In this case, the appropriate limit of Equation (8) (Ω g → 0) yields a collision speed between two particles of sizes a s and a l of Then, assuming the collision cross-section to be πa 2 l , dropping the term a s in Equation (19), and taking for simplicity the lower end of the mass spectrum to be 0, we can calculate a growth timescale for bodies initially at the top end of the mass distribution: This expression evaluates to 1.9 × 10 6 years in the core center and 8.4 × 10 6 years in the envelope.

Coagulation of Large Grains Due to Turbulence
In the limit of small Stokes numbers, one can show using Equations 6, 7, & 10 from Ormel & Cuzzi (2007) that the relative velocity between two grains of the same size is where S t l is the Stokes parameter of the larger grain, is the ratio of the size of the smaller grain to that of the larger, and ψ( ) varies between 1.95 for = 1 and 2.95 for = 0. To calculate the growth timescale, we assume that growth is dominated by collisions of similar-sized grains. We therefore assume a cross section of 4πa 2 l , and calculate the velocities using Equation (21) with ψ = 1.95. These approximations yield a timescale for growth of the largest grains from collisions arising from turbulence: r enters here through the assumption that it is equal to the turbulence injection scale. Assuming v g = c s , where c s = k B T/(2.33m p ), this expression evaluates to 9.4 × 10 4 years in the core center, and 2.0 × 10 6 years in the envelope. Given our belief that turbulent velocities near the core center are subsonic, the core center timescale represents a lower limit to the probable actual coagulation timescale.
Using Equation (22) in conjunction with the fact that a/ȧ = 3m/ṁ, we can solve for a l (t), assuming that a l (t 0 ) = a max : We show in Section 5.3 that this is a reasonable approximation, but yields a slight overestimate.

Core Center
We consider here the coagulation of dust at radius r 0 from the center. We consider four models with different values of v g and c d , the parameters determining the strength of the turbulence, and the ion-neutral drift velocity. The results are shown in Figure 4. The values of v g and c d in each panel are equal to those in the corresponding panel of Figure 2. Different color curves correspond to the size distribution at the time shown in the legend. The dashed black curve is the MRN distribution, used as the initial condition.
The simulation was run for time τ cl (see Section 4.1), equal in this case to 2.6 × 10 5 years. In all cases, for times less than τ cl /3 years, there is very little evolution of the upper end of the size distribution. The evolution of the high end of the size distribution depends most on the value of v g , which makes sense, as we see in Figure 3 that collisions between bodies with sizes around 1 µm depend primarily on the strength of the turbulence. In Section 4, we derived growth timescales of 9.4 × 10 4 years and 1.9 × 10 6 years (40% and 730% of τ cl ) for growth of the upper end of the size distribution due to turbulence and gravitational settling respectively, assuming sonic turbulence. This is consistent with the growth of the peak of the distribution to about a micron for the simulations with sonic turbulence, and only minor changes to the upper scale of the distribution for all simulations with subsonic turbulence. Note from Equation (22) that τ turb ∝ v −3/2 g . Examination of the bottom panels of Figure 4 shows that in the absence of turbulence, ambipolar diffusion can lead to a minor change in the upper edge of the size distribution, but cannot produce micron sized grains.
We find when c d = 1 that there is a rapid depletion on ∼ 10 4 year timescales of grains less than a c (defined by Equation 10), consistent with the estimate in Equation (18). By looking at the bottom right panel of Figure 4, we see that Brownian motion has only a minor effect, and only for the very smallest grains.

Envelope
We also consider coagulation in a lower density environment (n H 2 = 10 4 cm −3 ). This corresponds to a location at r = 20,800 AU in the model of Tafalla et al. (2002). Figure 5 is analogous to Figure 4, but made for the envelope. As discussed in Section 4, the cloud lifetime τ cl is higher by nearly a factor of 10 in this environment. This does not, however, lead to more growth, because τ turb is higher by a factor of ∼ 20, and τ AD higher by a factor of 200. τ grav is higher by only a factor of 4, but is still longer than τ cl . For this reason, there is somewhat less growth of the top of the dust distribution, and significantly less complete removal of the small grains.
The solid black curve is the curve derived by Weingartner & Draine (2001) for dense regions (their "case A" model with R v = 5.5 and 10 5 b c = 3.0). It is not clear precisely which region of the cloud these observations would be expected to best correspond to. We see none of the models are able to reproduce the absence of the smallest grains from Weingartner & Draine (2001). The ones which do the best are those with high values of c d . There is also only minimal growth of the large end of the spectrum, consistent with the model from Weingartner & Draine (2001), but suggesting that the coreshine observations from Pagani et al. (2010) must have come from dust in regions denser than 10 4 cm −3 or that our estimate of the cloud lifetime is too short (see next subsection).

Effect of CR Ionization Rate
CRs can affect the coagulation process in two ways. First, they are the primary cause of gas ionization, and therefore determine the ambipolar diffusion velocity v in defined in Equation (5). Since the ionization fraction is assumed to be proportional to √ ζ, this means that v in ∝ 1/ √ ζ. As the cloud lifetime, τ cl is  Here τ cl = 2.6 × 10 5 years. In each case, the initial distribution is an MRN distribution from 5 to 250 nm, which is shown in the black dashed line.
proportional to 1/v in (see Section 4.1), there is more time for coagulation to occur at higher ζ. We note, however, that τ AD also scales as 1/v in (see Section 4.2). Hence, to zeroth order, the degree to which small grains are eliminated in the time τ cl due to ambipolar diffusion drift is independent of ζ. Second, as discussed in Section 3.1, CRs also excite electronic states of H 2 which decay, producing UV photons, which lead to more positively charged dust grains (Prasad & Tarafdar 1983;Gredel et al. 1989;Ivlev et al. 2015).
To examine the effect of changing ζ on the evolution of the size distribution, we ran the coagulation simulation for different CR ionization rates. These runs were done for the envelope, where the CRs have more effect on the dust charge distribution. The results are shown in the top two panels of Figure 6. Different curves correspond to different assumed ionization rates, and simulation run times. The simulation run time is expressed in units of the cloud lifetime τ cl , given in Equation (15). Because v in ∼ ζ −1/2 , then the cloud lifetime is proportional to ζ 1 2 . The top panel assumes sonic turbulence and c d = 0. In this case, we see that the coagulation process is essentially unaffected by ζ, except that a higher ζ means a high cloud lifetime in our model. This shows that the Coulomb focussing/repulsion issue is completely unimportant, making only very minor changes to the size distributions at the smallest sizes. The middle panel assumes no turbulence, and c d = 1. In this case, naively one would expect little difference in the degree of coagulation within the cloud lifetime, since the collision velocities are dominated by ambipolar diffusion, so the cloud lifetime is inversely proportional to the scale for the collision velocities. In fact, we see that over the lifetime of the cloud, the high ζ model still gets rid of more of the small grains. This is because a higher CR rate decreases the mean grain charge for grains with sizes between 10 and 100 nm (see Table 2), and hence a c decreases too (see Equation (10)). Equations (15) and (18) demonstrate that τ AD /τ cl ∝ √ a c , and therefore also decreases with increasing ζ.
The ability to grow micron-sized grains depends sensitively on the cloud lifetime τ cl . Equation (23) shows that the peak size roughly scales as the square of the elapsed time for large times if growth is dominated by turbulence. The peak of the grain mass distribution (peak of m 2 dN/dm) is shown by the blue curve in the bottom panel of Figure 6 and the approximation given by Equation (23) is shown in orange. The orange curve gives a slight overestimate for two reasons. First, there is initially no move- Observation = 10 P18 , t = cl (0.1 P18 ) = 10 P18 , t = cl ( P18 ) = 10 P18 , t = cl (10 P18 ) Fig. 6. The top panel of this plot shows the size distribution for the case with sonic turbulence, and c d = 0, in the envelope. As shown in the legend, the different curves correspond to different amounts of time and different CR rates. Time in this case is expressed in terms of τ cl for each different CR ionization rate. τ cl (0.1ζ P18 ), τ cl (ζ P18 ), and τ cl (10ζ P18 ) are the cloud lifetimes, evaluated assuming CR ionization rates of 0.1ζ P18 , ζ P18 , and 10ζ P18 respectively, where ζ P18 is the ionization rate calculated from (Padovani et al. 2018). In the case of c d = 0 (top panel) we see that the curves at different times are nearly independent of the ionization rate. For this reason, the curves overlap nearly perfectly, and only one curve per time is visible. The middle panel assumes no turbulence, but assumes c d = 1. The blue curve in the bottom panel shows the peak of the distribution as a function of time (assuming ζ = ζ P18 , c d = 0, and v g = c s ). The orange curve in the bottom panel shows the approximation given by Equation (23). ment of the peak due to the sharp cutoff of the original size distribution. Additionally, Equation (23) assumes that all collisions are between objects at the peak size, so the cross-section is slightly overestimated. Nonetheless, it is clear that were the cloud to survive a few times longer than in our model, the peak size would be substantially larger.

Evolution in Both the Core Center and Envelope
In reality, material goes through a low-density phase better described by our envelope model before becoming part of the core center. To model this, we took the output of our envelope simula-tions as the initial input to the core center simulations. We found however that the pre-processing in the envelope made very little difference to the resulting distribution after coagulation within the core center. This makes sense, as the evolution is qualitatively similar, but lifetime of the core center is longer in units of the local collision timescale than the lifetime in the envelope in units of the envelope collision timescale.

Contribution of Small Grains to Ambipolar Diffusion Velocity
In Equation (5) we calculate the ambipolar drift speed, assuming that neutrals are slowed down only through collisions with ions.
In fact, neutrals must to some extent be slowed by collisions with small magnetically coupled dust grains as well. To calculate this, one can add a term a max a min to the RHS of Equation (5), where v dn (a) is the speed of a dust grain with radius a with respect to the neutral gas. This is not a significant effect in the envelope. Even if all the dust grains were completely coupled to the magnetic field, their addition to the RHS of equation (5) would be smaller than that of the ions by more than a factor of 5, assuming an MRN distribution.
On the other hand, this additional coupling may be significant in the core center where the ion density is significantly lower. Let us denote the ratio of the expression in Equation (24) to the RHS of Equation (5) as χ. An upper bound on χ can be found by assuming an MRN distribution, and that the grains are fully coupled, so v dn = v in . In this case χ = 5. In reality for an MRN distribution, χ = 3.5 because the larger grains are uncoupled. Nevertheless, this is still a significant correction. However, the problem resolves itself rapidly. After a time τ cl /10, the simulation with no turbulence, but c d = 1 has removed a sufficient number of the small grains so that χ = 0.4. Of course the removal is somewhat slower (by a factor of 1 + χ) in reality when the small magnetically coupled grains are still present, but it seems reasonable that for most of the evolution time of the cloud the magnetically coupled grain population would be small enough that we are justified in ignoring this correction (or lumping it together with the fudge factor c d ). Additionally, Figure 5 shows that there is a substantial reduction in the number of small grains even in the envelope stage, so by the time grains reach the core center, the density of magnetically coupled grains should be reduced by a factor of several.

Icy Mantles and Non-Compact Grains
In the interests of simplicity, the preceding analysis has left out two important effects which may enhance dust coagulation. First is the fact that at densities higher than ∼ 10 5 cm 3 , nearly all of the carbon and oxygen in the gas phase freezes out onto the dust grains (Caselli et al. 1999). This results in a substantial increase of the grain volume. To estimate this, we assume that all the carbon and oxygen in the ISM condenses onto the surface of the grains. We assume that the ratio of carbon and oxygen to hydrogen is the same as in the Sun. Using the abundances given in (Asplund et al. 2009), and assuming a density of 1 gram per cm 3 for this mantle material, we find that the total volume of the grains is increased by a factor of V rel = 2.8. This factor is substantially uncertain; for reference, Ossenkopf (1993) considers two cases in which the mantle component has 0.5 times and 4.5 times the volume of the refractory cores.  Figure 4, but assuming grains to be covered with icy mantles, and to be non-spherical/porous (which increases their projected surface area by a factor of 2). As in Figure 4, the different panels are made for different degrees of turbulence and ambipolar diffusion, parameterized by v g and c d , as labelled on the panels. The sharp peak around 20 nm is the result of our assumption that mantles of equal thickness were applied to all the grains in the MRN distribution. Different ways of assigning the mantle volume do not appreciably alter the results at later times. The wiggles seen in some of the curves near 20 nm are a simulation effect to do with the sharp cutoff of the initial mass distribution and the finite bin-width.
Additionally, realistic grains are not perfect compact spheres, and therefore have a higher projected area to mass ratio. Shen et al. (2008) consider model grains created through ballistic aggregation followed by compaction. They find, depending on the model, a ratio of angle-averaged projected area to the area of a solid sphere of the same material between 1.4 and 4. In this section we take this to be a factor of 2. To account for this, we increase our collision cross sections, and decrease our stopping time by a factor of 2. Figure 7 is analogous to Figure 4, but as discussed above, we have included the icy mantles and considered the effect of nonspherical/porous grains on the dynamics and collision rates. We use a couple of approximations in our treatment of the mantles. At t = 0, we assume that the refractory cores have an MRN size distribution. In addition, each grain gets a mantle of a constant thickness a mantle = 13nm, which is chosen so that the ratio of the mantle volume of the whole grain population to the core volume of the whole grain population is V rel = 2.8. This leads to the majority of the volume being in grains of size between 1 and 2 times the mantle thickness, as can be seen in Figure 7. Constant mantle thickness is the result one would get if there were no grain size dependence to the freeze-out process. In reality, the smaller grains may have their mantles evaporate due to transient CR heating (Leger et al. 1985;Zhao et al. 2018a). We find that the initial way the mantle material is distributed does not matter. At time τ cl /10, all memory of the spike in the size distribution at the smallest sizes has been removed for the cases with c d = 1, and even with just turbulence, the small sizes no longer dominate the grain volume. We verified that the way the mantle volume is initially distributed has a negligible effect on the size distribution at later times.
In determining the grain dynamics, we assume that for grains of size between a min and a max , For sizes a > a max the density smoothly approaches the mean density given bȳ Comparing Figure 7 with Figure 4, it is apparent that growth to above a micron within the estimated core lifetime is only possible if we consider the icy mantles and non-spherical/porous grain shapes. With these factors taken into consideration, we find that growth to about 10 microns occurs in our models with sonic turbulence. With just ambipolar diffusion, the differences between Figure 7 and Figure 4 are less extreme. The evolution proceeds faster with the icy mantles, but the end result in both cases is a narrow distribution centered near the top of the initial distribution.

Grain Fragmentation
A significant uncertainty which we have not considered is the possibility of collisional outcomes other than perfect sticking. There have been significant efforts both experimentally and theoretically to characterize the outcome of grain-grain collisions, but it remains extremely uncertain -see Blum & Wurm (2008) for a good review. It is likely that silicate grains would fragment given our collision speeds of up to 50 m s −1 . However, the situation may be different if the grains are coated with icy mantles. Wada et al. (2013) performed simulations of particle aggregates with various compositions, and found that icy aggregates grew in collisions with velocities below 80 (a/0.1µm) −5/6 m s −1 (assuming that the aggregates were composed of spherical monomers with radius 0.1 µm). For bare silicate grains, the critical velocity was an order of magnitude lower. Kimura et al. (2015) found a critical velocity of 4 m/s for 0.1 µm bare silicate spheres, but argued that coating in water ice does not help very much. On the other hand, they also studied a model in which the grains are coated in complex organic mantles, and found a critical sticking velocity of 66 m/s. Ormel et al. (2009), using a model of turbulent collision velocities similar to ours, concluded that ice-coated grains reach ∼ 100 µm before fragmentation prevents further growth, whereas for pure silicate grains this occurs at around a micron. On the experimental side, Gundlach & Blum (2015) find a critical velocity for 1.5 µm H 2 O spheres of 9.6 m s −1 . Smaller particles are thought to be more likely to stick. Extrapolating the experimental result from Gundlach & Blum (2015) down to 0.1 µm spheres (using the ∝ a −5/6 scaling) would suggest a critical sticking velocity of 90 m s −1 . Not all refractory materials, however, have such a high sticking threshold. Musiolik et al. (2016) found a critical velocity of 0.04 m s −1 for 100 µm CO 2 grains, and noted that it will be nearly an order of magnitude lower if these grains are scaled to 1.5 µm.
Given such a wide range of available results in the literature, and noting that we are primarily concerned with grains ≤ 0.1 µm, we neglect the fragmentation effect in the present paper and assume a perfect sticking of colliding grains.

Gas Temperature
It has been shown (Ivlev et al. 2019) that in the core center, the gas temperature is sensitive to the dust size distribution. At these densities, the gas heating is dominated by CRs, which heat the gas at a rate proportional to the ionization rate, and only weakly dependent on other parameters (Glassgold et al. 2012). The gas cools primarily due to collisions with dust grains. If the dust has coagulated into larger bodies, then there is less surface area for the gas to collide with, so the cooling is reduced. In addition to grain cooling, we model the radiative cooling by molecules, which enables us to extend our analysis to lower densities. We calculated the gas temperature using the code presented in Sipilä & Caselli (2018); Sipilä et al. (2019), which combines hydrodynamics with chemical and radiative transfer models in a selfconsistent way. For the present work, we assumed that the core has a static physical structure, so that the gas temperature is determined by balancing CR heating with cooling by the gas-dust collisional coupling and molecular line radiation. Figure 8 shows the gas temperature as a function of local density in our model cloud. This is made for three different dust coagulation models, as labelled on the panels. The top panel shows results for an MRN size distribution with no mantles. In the bottom two panels, the size distribution is determined by a coagulation simulation in which we assume compact spherical grains, no turbulence, and c d = 1. The middle panel shows a model in which grains have mantles as discussed in Section 6.3. The bottom panel is the same as the middle, but we have ignored the presence of icy mantles. In each panel, the black line shows the dust temperature determined by the radiative transfer code CRT (Juvela 2005). The curves labelled "High CR" and "Low CR" correspond to the "high" and "low" models in Padovani et al. (2018), which have ionization rate around 10 −16 s −1 and 2 × 10 −17 s −1 respectively at these densities. We did not run different dust coagulation simulations for the high and low ζ cases, but simply used different models for ζ to calculate the gas temperature. The dotted lines correspond to temperatures from Equation 18 from Ivlev et al. (2019), and the solid lines are the results when line cooling is included as well. Observations of gas temperatures in dense cores suggest that at densities of 10 5 cm −3 , the temperature is not higher than 12 Kelvin (Crapsi et al. 2007;Pagani et al. 2007). If true, this constraint would eliminate both of our evolved dust models if the "high" CR rate were accurate.
At higher densities the curves converge, and the differences between them become comparable to the uncertainty in the dust temperature. At a density of 10 6 cm −3 , Crapsi et al. (2007) find a gas temperature of around 6 K for the prestellar core L1544, and Pagani et al. (2007)  T, K no mantle Fig. 8. This shows the gas temperature as a function of position within the cloud. Each panel corresponds to a different dust coagulation model, as labelled. The black curve shows the dust temperature, as determined from the radiative transfer model. The blue and red curves show the gas temperature assuming the "high" and "low" CR ionization rates respectively from Padovani et al. (2018). Solid lines are from from the full cooling model, and dashed lines correspond to the dust-only cooling model in Ivlev et al. (2019). noting that there is some degree of uncertainty in the dust temperature, since it relies on a knowledge both of the optical properties of the absorbing dust, and also on the strength of the IR field impinging on the cloud. For these reasons, we believe the gas temperature at a density between 1 and 3 times 10 5 cm −3 has more ability to constrain dust growth and the CR rates than the temperatures at higher density.

Removal of Small Grains
The question of the effect of the grain size distribution on the ambipolar diffusion rate is important for more than just grain coagulation. At higher densities, the ionization fraction drops further, and small dust grains become very important for the coupling between the magnetic field and the neutral gas. It was shown in Zhao et al. (2016Zhao et al. ( , 2018a that removing the grains with sizes less than a few tens of nanometers could increase both the ambipolar and Hall diffusivity by a factor of 10-100. In simulations, this difference in the diffusivity results in a reduction in the magnetic flux dragged into the central disk-forming region, as well as the magnetic braking efficiency, which allows the formation of disks. Our results show that precisely these grains which are relatively well coupled to the magnetic field are removed in collisions with larger grains, thus increasing the magnetic diffusivities. Ciolek & Mouschovias (1996) noted that because of their coupling to the magnetic field, small grains will not drift to the center as fast as larger ones. This provides another reason to think that small grains would be underabundant in the core center. They predicted a reduction by a factor of ∼ 10 in the number of 10 nm or smaller grains in the center. In the absence of fragmentation, when collisions with larger grains are taken into consideration, we find a near complete elimination of such grains (see Figure 4). Ormel et al. (2009) consider a model of fragmen-tation and find that fragmentation only becomes an issue after the particles are several microns in size.
We have not explicitly taken the PAH population into account in this work, since their charging is uncertain. Due to their small size, they will be in gross violation of our assumption that the charging timescales are shorter than the stopping time and inverse gyrofrequency. Nevertheless, we believe the results shown here are sufficient to argue that they will be rapidly removed, similar to, or even faster than the smallest silicate grains. If they are predominantly charged, then they will be coupled to the magnetic field in a way similar to the smallest silicate grains, so they will be removed in collisions with larger grains at about the same rate as the smallest silicate grains we consider. If they were neutral, then their removal should be even faster, as they will be coupled to the neutral gas, and thus be moving near speed v in with respect to the predominantly negatively charged grains with size ≤ a c , which make up the majority of the surface area of the initial distribution. So, we predict that PAHs are efficiently adsorbed onto the icy mantles of larger dust grains, enriching their chemistry.

Conclusions
In this paper we discussed the role of different mechanisms that lead to grain-grain collisions. In particular, we drew attention to a mechanism not considered in previous coagulation studies -that of differential grain drift due to ambipolar diffusion. We showed that this is more efficient than turbulence in removing grains less than about 50 nm from the size distribution. This mechanism is not, however, a significant player in growing micron-sized grains. Turbulence remains the only mechanism we are aware of that can greatly change the upper edge of the size distribution. Turbulence is expected to be minimal near the core center, but may be present in the envelope. Assuming perfect sticking in collisions, and relative velocities coming from ambipolar diffusion, turbulence, gravitational settling and Brownian motion, we ran a series of coagulation simulations. We reached the following general conclusions: -The effect of grain charge on the collision cross-section due to Coulomb focussing/repulsion is not significant at densities of 10 6 and below, as there are mechanisms which lead to super-thermal grain motion. -An increased CR ionization rate increases the core lifetime against ambipolar diffusion, thus resulting in larger grains if significant turbulence is present. If the main growth mechanism is ambipolar diffusion, the differences are minimal, as the increased lifetime is largely cancelled by the correspondingly slower differential drift of dust grains. -The inclusion of icy mantles, and the recognition that the non-spherical shape of grains significantly increases their collision cross-section has a significant effect on the degree of grain growth expected within the cloud lifetime. With these effects considered, growth proceeds to 10 microns within the lifetime of the core if sonic turbulence is present at the core center. Even without turbulence, essentially all grains under 100 nanometers are eliminated, but grains larger than a micron do not form. -The resulting size distributions are nearly independent of the initial distribution of mantle thickness with grain size, but depend only on the total volume of mantle material relative to grain core, averaged over the size distribution. -Particularly at densities above 10 6 cm −3 , the removal of small grains can increase the ambipolar and Hall conductivities by over an order of magnitude. This greatly alleviates the magnetic braking problem which hinders disk formation in ideal MHD. -The gas temperature is a function both of the CR ionization rate and the dust grain size distribution. By carefully modeling the gas heating and cooling, we were able to show that the observed gas temperature near the center of dense cores such as L1544 is inconsistent with a CR ionization rate of 10 −16 s −1 if the grain growth expected from our ambipolar diffusion model has taken place, even without additional relative dust motion due to turbulence.