© K. Silsbee et al. 2020

1. Introduction

Knowledge of the dust grain-size distribution is crucial to interpretating 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 (e.g., 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⁸ cm⁻³. He ran coagulation simulations for regions with densities between 10⁵ and 10⁹ cm⁻³. 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⁵ years for densities of up to 10⁶ cm⁻³.

Ormel et al. (2009) used a detailed model for the outcome of grain–grain collisions, in which both grain size and structure are altered. 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 Könyves et al. (2015), among others. The lifetime of the central regions of a dense core (with densities above 10⁵ cm⁻³, 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⁷ years with a gas density of 10⁵ cm⁻³. In this longer time they found that large aggregates (tens of μm) form, and fragmentation becomes important.

Recent non-ideal magnetohydrodynamic simulations have pointed out the importance of the removal of very small grains for the formation of protoplanetary disks (Zhao et al. 2016, 2018a). 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 timescale for very small grain removal in prestellar cores, which represent the initial conditions in the process of star and planet formation. If this depletion timescale is significantly shorter than the dynamical timescale, 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 provide a large area on which surface chemical processes can proceed and can also 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 grains (mostly negatively charged) 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.

2. Cloud model

We take the cloud density profile that Tafalla et al. (2002) used to model the prestellar core L1544. This is a spherically symmetric cloud with the H₂ number density given by

$$\begin{array}{c}\hfill n(r)=\frac{n_{\mathrm{c}}}{1+{(r/r_0)}^{2.5}}·\end{array}$$(1)

Here n_c = 1.4 × 10⁶ cm⁻³ is the central density of the core, and r₀ = 0.014 pc = 2900 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₀, 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₀ where n = 10⁴ cm⁻³, denoted our “envelope” location. In what follows we assume that the primary ion is ${\text{H}}_{\text{3}}^{\text{+}}$ in the core center (Caselli et al. 2003) and HCO⁺ (Caselli et al. 2002) 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 Fig. 3 in Keto & Caselli (2010) we take the temperature to be 6.5 K in the core center and 10 K 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 ζ:

$$\begin{array}{c}\hfill n_i=2.1\times {10}^{-5}{\left(\frac{n_{{\mathrm{H}}_2}}{1{\mathrm{cm}}^{-3}}\right)}^{-0.56}{\left(\frac{\zeta }{{10}^{-16}{\mathrm{s}}^{-1}}\right)}^{0.5}·\end{array}$$(2)

The value of ζ 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⁻¹⁶ s⁻¹ in the core center and ζ = 3 × 10⁻¹⁶ s⁻¹ 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 ${\text{H}}_{\text{3}}^{\text{+}}$ observations (Neufeld & Wolfire 2017). A summary of the important parameters in each environment is provided in Table 1, with references provided at the bottom.

3. 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 because grains of different sizes and charge states couple differently to the neutral gas and to the magnetic and gravitational fields.

3.1. Grain charging

As we show below, the collision cross sections and the grain dynamics are both 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₂ 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 Fig. 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 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 fractions of grains that 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 a single charge except at high CR rates, where large grains can develop multiple positive charges. The charge distribution in the core center is dominated by plasma charging, and consists almost exclusively of grains with either 0 or −1 charge.

3.2. Brownian motion

We assume that each dust grain has a thermal velocity distribution, which leads to expected collision velocities

$$\begin{array}{c}\hfill \sqrt{\langle \mathrm{\Delta }{V}_{12,\mathrm{BM}}^2\rangle }=\sqrt{\frac{8k_{\mathrm{B}}T}{\pi \mu }},\end{array}$$(3)

where μ = m₁m₂/(m₁ + m₂) is the reduced mass of the two dust grains.

3.3. 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 Eqs. (5)–(10) in the limit of infinite Reynolds number. These equations specify the collision velocity as a function of the turbulent velocity v_g, and the Stokes numbers of the two grains, defined as St = τ_s/τ_L, where τ_L is the turbulent crossing time at the injection scale and τ_s is the stopping time of the grain, given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{s}}=\frac{3m}{4v_{\mathrm{th}}^{\ast }{\rho }_{\mathrm{g}}{\sigma }_{\mathrm{g}}},\end{array}$$(4)

where $v_{\mathrm{th}}^{\ast }=0.92\sqrt{4k_{\mathrm{B}}T/(\pi m_p)}$ is the thermal velocity scale of the gas (with a factor to account for the helium fraction), ρ_g = 2.8 m_pn is the gas density, and σ_g = πa² (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 cutoff at some k due to either ion-neutral friction or neutral viscosity (Xu et al. 2016). In addition, even if the cutoff is not sharp, there is evidence from simulations (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.

3.4. Systematic drifts

In this section we consider systematic drifts due to the size dependence 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 that 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 that 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

$$\begin{array}{c}\hfill c_{\mathrm{d}}g{\rho }_{\mathrm{g}}=n_i{v}_{\mathrm{in}}[{\mu }_{{\mathrm{H}}_2}{n}_{{\mathrm{H}}_2}{⟨\sigma v⟩}_{{\mathrm{H}}_2}+{\mu }_{\mathrm{He}}{n}_{\mathrm{He}}{⟨\sigma v⟩}_{\mathrm{He}}],\end{array}$$(5)

where g is the local gravitational field; μ_H2 and μ_He are the reduced masses of H₂ and helium, respectively, in collisions with the ions; and ⟨σ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\pi a_0^2\sqrt{(\alpha /a_0^3)I_H/\mu }$, where a₀ is the Bohr radius, I_H the Rydberg energy, and α is the polarizability of the neutral species. The value of $\alpha /a_0^3$ is 5.52 for H₂ and 1.39 for He. We 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 Eq. (5) (with c_d = 1) for v_in yields v_in = 5.2 × 10³ cm s⁻¹ in the core center and 4.2 × 10³ 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

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{g}}=\frac{|\overline{q}|B}{\mathit{mc}}·\end{array}$$(6)

Here $\overline{q}$ is the time-averaged grain charge, B the magnetic field strength, and c the speed of light. For the sizes where these timescales are comparable, nearly all grains are either singly negatively charged or neutral (see Table 2; we 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 ${\mathrm{\Omega }}_{\mathrm{g}}^{-1}$, are plotted in Fig. 1 as a function of grain size. We calculated τ_ch(0) and τ_ch(−1) as a function of grain size using Eqs. (3.1)–(3.5) from Draine & Sutin (1987). These equations ignore the contribution of UV photons from H₂ 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 nm, in both environments, the two charging timescales are shorter than both τ_s and ${\mathrm{\Omega }}_{\mathrm{g}}^{-1}$. 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 that the smallest grains were indeed depleted even faster, and that there was little difference for the larger grains.

Fig. 1. Timescales for processes relevant to grain dynamics: τs is the stopping time of the grain, given by Eq. (4); ${\mathrm{\Omega }}_{\mathrm{g}}^{-1}$ is the inverse grain gyrofrequency given by Eq. (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.

To estimate the grain velocities, we assume the following model. Consider a coordinate system in which $\widehat{r}$ 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 $\overline{q}$ are given by

$$\begin{array}{cc}& \frac{\mathrm{d}{v}_r}{\mathrm{d}t}=-g+{\mathrm{\Omega }}_{\mathrm{g}}{v}_{\perp }-\frac{v_r+v_{\mathrm{in}}}{{\tau }_{\mathrm{s}}},\hfill \\ \hfill & \frac{\mathrm{d}{v}_{\perp }}{\mathrm{d}t}=-{\mathrm{\Omega }}_{\mathrm{g}}{v}_r-\frac{v_{\perp }}{{\tau }_{\mathrm{s}}},\hfill \end{array}$$(7)

where v_⊥ is the velocity in the direction parallel to $\widehat{\mathit{B}}\times \widehat{\mathit{r}}$. In equilibrium, these equations can be solved to yield

$$\begin{array}{c}\hfill v_r=\frac{-g_{\mathrm{eff}}{\tau }_{\mathrm{s}}}{1+{\mathrm{\Omega }}_{\mathrm{g}}^2{\tau }_{\mathrm{s}}^2},v_{\perp }=\frac{g_{\mathrm{eff}}{\mathrm{\Omega }}_{\mathrm{g}}{\tau }_{\mathrm{s}}^2}{1+{\mathrm{\Omega }}_{\mathrm{g}}^2{\tau }_{\mathrm{s}}^2}\equiv -v_r{\mathrm{\Omega }}_{\mathrm{g}}{\tau }_{\mathrm{s}},\end{array}$$(8)

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

$$\begin{array}{c}\hfill \mathrm{\Delta }{V}_{12,\mathrm{drift}}=\sqrt{{(v_{r1}-v_{r2})}^2+{(v_{\perp 1}-v_{\perp 2})}^2}.\end{array}$$(9)

Small dust grains (those for which ${\mathrm{\Omega }}_{\mathrm{g}}^{-1}\ll {\tau }_{\mathrm{s}}$), are relatively well coupled to the magnetic field so they move with respect to the gas at velocity v_in in the $\widehat{r}$ direction. Large grains (${\mathrm{\Omega }}_{\mathrm{g}}^{-1}\gg {\tau }_{\mathrm{s}}$) do not feel the magnetic field, and drift in the $-\widehat{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

$$\begin{array}{c}\hfill a_{\mathrm{c}}=\sqrt{\frac{3|\overline{q}|B}{4\pi v_{\mathrm{th}}^{\ast }c{\rho }_{\mathrm{g}}}}·\end{array}$$(10)

We observe that the RHS of Eq. (10) also depends on a_c through $\overline{q}$. The mass of the grain does not enter Eq. (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

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}r}\left(\frac{B^2}{8\pi }\right)=-g{\rho }_{\mathrm{g}},\end{array}$$(11)

where we assume that B outside the cloud is negligible, and g(r) is calculated from the density distribution ρ_g(r) determined by Eq. (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₂ density of 3 × 10⁴ cm⁻³ in the cloud TMC-1. Multiplying by $\sqrt{3}$ to estimate the total magnetic field gives 203 μG. Our model gives 81 μG at this density.

3.5. Collision velocity from all mechanisms

We approximate the collision velocity between two dust particles as

$$\begin{array}{c}\hfill \mathrm{\Delta }{V}_{12}=\sqrt{\langle \mathrm{\Delta }{V}_{12,\mathrm{turb}}^2\rangle +\langle \mathrm{\Delta }{V}_{12,\mathrm{BM}}^2\rangle +\mathrm{\Delta }{V}_{12,\mathrm{drift}}^2}.\end{array}$$(12)

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 a distribution of collision velocities. Given the large uncertainties on the individual sources of relative motion, both these issues seem rather minor.

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 highest 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 where 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 gray and black lines across the panels, corresponding to same-size collision partners where the systematic drifts have no effect.

Fig. 2. Collision velocities calculated from Eq. (12) in the space of the sizes of the collision partners. The turbulent velocity vg and the constant cd, which describes the speed of the ion-neutral drift (see Eq. (5)), are different in each panel, as labeled. The isothermal sound speed cs is given by $\sqrt{k_{\mathrm{B}}T/(2.33m_p)}$. The labeled contours are at 1, 10, and 50 m per sec.

3.6. 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 shell as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{coll}}}{\mathrm{d}t}=\frac{N_i{N}_j}{V}{\displaystyle \sum _{q_i,q_j}{v}_{\mathit{ij}}(q_i,q_j){\sigma }_{\mathit{ij}}(q_i,q_j)f_{q_i}{f}_{q_j}.}\end{array}$$(13)

Here N_i and N_j are the numbers of particles in mass bin i and j, respectively, and v_ij(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. We note that for our standard assumption that the charging timescales are rapid compared to the other dynamical timescales, v_ij is calculated using ${\overline{q}}_i$ and ${\overline{q}}_j$; however, σ_ij is still calculated using the instantaneous charge states q_i and q_j. The parameter f_qi is the fraction of particles in mass bin i with charge q_i, and f_qj is defined analogously for mass bin j; V is the volume of the shell. The collision cross section is given by

$$\begin{array}{c}\hfill {\sigma }_{\mathit{ij}}(q_i,q_j)=\pi {(a_i+a_j)}^{2}max\{0,1-\frac{2q_i{q}_j}{\mu v_{\mathit{ij}}^2(a_i+a_j)}\},\end{array}$$(14)

where a_i and a_j are the radii of the two dust grains. We find in Sect. 5.3 that Coulomb attraction and repulsion play a negligible role in the environments we consider¹. This is to be expected given the typical collision velocities. As can be seen from Eq. (14), the critical parameter determining the Coulomb effect is $2q_i{q}_j/\left(\mu v_{ij}^2(a_i+a_j)\right)$. For 10 nm grains, this is less than unity for v >  2 m s⁻¹. It scales inversely with the fourth 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⁻¹, except in the case with no turbulence or ambipolar diffusion.

4. 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 τ_cl, the lifetime of the cloud against ambipolar diffusion, which we take as the runtime of our simulations.

4.1. Cloud lifetime

We assume that the lifetime of the cloud at scale r is roughly equal to the ambipolar diffusion time at that scale in the cloud (Tassis & Mouschovias 2004), i.e., the cloud lifetime τ_cl is given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{cl}}=r/v_{\mathrm{in}}.\end{array}$$(15)

This works out to be 2.6 × 10⁵ years in the core center, and 2.3 × 10⁶ years in the envelope, assuming the standard value for the ionization fraction given in Eq. (2).

4.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 distribution (Mathis et al. 1977) 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

$$\begin{array}{c}\hfill n_{\mathrm{d}}(a)=\frac{3{\rho }_{\mathrm{g}}{f}_{\mathrm{d}}{a}^{-7/2}}{8\pi {\rho }_{\mathrm{d}}\sqrt{a_{\max }}},\end{array}$$(16)

where f_d = 0.01 is the gas-to-dust mass ratio. Then, we consider a small grain of size a_s <  a_c. We 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 approximate the collision cross section as $\pi a_{\text{l}}^2$. Then the expected time for the smaller grain to experience a collision with a larger one is given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{AD}}^{-1}=v_{\mathrm{in}}{\displaystyle {\int }_{a_{\mathrm{c}}}^{a_{\max }}\pi a_{\mathrm{l}}^2{n}_{\mathrm{d}}(a_{\mathrm{l}})\mathrm{d}{a}_{\mathrm{l}}.}\end{array}$$(17)

This gives

$$\begin{array}{c}\hfill {\tau }_{\mathrm{AD}}=\frac{4{\rho }_{\mathrm{d}}\sqrt{a_{\max }{a}_{\mathrm{c}}}}{3{\rho }_g{f}_{\mathrm{d}}{v}_{\mathrm{in}}}·\end{array}$$(18)

For the parameters used here, we find τ_AD = 1.1 × 10⁴ years in the core center, and 2.4 × 10⁶ years in the envelope.

4.3. 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 Eq. (8) (Ω_g → 0) yields a collision speed between two particles of sizes a_s and a_l of

$$\begin{array}{c}\hfill v_{\mathrm{grav}}={\rho }_{\mathrm{d}}g(a_{\mathrm{l}}-a_{\mathrm{s}})/(v_{\mathrm{th}}^{\ast }{\rho }_{\mathrm{g}}).\end{array}$$(19)

Then, assuming the collision cross section to be $\pi a_{\text{l}}^2$, dropping the term a_s in Eq. (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:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{grav}}=M/\dot{M}=\frac{4v_{\mathrm{th}}^{\ast }}{3g{f}_{\mathrm{d}}}·\end{array}$$(20)

This expression is equal to 1.9 × 10⁶ years in the core center and 8.4 × 10⁶ years in the envelope.

4.4. Coagulation of large grains due to turbulence

In the limit of small Stokes numbers, Eqs. (6), (7), and (10) from Ormel & Cuzzi (2007) show that the relative velocity between two grains of the same size is

$$\begin{array}{c}\hfill \mathrm{\Delta }{V}_{12,\mathrm{drift}}^2=v_{\mathrm{g}}\sqrt{\psi (ϵ){\mathrm{St}}_{\mathrm{l}}},\end{array}$$(21)

where St_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\pi a_{\text{l}}^2$, and calculate the velocities using Eq. (21) with ψ = 1.95. These approximations yield a timescale for growth of the largest grains from collisions arising from turbulence:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{turb}}(a_{\mathrm{l}})={\left(216{\psi }^2{\pi }^2\right)}^{-1/4}\sqrt{\frac{{\rho }_{\mathrm{d}}{v}_{\mathrm{th}}{a}_{l}r}{f_{\mathrm{d}}^{\ast 2}{\rho }_{\mathrm{g}}{v}_{\mathrm{g}}^3}}·\end{array}$$(22)

Here r enters through the assumption that it is equal to the turbulence injection scale. Assuming v_g = c_s, where $c_{\mathrm{s}}=\sqrt{k_{\mathrm{B}}T/(2.33m_p)}$, this expression gives 9.4 × 10⁴ years in the core center, and 2.0 × 10⁶ 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 Eq. (22) in conjunction with the fact that a/ȧ = 3m/ṁ, we can solve for a_l(t), assuming that a_l(t₀) = a_max:

$$\begin{array}{c}\hfill a_{\mathrm{l}}(t)=a_{\max }{[1+\frac{t}{6{\tau }_{\mathrm{turb}}(a_{\max })}]}^2.\end{array}$$(23)

We show in Sect. 5.3 that this is a reasonable approximation, but yields a slight overestimate.

5. Results

5.1. Core center

We consider here the coagulation of dust at radius r₀ 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 Fig. 4. The values of v_g and c_d in each panel are equal to those in the corresponding panel of Fig. 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 Sect. 4.1), equal in this case to 2.6 × 10⁵ years. In all cases, for times shorter 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 Fig. 3 that collisions between bodies with sizes around 1 μm depend primarily on the strength of the turbulence. In Sect. 4 we derived growth timescales of 9.4 × 10⁴ years and 1.9 × 10⁶ 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. It should be noted from Eq. (22) that ${\tau }_{\mathrm{turb}}\propto v_{\mathrm{g}}^{-3/2}$. Examination of the bottom panels of Fig. 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.

Fig. 3. Dominant source of collision velocity in the space of the sizes of the collision partners. The turbulent velocity vg and the constant cd, which describes the speed of the ion-neutral drift (see Eq. (5)), are varied between the panels, as labeled. vg and cd in each panel take the same values as in the corresponding panel of Fig. 2. The mapping between color and dominant velocity source is shown in the legend on top.

We find when c_d = 1 that there is a rapid depletion on ∼10⁴ year timescales of grains less than a_c (defined by Eq. 10), consistent with the estimate in Eq. (18). The bottom right panel of Fig. 4, shows that Brownian motion has a minor effect, and only for the very smallest grains.

Fig. 4. Evolution of the grain-size distribution in the core center for different values of the parameters vg and cd, as labeled in the panels. The different color curves correspond to different amounts of time the simulation was run for in terms of the cloud lifetime τcl (see Eq. (15)). Here τcl = 2.6 × 105 years. In each case the initial distribution is an MRN distribution from 5 to 250 nm, which is shown as the black dashed line.

5.2. Envelope

We also consider coagulation in a lower density environment (n_H2 = 10⁴ cm⁻³). This corresponds to a location at r = 20 800 AU in the model of Tafalla et al. (2002). Figure 5 is analogous to Fig. 4, but made for the envelope. As discussed in Sect. 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. The value of τ_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.

Fig. 5. Same as Fig. 4, but for the envelope. The colored curves show the size distribution at different times. In this case τcl = 2.3 × 106 years. The black solid line shows the curve inferred in Weingartner & Draine (2001) for dense regions.

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⁵b_c = 3.0). It is not clear which region of the cloud these observations are 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⁴ cm⁻³ or that our estimate of the cloud lifetime is too short (see next subsection).

5.3. Effect of CR ionization rate

Cosmic rays 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 Eq. (5). Since the ionization fraction is assumed to be proportional to $\sqrt{\zeta }$, this means that $v_{\mathrm{in}}\propto 1/\sqrt{\zeta }$. As the cloud lifetime τ_cl is proportional to 1/v_in (see Sect. 4.1), there is more time for coagulation to occur at higher ζ. We note, however, that τ_AD also scales as 1/v_in (see Sect. 4.2). Hence, to the 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 Sect. 3.1, CRs also excite electronic states of H₂ that decay, producing UV photons that 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 Fig. 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 Eq. (15). Because v_in ∼ ζ^−1/2, the cloud lifetime is proportional to ${\zeta }^{\frac{1}{2}}$.

Fig. 6. Top panel: size distribution in the envelope for the case with sonic turbulence and cd = 0. The different curves correspond to different amounts of time and different CR rates (see legend at top). Time 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 cd = 0 (top panel) the curves at different times are nearly independent of the ionization rate and the curves overlap nearly perfectly (i.e., only one curve is visible for each time). Middle panel: without turbulence, but assuming cd = 1. Bottom panel: shown in blue is the peak of the distribution as a function of time (assuming ζ = ζP18, cd = 0, and vg = cs), and in orange the approximation given by Eq. (23).

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 long cloud lifetime in our model. This shows that the issue of Coulomb focussing–repulsion 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 we 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. Instead, we see that over the lifetime of the cloud, the high ζ model still gets rid of more of the small grains 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 Eq. (10)). Equations (15) and (18) demonstrate that ${\tau }_{\mathrm{AD}}/{\tau }_{\mathrm{cl}}\propto \sqrt{a_{\mathrm{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²dN/dm) is shown by the blue curve in the bottom panel of Fig. 6, and the approximation given by Eq. (23) is shown in orange. The orange curve gives a slight overestimate for two reasons. First, there is initially no movement of the peak due to the sharp cutoff of the original size distribution. Additionally, Eq. (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.

6. Additional effects

6.1. Evolution in the core center and the 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 simulations as the initial input to the core center simulations. We found, however, that the preprocessing 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 the 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.

6.2. Contribution of small grains to ambipolar diffusion velocity

In Eq. (5) we calculate the ambipolar drift speed, assuming that neutrals are slowed down only through collisions with ions. Neutrals must to some extent be slowed by collisions with small magnetically coupled dust grains as well. To calculate this, we can add the term

$$\begin{array}{c}\hfill {\displaystyle {\int }_{a_{\min }}^{a_{\max }}\pi a^{2}n(a)v_{\mathrm{dn}}(a)[m_{{\mathrm{H}}_2}{n}_{{\mathrm{H}}_2}⟨v_{{\mathrm{H}}_2}⟩+m_{\mathrm{He}}{n}_{\mathrm{He}}⟨v_{\mathrm{He}}⟩]\mathrm{d}a}\end{array}$$(24)

to the right-hand side of Eq. (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 Eq. (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 Eq. (24) to the right-hand side of Eq. (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, Fig. 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.

6.3. Icy mantles and non-compact grains

In the interests of simplicity, the preceding analysis leaves out two important effects that may enhance dust coagulation. The first is that at densities higher than ∼10⁵ cm³, 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 in 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 g per cm³ 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.

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 Fig. 4, but as discussed above includes the icy mantles and considers the effect of non-spherical and/or 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 = 13 nm, 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 most of the volume being in grains of size between 1 and 2 times the mantle thickness, as can be seen in Fig. 7. Constant mantle thickness is the result we 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. 2018b). 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.

Fig. 7. Same as Fig. 4, but assuming that grains are covered with icy mantles, and that they are non-spherical and/or porous (which increases their projected surface area by a factor of 2). As in Fig. 4, the panels show different degrees of turbulence and ambipolar diffusion, parameterized by vg and cd, as labeled. The sharp peak around 20 nm is the result of the 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.

In determining the grain dynamics, we assume that for grains of size between a_min and a_max,

$$\begin{array}{c}\hfill \rho ={\rho }_{\mathrm{core}}{(1-\frac{a_{\mathrm{mantle}}}{a})}^3+{\rho }_{\mathrm{mantle}}[1-{(1-\frac{a_{\mathrm{mantle}}}{a})}^3].\end{array}$$(25)

For sizes a >  a_max the density smoothly approaches the mean density given by

$$\begin{array}{c}\hfill \overline{\rho }=\frac{{\rho }_{\mathrm{core}}+V_{\mathrm{rel}}{\rho }_{\mathrm{mantle}}}{1+V_{\mathrm{rel}}}·\end{array}$$(26)

Comparing Fig. 7 with Fig. 4, it is apparent that growth to above one micron within the estimated core lifetime is only possible if we consider the icy mantles and non-spherical and/or porous grain shapes. With these factors taken into consideration, we find that growth to about ten microns occurs in our models with sonic turbulence. With just ambipolar diffusion, the differences between Figs. 4 and 7 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.

6.4. 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⁻¹. 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⁻¹ (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 the grains 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₂O spheres of 9.6 m s⁻¹. 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⁻¹. Not all refractory materials, however, have such a high sticking threshold. Musiolik et al. (2016) found a critical velocity of 0.04 m s⁻¹ for 100 μm CO₂ grains, and noted that it would be nearly an order of magnitude lower if these grains were 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.

7. Some implications of this work

7.1. Gas temperature

Ivlev et al. (2019) has shown that in the core center the gas temperature is sensitive to the dust grain-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 is 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) and Sipilä et al. (2019), which combines hydrodynamics with chemical and radiative transfer models in a self-consistent 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 labeled 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 Sect. 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 labeled “High CR” and “Low CR” correspond to the high and low models in Padovani et al. (2018), which have respective ionization rates of around 10⁻¹⁶ s⁻¹ and 2 × 10⁻¹⁷ s⁻¹ 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 Eq. (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⁵ cm⁻³ the temperature is not higher than 12 K (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.

Fig. 8. Gas temperature as a function of position within the cloud. Each panel corresponds to a different dust coagulation model, as labeled. 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). The solid lines are the full cooling model, and the dashed lines correspond to the dust-only cooling model in Ivlev et al. (2019).

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⁶ cm⁻³, Crapsi et al. (2007) find a gas temperature of around 6 K for the prestellar core L1544, and Pagani et al. (2007) estimate 7 K for the prestellar core L183. Our model dust temperature at this density is 7.1 K. The warmest gas temperature of any of our models is only 8.1 K. It is worth noting that there is some degree of uncertainty in the dust temperature since it relies on a knowledge of the optical properties of the absorbing dust and 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⁵ cm⁻³ has more ability to constrain dust growth and the CR rates than the temperatures at higher density.

7.2. 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 furthexr, 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. (2016, 2018b) that removing the grains with sizes of 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 the grains that 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 Fig. 4). Ormel et al. (2009) consider a model of fragmentation and find that fragmentation only becomes an issue after the particles are several microns in size.

We do not explicitly take the polycyclic aromatic hydrocarbon 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 would 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 most of the surface area of the initial distribution. So, we predict that polycyclic aromatic hydrocarbons are efficiently adsorbed onto the icy mantles of larger dust grains, enriching their chemistry.

8. 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 smaller 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 bound 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⁶ and below as there are mechanisms that lead to superthermal 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 and their porosity 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⁶ cm⁻³, 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 magnetohydrodynamics.

• 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⁻¹⁶ s⁻¹ if the grain growth expected from our ambipolar diffusion model has taken place, even without additional relative dust motion due to turbulence.

Acknowledgments

The authors wish to thank the referee, Satoshi Okuzumi, for his careful reading of the manuscript and insightful comments that considerably improved the paper. K.S. would like to thank Vincent Guillet for stimulating discussions. A.V.I. acknowledges support by the Russian Science Foundation (project 18-12- 00351).

References

All Tables

All Figures

Fig. 1. Timescales for processes relevant to grain dynamics: τs is the stopping time of the grain, given by Eq. (4); ${\mathrm{\Omega }}_{\mathrm{g}}^{-1}$ is the inverse grain gyrofrequency given by Eq. (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.

In the text

	Fig. 2. Collision velocities calculated from Eq. (12) in the space of the sizes of the collision partners. The turbulent velocity vg and the constant cd, which describes the speed of the ion-neutral drift (see Eq. (5)), are different in each panel, as labeled. The isothermal sound speed cs is given by $\sqrt{k_{\mathrm{B}}T/(2.33m_p)}$. The labeled contours are at 1, 10, and 50 m per sec.
In the text

Fig. 3. Dominant source of collision velocity in the space of the sizes of the collision partners. The turbulent velocity vg and the constant cd, which describes the speed of the ion-neutral drift (see Eq. (5)), are varied between the panels, as labeled. vg and cd in each panel take the same values as in the corresponding panel of Fig. 2. The mapping between color and dominant velocity source is shown in the legend on top.

In the text

Fig. 4. Evolution of the grain-size distribution in the core center for different values of the parameters vg and cd, as labeled in the panels. The different color curves correspond to different amounts of time the simulation was run for in terms of the cloud lifetime τcl (see Eq. (15)). Here τcl = 2.6 × 105 years. In each case the initial distribution is an MRN distribution from 5 to 250 nm, which is shown as the black dashed line.

In the text

	Fig. 5. Same as Fig. 4, but for the envelope. The colored curves show the size distribution at different times. In this case τcl = 2.3 × 106 years. The black solid line shows the curve inferred in Weingartner & Draine (2001) for dense regions.
In the text

Fig. 6. Top panel: size distribution in the envelope for the case with sonic turbulence and cd = 0. The different curves correspond to different amounts of time and different CR rates (see legend at top). Time 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 cd = 0 (top panel) the curves at different times are nearly independent of the ionization rate and the curves overlap nearly perfectly (i.e., only one curve is visible for each time). Middle panel: without turbulence, but assuming cd = 1. Bottom panel: shown in blue is the peak of the distribution as a function of time (assuming ζ = ζP18, cd = 0, and vg = cs), and in orange the approximation given by Eq. (23).

In the text

Fig. 7. Same as Fig. 4, but assuming that grains are covered with icy mantles, and that they are non-spherical and/or porous (which increases their projected surface area by a factor of 2). As in Fig. 4, the panels show different degrees of turbulence and ambipolar diffusion, parameterized by vg and cd, as labeled. The sharp peak around 20 nm is the result of the 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.

In the text

Fig. 8. Gas temperature as a function of position within the cloud. Each panel corresponds to a different dust coagulation model, as labeled. 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). The solid lines are the full cooling model, and the dashed lines correspond to the dust-only cooling model in Ivlev et al. (2019).

In the text