© V. Guillet et al. 2020

1 Introduction

While it is now well established that stars form through the collapse of prestellar cores, understanding the exact outcome of this process remains a challenge (Li et al. 2014). In particular, many studies have found that the properties of the centrifugally supported disks that form around protostars sensitively depend on the intensity of the magnetic field and even possibly on its orientation (Allen et al. 2003; Mellon & Li 2008; Hennebelle & Fromang 2008; Joos et al. 2012; Li et al. 2013; Gray et al. 2018). This is due to the magnetic braking that can efficiently transport angular momentum from the inner part of the collapsing cloud to the surrounding envelope. In the most extreme case, it has even been found that the formation of a centrifugally supported disk can be entirely suppressed, a process known as catastrophic magnetic braking (Allen et al. 2003; Mellon & Li 2008; Hennebelle & Fromang 2008).

The magnetic field evolution is a direct consequence of its coupling with the gas. If this coupling is perfect (ideal Magnetohydrodynamics; MHD), then magnetic intensity is typically expected to be ∝ ρ^κ, where ρ is the gas density and κ ≃ 1∕2−2∕3. On the other hand, in the extreme case where the magnetic field would be completely decoupled from the gas, the magnetic intensity would stay constant and the magnetic field would have no (or a much more limited) influence on the gas evolution. It is therefore fundamental to understand with enough accuracy how magnetic field and gas are coupled together. Since the gas within molecular clouds is weakly ionized, with a ionisation fraction on the order of, or even below 10⁻⁷ (e.g. Shu et al. 1987), the coupling between the gas, mainly the neutrals, and the magnetic field is imperfect. The neutrals can slip through the field lines which are attached to the ions, a process known as ambipolar diffusion. This process largely dominates over the other non-ideal MHD processes in the interstellar medium (ISM) such as the Hall effect and the Ohmic resistivity, which are due to the imperfect coupling between the magnetic field and the ions and the electrons, respectively.

In the context of dense prestellar cores, it is believed that dust grains are playing an important role. As the recombination rate of ions on the surface of grains scales with the density, the ionisation is several orders of magnitude lower in collapsing cores than in the rest of the ISM, and the impact of non-ideal MHD processes is more important. Many detailed calculations of the ionisation inside dense cores have been performed (e.g. Nakano et al. 2002; Kunz & Mouschovias 2009; Zhao et al. 2016; Wurster et al. 2016; Marchand et al. 2016; Dzyurkevich et al. 2017). The resulting resistivities depend on the exact assumptions regarding the ionisation rate, the chemistry network and the grain properties. In particular, it has been concluded that the abundance of small grains is particularly critical (Nishi et al. 1991; Zhao et al. 2016, 2018; Dzyurkevich et al. 2017; Grassi et al. 2019). For example Zhao et al. (2016) assuming an MRN size distribution (Mathis et al. 1977)found that the ambipolar and Hall resistivities reach a maximum when only grains of size above 0.1 and 0.04 μm are considered, respectively. Moreover when varying the size of the smallest grains considered from 0.005 to 1 μm, Zhao et al. (2016) found that the resistivities vary by one to two orders of magnitude (depending on the considered gas density). Therefore, asufficient knowledge of the grain distribution appears to be a crucial issue. Indeed several studies have investigated the impact that the resistivities have on circumstellar disks confirming that disks form more easily when the resistivity is high which happens when small grains are removed (e.g. Krasnopolsky et al. 2011; Zhao et al. 2016). Hennebelle et al. (2016) proposed an analytical model in which the size of the disk at the early stage is predicted to be $\propto {\eta }_{\text{AD}}^{2/9}$, where η_AD is the ambipolar resistivity. The influence of the Hall resistivity has also been investigated in a series of papers (Krasnopolsky et al. 2011; Wurster et al. 2016; Tsukamoto et al. 2015; Koga et al. 2019). Because of the quadratic dependence on the magnetic field of the Hall term, changing the sign of B does not lead to a physically identical situation and one must distinguish between the parallel and anti-parallel configurations depending on the respective orientation of B and the angular momentum. In the anti-parallel configuration the Hall effect tends to produce bigger disks while in the parallel one, it tends to produce smaller ones.

Knowing the grain distribution within dense cores therefore appears to be a major issue. As recalled above, many authors have assumed an MRN type distribution, with a maximal grain size ~ 0.25 μm. However, this particular size distribution is only valid for the diffuse phase of the ISM. In denser gas, the phenomenon of coreshine in the NIR (Pagani et al. 2010) and the observed increase in FIR and sub-mm emissivity (e.g. Stepnik et al. 2003; del Burgo & Laureijs 2005; Ysard et al. 2013) all point towards the coagulation of grains in the envelope of molecular clouds, affecting the grain structure and composition (Jones et al. 2013; Köhler et al. 2015) and increasing the grain size significantly (~ 1 μm, Steinacker et al. 2015) or only modestly (<0.5 μm, Ysard et al. 2016), depending on the dust models used. This indicates that dust coagulation must be an efficient process and that, given the grain size-dependence of the resistivities inferred in previous studies, it may play an major role in the coalescence of the cloud.

In this article, we make use of the Paris-Durham shock code (Flower & Pineau des Forêts 2003) amended for the calculation of dust charging, dynamics and evolution (Guillet et al. 2007, 2011), to study the coagulation of dust grains in a collapsing core, as well as its feedback on the evolution of the Ohmic, Hall and ambipolar diffusion MHD resistivities of the medium. The article is structured as follows. In Sect. 2, we present our model for the grain dynamics and coagulation in a protostellar collapse. Section 3 deals with the calculation of grain charge, ionisation equilibrium and MHD resistivities. In Sect. 4, we present our results for a collapse without dust, a collapse with dust but no coagulation and a collapse with dust and coagulation. We discuss in Sect. 5 some limitations in our study and recall our main results in Sect. 6.

2 The dust evolution model

The Paris-Durham shock code is a fortran code initially designed to solve the structure of stationary shocks and predict the intensity of their emission lines (Flower et al. 1985). From the beginning, it included more than 100 chemical species and a large network of chemical reactions (>1000). The impact of dust grains on the chemical compositionand dynamics of the shock was later included (Flower & Pineau des Forêts 2003), but only through a single, effective, grain fluid. In Flower et al. (2005), the code was adapted to study grain–grain coagulation and its feedback on the ionisation and depletion during the gravitational collapse of an elementary cell composed of gas and dust.

The Dust Dynamics and Processing code (DUSTDaP) was developed from the Paris-Durham shock code to study the detailed charge and dynamics of a full size distribution of dust grains in shocks (Guillet et al. 2007). It includes the physics of grain charging as well as grain–grain destruction, namely the shattering and vaporisation of grain cores, encountered in shocks propagating through dense clouds (Guillet et al. 2009, 2011). We detail below how this code was amended to follow the evolution of the dust size distribution by grain–grain coagulation during the isothermal collapse of 1 cm⁻³ of gas.

2.1 Initial conditions in the parent cloud

Our gas-phase chemistry derived from Flower et al. (1985) and Le Gal et al. (2014) includes 134 gas-phase species and ~ 700 chemical reactions, with ion-neutral, neutral-neutral, and dissociative recombination reactions involving species containing H, He, C, N, O and S. The initial distribution of the elemental abundances across the gas phase and the solid phase is given in Table 1 of Flower et al. (2005) with the exception of S and Fe, whose fractional elemental abundances were 1.47 × 10⁻⁵ (Hily-Blant et al. 2020) and 1.50 × 10⁻⁸, respectively.

In dense clouds, grains are covered by icy mantles (Hagen et al. 1983). The thickness of this mantle is independent of the grain size and can be computed knowing the size distribution of grain cores (see Appendix B from Guillet et al. 2007). For an MRN size distribution (α = −3.5) with core radii ranging over [5 : 250] nm, the mantle is 8.8 nm thick using the Table 2 from Flower & Pineau des Forêts (2003) for the abundances of chemical species in icy mantles, assuming a specific density of 1 g cm⁻³ for ices and 3 g cm⁻³ for cores. The initial dust-to-gas mass ratio is initially of 0.9%, and will increase up to 1.2% through accretion during the collapse.

The initial proton density is n_H,0 = 10⁴ cm⁻³ and the cosmic ray ionisation rate ζ = 5 × 10⁻¹⁷ s⁻¹. Under these conditions, we determine the steady-state abundances of the chemical species in the gas phase (neglecting any further depletion) which are then used as initial conditions of the collapsing core.

2.2 Recipe for core collapse

The free-fall timescale of a spherical cloud of uniform mass density ρ₀ is: $${\tau }_{\text{ff}}=\sqrt{\frac{3\pi }{32\mathcal{G}{\rho }_0}},$$(1)

where $\mathcal{G}$ is the gravitational constant. For our initial proton density (n_H,0 = 10⁴ cm⁻³) τ_ff ≃430 kyr. We will consider the collapse of a spherical core in a free-fall time.

If we state x = R(t)∕R_C the ratio of the core radius at time t to its initial radius (0 ≤ x < 1), the dynamics of the free-fall follows the equation (Flower et al. 2005): $$\frac{\text{d}x}{\text{d}t}=-\frac{\pi }{2{\tau }_{\text{ff}}}\sqrt{\frac{1}{x}-1}.$$(2)

We must now make an assumption on the dependence of the gas density with x. In Appendix A, we demonstrate that the assumption of an uniform compression of all fluids as per Flower et al. (2005), meaning that the mass density $\rho (t)={\rho }_0{\left(R_{\text{C}}/R(t)\right)}^3$ is a function of time only, is a good approximation of what happens in the Larson compression scenario (Larson 1969) regarding the level of coagulation achieved up to a given local density, even if it is a bad description for the density profile of the collapsing core expected at a given time. Therefore, we state $$\frac{1}{n_{\text{H}}}\frac{\text{d}{n}_{\text{H}}}{\text{d}t}=\frac{1}{\rho }\frac{\text{d}\rho }{\text{d}t}=-\frac{3}{x}\frac{\text{d}x}{\text{d}t}.$$(3)

Following Marchand et al. (2016), the magnetic field is assumed to be transverse to the direction of the collapse, and its intensity (inGauss) is assumed to increase with the gas density during the collapse as follows $$B(\text{G})=b\sqrt{n_{\text{H}}/n_{\text{H},0}},$$(4)

with b = 30 μG taken as a reference value at the initial density n_H,0. This corresponds to a mass-to-flux ratio of 2.3 (Li et al. 2011). For clarity, the evolution with time of the cloud density and magnetic field intensity are plotted in Fig. 1.

Fig. 1 Evolution of the cloud proton density nH and magnetic field intensity B with timein our model. The collapse time is the free-fall time computed at the initial cloud density nH,0 = 104 cm−3: τff ≃ 430 × 103 yr.

2.3 Grain–grain coagulation

In order to study the evolution of the grain size distribution during cloud core collapse, we updated the DUSTDaP code to include the coagulation of grains in low velocity grain–grain collisions. The evolution of the grain size distribution by grain–grain coagulation is controlled by the Smoluchowski equation (Smoluchowski 1916; Mizuno et al. 1988): $$\begin{array}{l}\frac{\text{d}\rho (m,t)}{\text{d}t}=-{\displaystyle {\int }_0^{\infty }}mK(m,m^{\prime })n(m,t)n(m^{\prime },t)\text{d}{m}^{\prime }\\ +\frac{1}{2}{\displaystyle {\int }_0^{m}m}K(m^{\prime }-m,m^{\prime })n(m^{\prime }-m,t)n(m^{\prime },t)\text{d}{m}^{\prime },\end{array}$$(5)

where m is the grain mass, n(m, t) (respectively ρ(m, t)) the density (respectively mass density) of grains of mass m at time t, and K(m, m′) (cm ³ s⁻¹) the collision kernel between grains of mass m and m′.

Our numerical model for grain–grain coagulation is rudimentary: grains are assumed to be compact and spherical, and the product of grain–grain coagulation as well (Hirashita 2012). The dust size distribution is modelled by N_bins logarithmic bins (Guillet et al. 2007). For the code to run fast, the number of bins cannot exceed a few tens because of the numerous integration variables associated to the grain charge distribution for each bin. In Appendix B, we detail the implementation and tests of our coagulation algorithm, and demonstrate that the convergence of the algorithm is obtained as soon as there is more than ~7−8 logarithmic bins per decade in grain radius, or equivalently 13 for the MRN size distribution.

2.4 Handling of ice mantles covering grain cores

The presence of icy mantles on the surface of grains is very important for grain–grain coagulation. Firstly, ice mantles greatly improve the sticking probability in grain–grain collisions (Chokshi et al. 1993). Secondly, this mantle significantly enhances the radius of the smaller grains (Guillet et al. 2007), and therefore their coagulation rate.

Coagulation will modify the uniformity of mantle thickness by transferring the large volume of ices carried by small grains to the larger grains. The mass of icy mantles in each bin is therefore a variable that must be integrated. The transfer of the mass of icy mantles through coagulation is handled together with the transfer of the mass of grain cores. At each step, knowing the mass of icy mantles and the properties of grain cores (average radius, cross-section and mass, number density) in each bin, we are able to compute the mantle thickness of grains using the equation detailed in Appendix B from Guillet et al. (2007), and derive the effective radius, cross-section and mass of the core-mantle grains in each size bin.

Once the collapse has started, ices continue to accrete onto the surface of grains, increasing the mantle thickness uniformly. We note that mantle accretion cannot modify the number of grains.

2.5 Velocity limit for grain–grain coagulation

Laboratory experiments (Poppe & Blum 1997) and theoretical studies (Chokshi et al. 1993) led to a first conclusion that grain coagulate if their relative velocity does not exceed a velocity limit, and bounce on each other otherwise. This velocity limit for coagulation is a function of the grain size (Chokshi et al. 1993) $$\mathcal{V}=V_{0.1\mu \text{m}}{\left(\frac{\widehat{a}}{0.1\mu \text{m}}\right)}^{-5/6},$$(6)

where â = (a₁ × a₂)∕(a₁ + a₂) is the reduced grain radius of the two colliding grains, and V_{0.1 μm} the velocity limit for two grains of 0.1 μm, which lies between 0.1 and 0.4 km s⁻¹ (Chokshi et al. 1993; Poppe & Blum 1997). We will take the upper value V_{0.1 μm} = 0.4 km s⁻¹ as per Poppe & Blum (1997), to take into account the increase in the velocity limit due to the presence of icy mantles on the surface of grain cores.

2.6 Grain dynamics

Grain dynamics is the fundamental ingredient of grain growth as it controls the collisional rates between grains (e.g. Ossenkopf 1993). In the context of core collapse, grain velocity dispersion can originate from the thermal agitation of the grain, from the acceleration of grains by HD or MHD turbulence, from the gravitational force, or from the electromagnetic force.

2.6.1 Grain acceleration by Hydrodynamic turbulence

Since the pioneering work by Voelk et al. (1980), the acceleration of dust grains by HD turbulence has been modelled by the stochastic acceleration of grains in a Kolmogorov cascade of turbulent eddies. In this framework, grains are accelerated by the eddies to which they can couple.

We use the model of grain acceleration by HD turbulence, as revised by Ormel & Cuzzi (2007). The injection scale of the turbulence is assumed to be equal to the Jeans length, and the injection velocity to the isothermal sound speed $c_{\text{s}}=\sqrt{k_{\text{B}}T/\mu }$, with μ the mean mass of gas particles. This defines an upper timescale for the turbulence cascade $${\tau }_{\text{L}}=\frac{L_{\text{Jeans}}}{c_{\text{s}}}=\frac{1}{2}\sqrt{\frac{\pi }{\mathcal{G}\rho }}.$$(7)

The dissipation scale of the turbulence is controlled by the Reynolds number (Ormel et al. 2009) $$\text{Re}=6.2\times {10}^7\sqrt{\frac{\rho /\mu }{{10}^5{\text{cm}}^{-3}}},$$(8)

corresponding to a timescale ${\tau }_{\eta }={\tau }_{\text{L}}/\sqrt{\text{Re}}$.

We chose not to use the model of grain acceleration by MHD turbulence of Yan et al. (2004) because the relations established in this article are mostly adapted to diffuse and moderately dense medium (n_H ≤ 10⁴ cm⁻³) and not the high densities encountered in our numerical collapsing core (10⁴ ≤ n_H ≤ 10¹² cm⁻³). In Sect. 5, we discuss the limitations of the HD models of grain dynamics and how the use of a model of grain acceleration by MHD turbulence would affect our conclusions.

In the theory of grain acceleration by hydrodynamical turbulence, the grain dynamic is controlled by the stopping time τ_drag of the grain by collisions with gas particles. Large grains tend to couple to the large eddies, thereby acquiring strong kicks that accelerate them through a random walk in the velocity space. In the Epstein regime valid for sub-sonic grain velocities, the grain stopping time through collisions with the gas is $$\tau =\frac{3m}{4\rho v_{\text{th}}\widehat{\sigma }}\propto a,$$(9)

where a, $\widehat{\sigma }$ and m are the grain radius, cross-section and mass, respectively, ρ is the mass density of the gas and $v_{\text{th}}=\sqrt{8k_{\text{B}}T/\pi \mu }$ the mean velocity of gas particles¹.

We now consider in detail the expressions for the relative velocity between grains of different sizes. We call target (stopping time τ_T) the larger grain, and projectile (stopping time τ_P) the smaller grain of the two. Three different regimes can be distinguished for the rms relative velocity ΔV_P,T between the projectile and the target (Ormel & Cuzzi 2007): (1) a regime of tight coupling where the stopping timescale of the target τ_T is smaller than the dissipation timescale of the turbulence τ_η, meaning that the target is so small that it cannot escape even from the smallest eddies; (2) a regime for heavy particles where the target is so large that its stopping timescale τ_T is much longer than injection timescale of the turbulence τ_L, thereby reducing the efficiency of the kicks by the largest eddies as the grain size increase; and (3) an intermediate regime where grains are optimally accelerated by an eddy of a particular size. Regarding the grains sizes involved in our study, small grains are in the tightly coupled particles regime, while large grains are in the intermediate regime. No grains from our size distributions fall into the heavy particles regime at any density as the tightly coupled and intermediate regimes becomes more and more important as the density increases. The level of turbulence, as well as the threshold between these two regimes, depends on the gas density through the stopping time (Eq. (9)). For clarity, we recall here the expressions for the mean quadratic relative velocity between the projectile and the target for the two first regimes (Ormel & Cuzzi 2007):

Fig. 2 Grain velocities as a function of the grain radius for a homogeneous grain of specific density 3 g cm −3. The model of grain acceleration by turbulence Ormel & Cuzzi (2007) is presented at different densities. Thick lines represent the relative rms velocity between the grain and the gas, while thin lines represent the relative rms velocity between grains of the same size. The contribution of Brownian motion to the velocity of grains is almost always negligible. The velocity limit for coagulation is presented in two versions: the theory from Chokshi et al. (1993) and the experiments from Poppe & Blum (1997). The sound velocity (T = 10 K) in this isothermal phase is added for reference.

Tightly coupled particles (τ_t < τ_η)

$$\Delta V_{\text{P},\text{T}}^2=\frac{3}{2}{c}_{\text{s}}^2\sqrt{\text{Re}}{\left(\frac{{\tau }_{\text{T}}}{{\tau }_{\text{L}}}-\frac{{\tau }_{\text{P}}}{{\tau }_{\text{L}}}\right)}^2.$$(10)

Intermediate regime (τ_η ≤ τ_t < τ_l)

$$\begin{array}{l}\Delta V_{\text{P},\text{T}}^2=\frac{3}{2}{c}_{\text{s}}^{2}f\left(\frac{{\tau }_{\text{P}}}{{\tau }_{\text{L}}}\right)\frac{{\tau }_{\text{T}}}{{\tau }_{\text{L}}},\\ f(x)=3.2-(1+x)+\frac{2}{1+x}\left(\frac{1}{2.6}+\frac{x^3}{1.6+x}\right),\end{array}$$

with 0 ≤ τ_P∕τ_L ≤ 1 and 1.97 ≤ f ≤ 2.97 (Ormel et al. 2009).

Figure 2 shows how the rms relative velocity between grains depends on the grain radius, for increasing values of the gas proton density. The grain rms velocity increases with the grain radius and decreases with the gas density, and is in most cases strongly sub-sonic. If the projectile and target are of the same size, they can still present a relative velocity in the intermediate regime (thick lines), but not in the tightly coupled regime (thin lines). A target and a projectile hitting each other with a relative velocity larger than velocity limit for coagulation (dashed blue lines) wouldnot stick together but rather rebound, the presence of ice coating increasing this limit velocity (see Sect. 2.5). The relative velocities generated by the thermal motion of grains (supposedly at the gas temperature T = 10 K) decrease strongly with the grain size: $v_{\text{thermal}}\simeq 0.1{(a/0.1\mu \text{m})}^{-3/2}$ cm s⁻¹, for a grain specific density of 3 g cm⁻³. They are small compared to the relative velocities generated by the gas turbulence, except at high densities.

Fig. 3 Grainvelocities (left axis) and Hall factor (right axis) as a function of the grain radius for a gas proton density nH of 104 cm −3 (left), 106 cm−3 (middle) and 108 cm−3 (right). The ambipolar diffusion velocity VAD is fixed and equal to cs∕10 = 20 m s−1. We present the components V∥ parallel to the infall (magenta) and V⊥ perpendicular to the infall and to the magnetic field direction (green). For comparison, the turbulent relative velocities ΔV from Fig. 2 are overplotted in dotted blue lines. The grain specific density is 3 g cm −3.

2.6.2 Grain drift velocity through ambipolar diffusion

Turbulence is not the only way to accelerate or decouple the different grain sizes from each other in dense clouds. Ambipolar diffusion, which appears as soon as the ionisation is too weak to couple the gas to the magnetic field, can also affect the grain dynamics, as it was demonstrated in the case of MHD C-shocks (Guillet et al. 2007).

Grains are charged in the ISM, and as a consequence gyrate around magnetic fields (Spitzer 1941). In dense clouds, their charge is usually negative and close to − e (Flower & Pineau des Forêts 2003), with e the elementary electric charge. Grains are also coupled to the gas through the impact of gas particles. The Hall factor Γ = ω τ, defined as the product of the grain stopping time τ by the grain gyro-frequency ω, characterises the degree of coupling of the grain with the magnetic field and the gas. The Hall factor depends strongly on the grain size (Γ ∝ a⁻², Guillet et al. 2007). Small dust grains (those with Γ ≫ 1), are strongly coupled to the magnetic field and follow on average the dynamics of ions, spiralling around magnetic field lines and participating with ions to the coupling of the gas with the magnetic field through collisions with gas particles. Large grains (those with Γ ≪ 1) follow the gas on average, being almost insensitive to the magnetic field despite their electric charge.

In a model where the magnetic field direction (along b) is transverse to the direction e_r of the collapse, and ignoring the inertia of dust grains, we can express the drift velocity v_k of each grain size k as a function of the ambipolar diffusion velocity V_AD ≡V_n −V_i (Eq. (21) from Guillet et al. 2007) $$v_k=V_{\text{n}}-\frac{{\Gamma }_k^2}{1+{\Gamma }_k^2}{V}_{\text{AD}}+\frac{{\Gamma }_k}{1+{\Gamma }_k^2}{V}_{\text{AD}}\times b,$$(13)

where V_i is the velocity of ions and V_n that of neutrals along the direction e_r of the collapse. We note that the supplementary relative velocities generated by the turbulence between dust grains and gas particles were ignored in this derivation.

The grain drift velocity induced by ambipolar diffusion is not a stochastic velocity, but a systematic velocity. It has two components: one, V_∥, along the collapse direction e_r, and the other, V_⊥, along e_ϕ which is perpendicular to the direction of the collapse and to the magnetic field direction. Figure 3 presents how the component V_∥ and V_⊥ of the grain drift velocity depend on the grain radius, for increasing values of the proton density n_H, assuming a constant ambipolar diffusion velocity V_AD = c_s∕10 = 20 m s⁻¹, as per Basu & Mouschovias (1995). The parallel and perpendicular components of the relative velocity between the grain and the gas do not scale the same way with the Hall factor Γ of the grain (see Eq. (13)): V_∥ is maximal and equal in amplitude to V_AD when Γ ≫ 1, while V_⊥ is maximal and equal in amplitude to V_AD∕2 when Γ = 1. As Γ decreases with increasing density, the systematic velocities generated by ambipolar diffusion tend to disappear with the collapse.

2.7 Calculation of the ambipolar diffusion velocity V_AD

The amplitude V_AD of the ambipolar diffusion velocity depends on the coupling of charged particles with the gas, and therefore on the ionisation degree of the gas and on the size distribution and charge of dust grains. In high-density regions protected from UV radiations, the ionisation can only result from the energetic electrons produced by the interaction of cosmic rays with the gas (Padovani et al. 2018). As the density increases, the ionisation fraction n_i ∕n_H decreases because the recombination rate of ions with electrons or charged grains scales as $n_{\text{H}}^2$, while the ionisation rate scales with ζ n_H, with ζ the cosmic ionisation rate. The low ionisation does not allow a strong coupling between the gas and the magnetic field, and a velocity difference appears at small scales between the ion and neutral fluids, a phenomenon called ambipolar diffusion. Studying the Class 0 protostar B335 with ALMA observations, Yen et al. (2018) inferred an upper limit of 0.3 km s⁻¹ for the ion-neutral drift in this particular object.

The amplitude V_AD of the ion-neutral drift velocity can be calculated by equalising the Lorentz Force acting on charged particles (namely ions and dust grains) and the drag force exerted by the gas particles onto these charged particles. This yields: $$\frac{|B\times \left(\nabla \times B\right)|}{4\pi }={\rho }_{\text{n}}\left(V_{\text{AD}}{n}_{\text{i}}{\widehat{\sigma }}_{\text{in}}+\sqrt{\frac{8k_{\text{B}}T}{\pi {\mu }_{\text{n}}}}{\sum }_k{n}_k{\widehat{\sigma }}_k|V_{\text{n}}-v_k|\right),$$(14)

where ρ_n is the mass density of the gas, T is the gas temperature, n_k and ${\widehat{\sigma }}_k$ are the number grain density and mean grain cross-section in bin k, respectively, $${\widehat{\sigma }}_{\text{in}}=2.41\pi e\sqrt{\frac{{\mu }_{\text{i}}{\mu }_{\text{n}}}{{\mu }_{\text{i}}+{\mu }_{\text{n}}}{\alpha }_{\text{n}}}$$(15)

is the ion-electron collisional cross-section, e the electric charge of ions, μ_i and μ_n the mean mass of ions and neutral particles, respectively, and α_n the polarizability of neutrals (Flower & Pineau des Forêts 2003),

Based on Fig. 3, we neglect the second term in the RHS of Eq. (13) to remain linear² in V_AD and obtain $$|V_{\text{n}}-v_k|\simeq \frac{{\Gamma }_k^2}{1+{\Gamma }_k^2}{V}_{\text{AD}}.$$(16)

As a consequence, V_AD can be approximated by: $$V_{\text{AD}}\simeq \frac{|B\times \left(\nabla \times B\right)|}{4\pi {\rho }_{\text{n}}\left(n_{\text{i}}{\widehat{\sigma }}_{\text{in}}+\sqrt{\frac{8k_{\text{B}}T}{\pi {\mu }_{\text{n}}}}{\sum }_k{n}_k{\widehat{\sigma }}_k\frac{{\Gamma }_k^2}{1+{\Gamma }_k^2}\right)}.$$(17)

The value of V_AD scales as the square b² of the intensity of the magnetic field in the parent cloud (Eq. (4)). It also depends, through the gradient, on the geometry of the magnetic field, as well as on its evolution with the compression of the gas during the infall. For simplicity, we will start the collapse with $V_{\text{AD}}=V_{\text{AD}}^{\text{ref}}=c_{\text{s}}/10=20$m s⁻¹ at n_H = n_H,0 (Basu & Mouschovias 1995), and assume $B\times \left(\nabla \times B\right)\propto B^2/L_{\text{Jeans}}\propto n_{\text{H}}^{3/2}$. This scaling results in two asymptotic regimes: If the ion term dominates over the dust term at the denominator of Eq. (17), then $n_{\text{i}}\propto \sqrt{n_{\text{H}}}$ and V_AD will be approximately constant; on the contrary, if the dust term dominates, then $V_{\text{AD}}\propto 1/\sqrt{n_{\text{H}}}$ as long as the dust size distribution does not change and remains strongly coupled to the magnetic field (Γ_k ≫ 1).

2.8 The coagulation kernel

It is the coagulation kernel K, and therefore in our case the mean relative velocity $\langle \Delta v\rangle$ between theprojectile and the target grains, that control the coagulation rate (see Sect. 2.3 and Eq. (5)). Only those collisions at a relative velocity lower than the velocity limit $\mathcal{V}$ will lead to coagulation, as grains are expected to bounce, and not to stick, at higher velocities (see Sect. 2.5).

2.8.1 Turbulence alone

To express the coagulation kernel in the case where only turbulence accelerates dust grains, we follow the approach by Flower et al. (2005). We assume that the grains relative velocities along x, and y and z are Gaussian variables of variance ${\text{ΔV}}_{\text{P},\text{T}}^2/3$ where ΔV_P,T is the rms relative velocity between the target and the projectile (see Sect. 2.6.1), and only account for collisions with a relative velocity smaller than the velocity limit $\mathcal{V}$. In this frame, the demonstration by Flower et al. (2005) can be extended to the general case of grains of different sizes (and therefore velocities), with the same result: $$\frac{\langle \Delta v\rangle }{{\text{ΔV}}_{\text{P},\text{T}}}=\frac{1}{{\text{ΔV}}_{\text{P},\text{T}}}{\displaystyle {\int }_0^{\mathcal{V}}v}f(v)\text{d}v=\sqrt{\frac{8}{3\pi }}\left(1-\left(1+{\chi }^2\right){\text{e}}^{-{\chi }^2}\right),$$(18)

with $\chi \equiv \sqrt{3/2}\mathcal{V}/{\text{ΔV}}_{\text{P},\text{T}}$.

2.8.2 Turbulence and systematic velocity added in quadrature

If we assume for simplicity that the presence of ambipolar diffusion does not change the nature (hydrodynamic) of the turbulent cascade, Eq. (18) can be generalized to the case of turbulence velocities in the presence of a systematic differential velocity μ added in quadrature to the turbulent velocity of standard deviation ΔV_P,T (see the demonstration in Appendix C) $$\begin{array}{lll}\frac{\langle \Delta v\rangle }{{\text{ΔV}}_{\text{P},\text{T}}}\hfill & =\hfill & \frac{1}{\sqrt{3\pi /2}}\left({\text{e}}^{-{\xi }^2}-\frac{{\text{e}}^{-{\left(\chi -\xi \right)}^2}\left[1+\frac{\chi }{\xi }\right]+{\text{e}}^{-{\left(\chi +\xi \right)}^2}\left[1-\frac{\chi }{\xi }\right]}{2}\right)\hfill \\ \hfill & \hfill & +\frac{1}{\sqrt{3/2}}\left(\xi +\frac{1}{2\xi }\right)h(\chi ,\xi ),\hfill \end{array}$$(19)

with $\xi \equiv \sqrt{3/2}\mu /{\text{ΔV}}_{\text{P},\text{T}}$, $h(\chi ,\xi )\equiv \text{erf}\left(\xi \right)-\frac{\text{erf}\left(\chi +\xi \right)-\text{erf}\left(\chi -\xi \right)}{2}$ and erf the error function.

3 Grain charge, ionisation balance and plasma MHD resistivities

The charging of grains is, in dense clouds, a prerequisite to solve the ionisation of the plasma.

3.1 Charging processes

As per Guillet et al. (2007), the grain charge is computed accounting for: (1) the collection of thermal electrons, with sticking coefficient s_e = 0.5 (a value which is not taken to be 1 owing to the non-planar nature of dust grains, see Draine & Sutin 1987); (2) the collection (and subsequent recombination) of thermal ions, with sticking coefficient s_i = 1; 3) the photoelectric detachment of electrons by Cosmic Rays-induced secondary photons (see Sect. 3.2). If we momentarily ignore the photoelectric detachment of electrons, the mean grain charge is primarily controlled by the ratio of electrons and ions fluxes, by the ratio m_e n_e∕(m_i n_i), where m_e is the electron mass and m_i the mean ion mass.

The DUSTDaP code can integrate the charge distribution f(Z) as well as the mean charge ⟨Z⟩ of a grain (Guillet et al. 2007). The charge distribution of large grains (a > 0.25 μm) is almost Gaussian and can therefore be approximated by its mean charge ⟨Z⟩. Small grains (a < 0.25 μm), however, have a non-Gaussian charge distribution, centred around Z = −1, and intrinsically discrete in nature. Being the dominant charge carriers among grains, their full charge distribution must be integrated (Guillet et al. 2007).

Figure 4 presents the charge distribution of grains before the collapse starts (n_H = 10⁴ cm⁻³). Grains are on average negatively charged, with a mean charge number $\langle Z\rangle \simeq -1$. The charge distribution is not Gaussian, being skewed towards positive values due to the photoemission of electrons by impinging secondary UV photons (see the next section). The larger the grain, the more positive its charge number Z. This is due to the fact that the coulombian enhancement of the grain cross-section $\tilde{\sigma }=1+Z{e}^2/a{k}_{\text{B}}T$ (Spitzer 1941; Draine & Sutin 1987) for the capture of an electron at temperature T by a grain of radius a and charge Z decreases with a when the grain is positively charged, while the photoemission yield is almost independent of a and Z (Weingartner & Draine 2001).

The density of ions is calculated as the sum of the density of positively charged (molecular and atomic) gas species. Knowing the charge of each grain size, the density of free electrons can then be derived from the electroneutrality equation: n_i − n_e + ∑_kn_kZ_k = 0, where n_k and Z_k are the grain density and the mean grain charge in bin k, respectively.

Fig. 4 Grain charge distribution f(Z) for four representative grains sizes a in the cloud envelope (nH = 104 cm−3, T = 10 K, ζ = 5 × 10−17 s−1). The non-Gaussian tails at positive charges is the result of photoemission by CR-induced secondary UV photons (see Sect. 3.2).

3.2 Update of the ionisation rates by cosmic rays

The supra-thermal electrons generated by cosmic rays excite H₂ molecules, which relaxe by producing UV photons, which are called secondary to differentiate them from primary starlight UV photons. Since Flower & Pineau des Forêts (2003), we updated the ionisation rates by cosmic rays using the new tables provided by Heays et al. (2017). The cross-section per unit volume to secondary photons for gas and dust particles are: $$\begin{array}{lll}{\langle n\widehat{\sigma }\rangle }_{\text{gas}}^{\text{secpho}}\hfill & =\hfill & {\sum }_{j=1}^{N_{\text{species}}}\frac{{\gamma }_j}{0.15}{\left(\frac{T}{300}\right)}^{{\alpha }_j}\times 2\times {10}^{-21},\hfill \\ {\langle n\widehat{\sigma }\rangle }_{\text{dust}}^{\text{secpho}}\hfill & =\hfill & {\sum }_{k=1}^{N_{\text{bins}}}{n}_k{\widehat{\sigma }}_k{Q}_{\text{abs}}(a_k,\text{10}\text{eV}).\hfill \end{array}$$

The total cross-section to secondary photons, ${\langle n\widehat{\sigma }\rangle }_{\text{gas}}^{\text{secpho}}+{\langle n\widehat{\sigma }\rangle }_{\text{dust}}^{\text{secpho}}$, depends on the size distribution of dust grains, as well as on the gas content. It is, however, dominated by dust, as usually assumed (Gredel et al. 1989). Since the ionisation rates were provided by these authors under the assumption of an MRN size distribution, each ionisation rate γ_j for the gas specie j (Eq. (20)) must be multiplied by the correcting factor $2\times {10}^{-21}{n}_{\text{H}}/\left({\langle n\widehat{\sigma }\rangle }_{\text{gas}}^{\text{secpho}}+{\langle n\widehat{\sigma }\rangle }_{\text{dust}}^{\text{secpho}}\right)$, to adapt to the evolving size distribution in the collapsing core. Similarly, for the grain charging, Eq. (4) of Guillet et al. (2007) must be adapted by replacing $\langle n\widehat{\sigma }{Q}_{\text{abs}}\rangle$ at the denominator by ${\langle n\widehat{\sigma }\rangle }_{\text{gas}}^{\text{secpho}}+{\langle n\widehat{\sigma }\rangle }_{\text{dust}}^{\text{secpho}}$.

Fig. 5 Characteristical timescales as a function of the proton density, for the larger and smaller dust grain of the initial size distribution (of radius 0.22 and 0.015 μm respectively): gyration timescale (τgyr), stopping timescale τdrag, and charge fluctuation timescale (τZ), compared to the free-fall timescale τff.

3.3 Charge fluctuations

Grain charge fluctuates due to the stochastic nature of the grain charging processes. The timescale of charge fluctuations is larger for small grains than for large grains, owing to the reduced geometric cross-section of the latter. According to Yan et al. (2004), it is $${\tau }_Z=\frac{\text{Var}(Z)}{{\sum }_{Z}f(Z)\times \left(J_e(Z)+J_i(Z)+J_{\text{pe}}(Z)\right)},$$(22)

where Z is the grain charge number, Var(Z) is the variance of the charge distribution, and J_e(Z), J_i (Z), and J_pe (Z) are the rates forelectron collection, ion recombination on the surface of grains, and photoelectric detachment by secondary photons, respectively (see Sect. 3.1).

In Fig. 5, we show that the charge fluctuation timescale is smaller than the dynamical time of the collapse τ_ff, and also smaller than the grain Larmor timescale τ_gyr = 1∕ω for all grain radius considered here and at all gas densities. The charge fluctuation timescale is also smaller than the stopping timescale τ_drag for large grains, but not for small grains (a = 10 nm) at densities higher than 10⁸ cm⁻³. As we shall see, those grains have already been removed by coagulation when this condition is not met anymore. We will therefore ignore charge fluctuations in our study.

3.4 Hall factor for electrons, ions and grains

Let j be a particular species (ion, electron, or grain), τ_j its stopping time, ω_j its Larmor pulsation, and Γ_j = ω_j τ_j its Hall factor. Following Marchand et al. (2016), we have: $$\begin{array}{l}{\tau }_j=\frac{1}{a_{j\text{He}}}\frac{m_j+m_{{\text{H}}_2}}{m_{{\text{H}}_2}}\frac{1}{n_{\text{H}}{\langle \widehat{\sigma }v\rangle }_j},\\ {\omega }_j=\frac{q_{j}B/c}{m_j},\\ {\Gamma }_j={\omega }_j{\tau }_j.\end{array}$$

The factor a_j,He accounts for collisions with helium atoms and is equal to 1.14 for ions, 1.16 for electrons and 1.28 for grains (Desch & Mouschovias 2001).

The Hall factor Γ_j expresses the degree of coupling of the particle to the magnetic field. Its value is very high for electrons and ions (Γ_e ≫ Γ_i ≫ 1). For grains in our core collapse, ${\Gamma }_j\propto B/\left(a^2{n}_{\text{H}}\right)$ where a is the grain size. Grains therefore tend to decouple from the magnetic field as n_H increases during the collapse. Grain coagulation increases this trend by increasing the grain size a, and so does ice accretion by increasing a and decreasing the grain specific density μ_dust.

3.5 Plasma conductivities and resistivities

We note³ σ the conductivity of a charged fluid. The conductivity σ_j of the species j is by definition $${\sigma }_j=n_j\frac{q_j^2}{m_j}{\tau }_j.$$(26)

We note theuseful relation $$\frac{{\sigma }_j}{{\Gamma }_j}=\frac{n_j{q}_j}{B/c}.$$(27)

The parallel, perpendicular and Hall conductivites of the plasma have the following expressions (e.g. Marchand et al. 2016) $$\begin{array}{l}{\sigma }_{\parallel }={\sum }_j{\sigma }_j,\\ {\sigma }_{\perp }={\sum }_j{\sigma }_j\frac{1}{1+{\Gamma }_j^2},\\ {\sigma }_{\text{H}}=-{\sum }_j{\sigma }_j\frac{{\Gamma }_j}{1+{\Gamma }_j^2},\end{array}$$

where the sum is done over all charged species, namely electrons, ions and grains of all sizes.

The Ohmic (parallel) conductivity is dominated by electrons, the most mobile particles. The perpendicular conductivity is dominated by the particles with low values of Γ, that is to say by the particles the least coupled to the magnetic field, which are here dust grains.

We now derive approximate expressions for the Hall conductivity σ_H building on the fact that electrons and ions are strongly coupled to the magnetic field with Γ_e ≫ Γ_i ≫ 1. In the absence of dust grains, using the electroneutrality equation Eq. (30) becomes $${\left({\sigma }_{\text{H}}\right)}^{\text{nodust}}\simeq \frac{n_{\text{i}}e}{B/c}\frac{1}{{\Gamma }_i^2}\propto n_{\text{i}}\frac{n_{\text{H}}^2}{B^3}.$$(31)

In the presence of charged grains, $${\sigma }_{\text{H}}\simeq \frac{{\sum }_k{n}_k{Z}_k/\left(1+{\Gamma }_k^2\right)}{B/c}.$$(32)

The Hall conductivity is therefore controlled by the abundance and sizes of dust grains.

The Ohmic, Hall and ambipolar diffusion resistivities write (Marchand et al. 2016) $${\eta }_{\Omega }=\frac{1}{{\sigma }_{\parallel }},$$(33) $$\begin{array}{l}{\eta }_{\text{H}}=\frac{{\sigma }_{\text{H}}}{{\sigma }_{\text{H}}^2+{\sigma }_{\perp }^2},\\ {\eta }_{\text{AD}}=\frac{{\sigma }_{\perp }}{{\sigma }_{\text{H}}^2+{\sigma }_{\perp }^2}-\frac{1}{{\sigma }_{\parallel }}.\end{array}$$

4 Results

We now present our results for the evolution of the grain size distribution, plasma ionisation, and MHD conductivities and resistivities during the collapse of a spherical core for different scenarios of dust dynamics and evolution. In our standard model, the initial density is n_H,0 = 10⁴ cm⁻³, the initial size distribution is the MRN covered by a 8.8 nm ice mantle, representing a total of 0.9% of the gas mass (see Sect. 2.1). The collapse timescale is the free-fall timescale τ_ff computed at a density n_H,0, and the cosmic ray ionisation rate is ζ = 5 × 10⁻¹⁷ s⁻¹, taken to be constant through the collapse. Before entering into the details of the feedback of dust coagulation on the MHD properties of the collapsing core, we recall what would happen without dust.

4.1 Collapse without dust

Figure 6 presents the variation of the ionisation fraction, MHD conductivities and resistivities with the gas density during the collapse, for a plasma without dust. The ionisation fraction (left panel) results from an equilibrium between the recombination of ions in the gas phase (at a rate $\propto n_{\text{i}}^2$) and the ionisation by CRs (at a rate ∝ ζ n_H), leading to $${\left(\frac{n_{\text{i}}}{n_{\text{H}}}\right)}^{\text{nodust}}\propto \sqrt{\frac{\zeta }{n_{\text{H}}}}.$$(36)

Without grains, gas species cannot deplete, and heavy ions like S⁺ and Fe⁺ remain the main ions through the collapse.

The centre panel of Fig. 6 shows the monotonic evolution of the perpendicular and Hall plasma conductivities. As Γ_e ≫ Γ_ions ≫ 1, the plasma perpendicular conductivity is that of ions, and is therefore proportional to the ion density n_i. Combining Eqs. (4) and (31), the Hall conductivity σ_H is found to be proportional to the gas density n_H. These scaling are well reproduced in the centre panel of Fig. 6. The Hall conductivity is positive because the particles the less coupled to the magnetic field (here ions) are positively charged. We note that the perpendicular conductivity is higher than the Hall conductivity at all densities.

Regarding the plasma resistivities (Fig. 6, right panel), the ambipolar diffusion resistivity is dominating the Hall and Ohmic at all densities below 10¹⁰ cm⁻³. Our model reproduces the expected relations inferred from Eqs. (28)–(35) in the absence of dust grains: η_Ω ∝ n_H∕n_e, ${\eta }_{\text{H}}\propto n_{\text{H}}^2/(B^3{n}_{\text{i}})$ (constant in our collapse model because $B\propto \sqrt{n_{\text{H}}}$ and $n_{\text{i}}\propto \sqrt{n_{\text{H}}}$), and η_AD ∝ 1∕n_i.

4.2 Collapse without coagulation

We now compare the predictions for two kind of models with dust but without coagulation: (1) an MRN size distribution (5− 250 nm), and (2) a truncated-MRN (100−250 nm) which has ~10³ times fewer grains than the MRN. All grains in the parent cloud are covered by icy mantles, 8.8 nm thick for the former, and 50 nm thickfor the latter. The truncated-MRN, which was proposed by Zhao et al. (2016), can be considered as a dust model where the removal of small grains would have already happened in the parent cloud before the collapse.

Figure 7 presents our results for the MRN (bottom row) and truncated-MRN (top row) case. In the presence of dust grains, the ionisation fraction (left column) is not only determined by the CR ionisation rate and the density, but also by the total surface carried by dust grains which serve as a catalyst for ion recombination. This explains why the ionisation fraction is one order of magnitude smaller in the MRN case than in the truncated-MRN case, which is itself much smaller than the no-dust case. The fractional abundance of negative grains remains almost constant with the density, while the ionisation fraction decreases with the density. As a consequence, there appears a shortage of free electrons leading to a situation where dust grains become the main charge carrier and dominate the ionisation, a situation known as a dust-dust plasma where $n_{\text{e}}=n_{\text{i}}/\sqrt{m_{\text{i}}/m_{\text{e}}}$ remains constant (Ivlev et al. 2016). The higher the number of grains, the earlier this shortage appears: from n_H = 10⁶ cm⁻³ for the MRN case and from n_H = 10¹¹ cm⁻³ for the truncated-MRN case. In such a situation, the mean grain charge is not constant anymore and progressively converges towards zero.

The total perpendicular conductivity of the plasma, sum of that of ions and dust grains (Fig. 7, centre column), is a key ingredient to understand the evolution of non-ideal MHD effects during the collapse (see Eqs. (34) and (35)). Before the collapse starts, the total perpendicular conductivity of the plasma is lower in the MRN case than in the truncated-MRN case, which is itself lower than in the no-dust case. This is consequence of the fact that the total conductivity is dominated by ions, the density of which is anti-correlated with the abundance of small grains. To analyse the evolution of the plasma conductivities as the density increases during the collapse, we must distinguish two asymptotic regimes in the coupling of dust grains to the magnetic field determined by their Hall factor Γ ∝B∕(n_H a²):

Fig. 6 Results for a model without dust. Left: fractional abundances of ions ni ∕nH and electrons ne∕nH together with the fractional abundance of the dominant ions Fe+ and S+. Centre: perpendicular (σ⊥) and Hall (σH) conductivities. Right: Ohmic (ηO), Hall (ηH) and ambipolar diffusion (ηAD) resistivities. Asymptotic trends are over-plotted in dotted lines (see text for a derivation of those expressions). Horizontal lines are added to help the reader estimate the plasma resistivities for a given proton density.

Fig. 7 Results for our models without dust evolution, with truncated-MRN (100−250 nm) (top) and MRN (5−250 nm) (bottom) size distributions, both covered by icy mantles. Left: evolution of the fractional abundances of charged species: ions, electrons, negative charges (g-) and positive charges (g+) on grains. Centre: evolution of the dust parallel conductivity (σg∥), dust perpendicular conductivity (σg⊥), ion perpendicular conductivity (σi⊥) and Hall conductivity (σH). Scalings in powers of nH are over-plotted in dotted lines for an easier comparison of our numerical results with the asymptotic trends described in Sect. 4.1). Right: evolution of the Ohmic (ηΩ), Hall (ηH), and ambipolar diffusion (ηAD) resistivities.

Strong coupling (Γ ≫ 1, small grains or low density)

$$\begin{array}{lll}{\left({\sigma }_{g\perp }\right)}^{\Gamma \gg 1}\hfill & \propto \hfill & \frac{n_{\text{H}}}{B^2}{\sum }_k{n}_k{\widehat{\sigma }}_k{Z}_k\hfill \\ {\left({\sigma }_{\text{H}}\right)}^{\Gamma \gg 1}\hfill & \propto \hfill & \frac{1}{B}{\sum }_k\frac{n_k{Z}_k}{{\Gamma }_k^2}\propto \frac{n_{\text{H}}^2}{B^3}{\sum }_k\frac{n_k{\widehat{\sigma }}_k{}^2}{Z_k}\hfill \end{array}$$

Weak coupling (Γ ≪ 1, large grains or high density)

$$\begin{array}{lll}{\left({\sigma }_{g\perp }\right)}^{\Gamma \ll 1}\hfill & \simeq \hfill & {\sigma }_{g\parallel }\propto \frac{1}{n_{\text{H}}}{\sum }_k\frac{n_k{Z}_k^2}{{\widehat{\sigma }}_k}\hfill \\ {\left({\sigma }_{\text{H}}\right)}^{\Gamma \ll 1}\hfill & \propto \hfill & \frac{1}{B}{\sum }_k{n}_k{Z}_k\propto \frac{n_{\text{e}}-n_{\text{i}}}{B}.\hfill \end{array}$$

The transition from the strong to the weak coupling regime happens when the mean, density-weighted, grain Hall factor is close to 1. When dust grains are strongly coupled to the magnetic field (Γ_k ≫ 1), the total perpendicular conductivityof dust grains σ_g⊥ approximately scales as n_H while the ion perpendicular still scales as n_i, more slowlywith the density. When dust grains are weakly coupled to the magnetic field (Γ_k ≪ 1), the dust perpendicular conductivityσ_g⊥ is equal to the dust parallel conductivity σ_g∥, and therefore only depends on the grain mean charge if the size distribution is kept fixed.

These two different regimes can explain the differences in the evolution of the plasma conductivities observed between the MRN (where initially Γ ≫ 1) and truncated-MRN (where Γ ≪ 1) cases. Overall, the conductivities of the plasma increase faster with the density in the MRN case than in the truncated-MRN case, becoming higher for n_H ≥ 10⁶ cm⁻³ despite a lower initial value in the parent cloud. The large (a ≥ 0.1 μm) grains of the truncated-MRN case are already weakly coupled to the magnetic field at the beginning of the collapse. Furthermore, their charge remains almost constant during the collapse because the ratio n_e ∕n_i is itself constant (see left panel and Sect. 3.1). As a consequence, the dust perpendicular conductivity remains constant in the truncated-MRN case and is totally negligible at all densities compared to the ions perpendicular conductivity. The conductivities and resistivities are therefore very similar to the no-dust case regarding their evolution trends with the gas density, only differing in their amplitude which is a function of the ionisation degree of the plasma. On the contrary in the MRN case that entails a high abundance of small, strongly coupled, grains, the perpendicular conductivity initially dominated by dust grains first increases proportionally to the density. Then, as soon as the depletion of electrons appears, the mean grain charge converges towards zero which severely decreases the perpendicular conductivity of dust grains (Fig. 7 for densities between 10⁸ and 10¹⁰ cm⁻³). The dust perpendicular conductivity eventually becomes lower than the ions perpendicular conductivity, reaching a plateau when the charge distribution is stabilized (Fig. 7, left panel, n_H ≥ 10¹¹ cm⁻³).

The Hall conductivity, like the perpendicular conductivity, is affected by the degree of coupling of dust grains with the magnetic field. When dust grains are strongly coupled to the magnetic field, the Hall conductivity strongly increases with the density, scaling approximately as $n_{\text{H}}^{3/2}$ (Eq. (38) ignoring variations in the grain charge Z_k). This is what is observed at low density in the MRN case. When grains decouple from the magnetic field at higher density, the Hall conductivity does not increase so fast, or even decrease. In the truncated-MRN case we get ${\sigma }_{\text{H}}\propto n_{\text{H}}^{1/2}$ (Eq. (40) with n_k Z_k ∝ n_H), while in the MRN case the dusty plasma regime yields ${\sigma }_{\text{H}}\propto n_{\text{H}}^{-1/2}$ (Eq. (40) with n_e − n_i remaining constant, see Fig. 7 centre bottom panel).

The right panels of Fig. 6 shows the evolution of the plasma MHD resistivities for the MRN and truncated-MRN cases. Before the collapse starts, the Hall and ambipolar diffusion resistivities are higher in the MRN case than in the truncated-MRN case (and are higher in the truncated-MRN case than in the no-dust case) because, as explained above, the presence of small grains tends to decrease the ion density and therefore the ion conductivity. As the collapse proceeds the situation, however, reverses. Regarding the scaling of the resistivities with the density, the truncated-MRN case is similar to the no-dust case (Fig. 6 and Sect. 4.1) because the perpendicular conductivity is dominated by ions all through the collapse (σ_⊥∝ n_i). In the MRN-case at densities lower than 10⁹ cm⁻³ (σ_⊥∝ n_H), the perpendicular conductivity is dominated by dust grains. As a consequence, the Hall and ambipolar diffusion resistivities decrease more steeply with the density: η_AD ≃ 1∕σ_⊥∝ 1∕n_H and ${\eta }_{\text{H}}\simeq {\sigma }_{\text{H}}/{\sigma }_{\perp }^2\propto 1/\sqrt{n_{\text{H}}}$. These distinct scalings result in values of the MHD resistivities at intermediate densities (10⁸ cm⁻³) that can be a factor of 10 to 100 lower in the MRN case than in the truncated-MRN case despite a higher value in the parent cloud, with the value for the Hall resistivity eventually exceeding that of ambipolar diffusion for densities higher than ~ 10⁸ cm⁻³ in the MRN case.

Overall, these results are in good agreement with those of Marchand et al. (2016), an other model based on an MRN size distribution without dust evolution (see their Figs. 3 and 5), and also those of Zhao et al. (2016) who found a reduced magnetic braking for a truncated-MRN distribution compared to an MRN distribution.

4.3 Collapse with accretion and coagulation

We now turn to the effect of grain evolution, namely ice accretion and grain–grain coagulation, onto the grain size distribution and onto the MHD properties of the plasma. Figure 8 presents our results for three different scenarios of the dust velocity field: (1) turbulence only, (2) ambipolar diffusion only, and (3) turbulence and ambipolar diffusion (see Sect. 2.6 for a description of grain dynamics by turbulence and ambipolar diffusion).

The first row of Fig. 8 presents the evolution of the size distribution of grain cores (i.e. of the refractory part of the grain, covered by icy mantles) during core collapse when dust coagulation and ice accretion are activated. In the early stage of the infall (Fig. 8, n_H < 10⁶ cm⁻³), coagulation removes small grains by sticking them onto large grains. This does not, however, significantly change the upper limit of the grain size distribution. The removal of small grains is faster when the effect of ambipolar diffusion on grain velocities is included (centre and right) because in that case, unlike with turbulence, small grains decouple from the gas, and therefore also from the large grains which follow the gas. Depending on the density, this decoupling is more intense along B or perpendicular to B (see Fig. 3). We note that, to compare our results with the truncated-MRN, the disappearing of grains smaller than 0.1 μm only happens late in the collapse (n_H ≥ 10¹⁰ cm⁻³). When the acceleration of dust grains by turbulent eddies is included (left and right), grain growth proceeds rapidly, the grain size distribution becoming more and more dominated by a single size approaching ~ 10 μm at n_H = 10¹² cm⁻³. To summarise, turbulence makes grain grow in size at later times but only slowly removes small grains, while ambipolar diffusion efficiently removes smaller grains at early times but does not allow the grain size to increase significantly. In Sect. 5.4, we compare our results for the grain growth rate with the others studies.

The second row of Fig. 8 presents how the total geometrical cross-section of grains evolves during the collapse. The total dust cross-section is a key ingredient of the plasma physics because it controls the recombination rate of ions and the perpendicular conductivity of dust grains (see Eq. (37)), which usually dominates the total perpendicular conductivity. Starting with values characteristic for the MRN size distribution already covered by a 8.8 nm ice mantle (see Sect. 2.1), the total dust cross-section first slightly increases through the growth of icy mantles on the surface of grains (n_H≤ 2 × 10⁴ cm⁻³), a process thatis fast enough to complete before any significant coagulation through grain–grain collisions has occurred. Once started, grain–grain coagulation systematically decreases the total dust cross-section. When drift velocities by ambipolar diffusion are included, the faster removal of small grains accelerates this process at low densities (n_H ≤ 10⁶ cm⁻³). Still, it is only when the acceleration of grains by turbulence is included that the total dust cross-section keeps on decreasing all through the collapse (n_H ≥ 10⁶ cm⁻³). The dust size distribution reaches total cross-sections per H that are comparable to that for our truncated-MRN only late in the collapse (n_H ~ 10⁹−10¹⁰ cm⁻³).

The third row of Fig. 8 shows how the total number of negative charges carried by grains is affected by grain–grain coagulation. Unlike in the case without coagulation (Fig. 7), grains never become the dominant charge carrier characteristic of a dusty plasma (see Sect. 4.2). The ionisation is always dominated by ions, though we can observe at high densities (from n_H ≥ 10⁷ cm⁻³) a small depletion of free electrons, characteristic for a dust-ion plasma (Ivlev et al. 2016), when ambipolar diffusion drift velocities are ignored.

The fourth row presents the evolution of the parallel, perpendicular and Hall conductivities of the plasma. The total parallel conductivity (not shown) is high and always controlled by the electron density. As explained in the previous section, the perpendicular conductivity is dominated by dust grains at low densities, and ions at high densities. The density threshold between these two regimes depends on the abundance of the smaller grains, which are the last grains to decouple from the magnetic field as the density increases. For an initial MRN distribution, this transition happens at low density when the removal of small grains is efficient, as it is the case when grain drift velocities are included (n_H ≃ 10⁶−10⁷ cm⁻³), and at high densities (n_H ≃ 10⁹−10¹⁰ cm⁻³) when only turbulence is considered or, similarly, when coagulation is ignored (Fig. 7). From this analysis, we conclude that an efficient removal of small grains stops the increase in the Hall and ambipolar diffusion conductivities with the density, while the growth of large grains does not in itself affect these conductivities. When the density is high enough, or the mean grain size high enough, the Hall conductivity can change from being negative to being positive, implying that the plasma Hall conductivity is then dominated by ions.

The bottom row presents the evolution of the Ohmic, Hall and ambipolar diffusion resistivities with the density. The removal of small grains limits the drop of the Hall and ambipolar diffusion resistivities with the density, increasing their values by a factor of 10–100 compared with models where the coagulation is only triggered by turbulent velocities (left panel), or with models with an MRN distribution and no coagulation at all (Fig. 7, bottom row). Comparing the amplitude of MHD resistivities in the MRN-case with coagulation and ambipolar diffusion drift velocities (Fig. 8) and in the truncated-MRN case without coagulation (Fig. 7), the ambipolar diffusion resistivities have similar values at intermediate densities (n_H~ 10⁸ cm⁻³) while the Hall resistivity remains at least one order of magnitude higher in the former than in the latter case at densities lower than ~ 10⁹ cm⁻³. At a density of 10⁸ cm⁻³, the MHD properties of the model with coagulation and of the truncated-MRN model still differ significantly, being dominated by dust for the former, and by ions for the latter. To obtain high values for the Hall and ambipolar diffusion resistivities necessary to the formation ofa larger disk, it may therefore not be necessary to assume, as per the truncated-MRN model, that all grains smaller than 0.1 μm have been removed by coagulation in the parent cloud before the collapse: removing only parts of those grains through grain–grain coagulation should also achieve it if this happens early in the collapse, as is the case when grain drift velocities generated by ambipolar diffusion are included.

Figure 9 summarises our results concerning the influence of the efficiency in the removal of small grainson the profiles of the plasma resistivities. We present our results for 5 models starting with the same MRN distribution covered by ices, differing only by the initial (i.e. in the parent cloud) value $V_{\text{AD}}^{\text{ref}}$ of the ambipolar diffusion velocity, from 0 to 20 m s⁻¹. The value $V_{\text{AD}}^{\text{ref}}$ of the ambipolar diffusion velocity is used here as a way to control the efficiency in the removal of small grains. It can also be considered as a way to quantify the balance in the competition between coagulation and fragmentation of grains in higher velocity impacts, that we remind is not included in this work. From Fig. 9, we see that the more efficient the removal of small grains, the higher the Hall and ambipolar diffusion resistivities at intermediate densities (10⁶ < n_H < 10⁹ cm⁻³), here again confirming the results obtained by Zhao et al. (2016). This trend, opposite to what was observed in the parent cloud (see Sect. 4.2), is observed whenever the plasma perpendicular conductivity is dominated by dust grains, and not ions. Our conclusions regarding the evolution of the magnitude of the MHD resistivities during the collapse with dust coagulation therefore drastically depend on the detailed dynamics of small grains, and on the competition between coagulation and fragmentation (see Sect. 5.1 for a discussion of the modelling of grain–grain fragmentation in core collapse).

Fig. 8 Results for our standard model with grain accelerated by turbulence only (left) ambipolar diffusion only (centre), and both mechanisms (right). From top to bottom: dust size distribution at increasing densities, total dust cross-section per H (with our results for the truncated-MRN model over-plotted), ionisation, conductivities and resistivities. See Fig. 7 for a description of the last three rows.

Fig. 9 Impact of the initial value of the ambipolar diffusion velocity on the evolution of the MHD resistivities with the density. Both turbulent and drift velocities are included.

Fig. 10 Impact of the CR ionisation rate ζ on the evolution of the MHD resistivities with the density for a model with relative velocities generated by turbulence (top), and turbulence and ambipolar diffusion with $V_{\text{AD}}^{\text{ref}}$ = 20 m s−1 (bottom).

4.4 Influence of the cosmic ray ionisation rate

The ionisation level, and therefore the MHD conductivities and resistivities, are not only controlled by the grain size distribution, but also by the cosmic ray ionisation rate ζ. In Fig. 10, we present the evolution of the ionisation level and MHD resistivities for four values⁴ of the cosmic ray ionisation rate ranging from 5 × 10⁻¹⁹ to 5 × 10⁻¹⁶ s⁻¹, and two coagulation scenarios (turbulent velocities versus turbulent and ambipolar diffusion velocities) corresponding to an inefficient versus efficient removal of small grains, respectively. We recall that, for simplicity, the parameter ζ is kept constant through the collapse, a limit in our model that we discuss in Sect. 5.2. Figure 10 shows that, as expected, low values of the cosmic ray ionisation rate tend to produce high Hall and ambipolar diffusion MHD resistivities when the dust perpendicular conductivity is negligible compared to that of ions, but have almost no impact when the dust perpendicular conductivity dominates over that of ions.

By modifying the ionisation level, the variation of ζ in turn affects the coupling of the magnetic field with the gas, and therefore the ambipolar diffusion velocity. Figure 11 presents the evolution of V_AD with the density for different values of the cosmic ray ionisation rate and, for clarity, for the same initial value $V_{\text{AD}}^{\text{ref}}=20$m s⁻¹ of the ambipolar diffusion velocity in the parent cloud, meaning that we are interested by the relative, not absolute, evolution of V_AD for different values of ζ. We now build on our analysis of the dependence of V_AD presented in Sect. 2.6.2. For ζ = 5 × 10⁻¹⁶ s⁻¹, the ionisation fraction is so high that the dust term at the denominator of Eq. (17) is negligible compared to the ion term and V_AD is therefore almost constant. For lower values of ζ, the dust term tends to dominate over the ion term, leading to $V_{\text{AD}}\propto 1/\sqrt{n_{\text{H}}}$ at low density. Once the total dust cross-sections has significantly dropped through grain–grain coagulation and dust grains have started to decouple from the magnetic field lines (Γ_k ≃ 1, Eq. (17)), V_AD starts to increase. At high density, the dust term has become negligible and $V_{\text{AD}}\propto 1/\sqrt{\zeta }$. The scenario without dust evolution follows the same trends at low density. At high density, the ion density remains constant and grains become neutral on average (Fig. 7), leading to the asymptotic behavior $V_{\text{AD}}\propto \sqrt{n_{\text{H}}}$.

Fig. 11 Ambipolar diffusion velocity VAD as a function of the density, for different values of the CR ionisation rate ζ and, for clarity, the same initial value in the parent cloud $V_{\text{AD}}^{\text{ref}}=20$ m s−1. Calculations were done for our scenario with turbulent and ambipolar diffusion drift velocities. The scenario without dust evolution is added for reference.

5 Discussion

In this section, we discuss some limitations in our modelling that may affect our conclusions, and some perspectives of our work.

5.1 Impact of fragmentation in grain–grain collisions

For reason of simplicity, our modelling of the evolution of the dust size distribution ignores the possibility that the colliding grains would fragment each other, partially or totally, in the collision. This is a strong hypothesis that is probably not justified. While the fragmentation of solid particles necessitates relative velocities on the order of a few km s⁻¹ (Tielens et al. 1994) which most probably do not exist in dense clouds (Ormel & Cuzzi 2007; Yan et al. 2004), the fragmentation of porous aggregates can happen at much lower velocities (a few m s⁻¹) which are relevant for study of grain growth in dense clouds (Weidenschilling & Ruzmaikina 1993; Ormel et al. 2009). We note that, unlike for solids, the outcome of the collision of two aggregates is not a function of the relative velocity but of the energy involved in the collision and on the sizes of the aggregates (Dominik & Tielens 1997).

If the fragmentation process is efficient under the conditions that govern the dynamics of dust grains in collapsing clouds, it will strongly enhance the abundance of the smaller grains of the size distribution, and therefore drastically affect the ionisation balance and the MHD resistivities of the plasma (see Sect. 4). The effect of fragmentation of small aggregates made of very small monomers is expected to be so important for the resistivities of the collapsing cloud that it deserves a study in itself. Asymptotically, we expect the competition of fragmentation and coagulation to lead to evolutionary paths somewhere in between the “no evolution” and coagulation scenarios presented in this article.

There is currently no model available describing the outcome of a collision between two aggregates composed of small (between 5 and 250 nm) monomers, in particular for the proportion of the grain mass that would be converted into very small fragments. The DUSTDaP code can handle the fragmentation of grains, but the physics encoded is not relevant to the fragmentation of aggregates but to the fragmentation of solids (threshold velocity for fragmentation on the order of a km s⁻¹, Tielens et al. 1994), as it must happen in interstellar shocks (Guillet et al. 2009, 2011). This possibility is not used in this article. A first approach would be to keep this framework but use a low value for the threshold velocity for fragmentation (on the order of a few m s⁻¹), or use a threshold in the impact energy as recommended by Dominik & Tielens (1997).

5.2 Decrease in the cosmic ray ionisation rate ζ with the density

For reason of simplicity again, we kept the cosmic ray ionisation rate, ζ, constant throughout the collapse, though it is known from observations that is tends to decrease with increasing column density (Padovani et al. 2013). The overall effect of an hypothetic decreasing value of ζ with the local densityis to decrease the ionisation fraction and to increase the MHD resistivities (Fig. 10), though this is only an approximation as ζ depends on the total column density encountered by cosmic rays during their propagation, and not on the local density. The increase in the resistivities that we infer from the follow-up of dust coagulation in cloud collapse should therefore be reinforced by a decrease in ζ during the collapse. We note, however, that the DUSTDaP code computes all variables out-of-equilibrium, making it difficult to estimate from Fig. 10 the out-of-equilibrium response of the plasma to a drop of ζ over a timescale much smaller than the free-fall timescale. This aspect will be investigated in a future work.

5.3 Grain dynamics

The impact of grain–grain coagulation on the dust size distribution is very much dependent on the model assumed for the dynamics of dust grains. We have shown for example that the grain velocities induced by hydrodynamical turbulence and those generated by ambipolar diffusion result in very different size distribution, grain growth rate, and feedback on the plasma properties (Fig. 8). Using numerical simulations Pan & Padoan (2015) have reevaluated the models of grain turbulent dynamics for protostellar clouds and protoplanetary disks elaborated since the pioneering work of Voelk et al. (1980). They demonstrate that these models generally over-predict the rms velocity of dust grains by a factor of two. They also show that such dynamics ignore the enhancement of the coagulation rate induced by the increase in the dust-to-gas ratio (an effect caused by the clustering of dust grains, confirmed by other numerical studies, see Hopkins & Lee 2016; Lee et al. 2017; Lebreuilly et al. 2019). Luckily, these two sources of errors cancel out for certain grain sizes, revealing a overall good agreement between the Pan & Padoan (2015) simulation and the Ormel & Cuzzi (2007) model.

Another source of uncertainty in the grain dynamics is the nature of the turbulence assumed for the model (or for the simulation): hydrodynamical (HD) or magnetohydrodynamical (MHD). The model of grain dynamics used here (Ormel & Cuzzi 2007) is based on HD turbulence. Yan et al. (2004) developed a model of grain dynamics under MHD turbulence. In this model, grains are still accelerated by gas vortices, but these vortices are not isotropic anymore because of the presence of the magnetic field, affecting the dependence of the rms grain velocity with the grain radius and gas density. Magnetic acceleration processes such as gyro-resonance (a coupling between magnetic waves and grain gyration, Yan et al. 2004) also bring their contribution to the total rms velocity of dust grains. The main difference between these MHD (Yan et al. 2004) and HD (Ormel & Cuzzi 2007) models of grain dynamics is the magnitude of the rms velocity, not so much its dependence on the grain radius and gas density. Indeed, the injection velocity of the MHD turbulence is assumed to be the Alfven velocity in Yan et al. (2004), which is almost one order of magnitude higher than the sound velocity used in the case of HD turbulence (see Sect. 2.6.1). This factor difference has a straightforward, proportional, effect on the coagulation rate. The physics developed in Yan et al. (2004) is mainly focused on the warm and diffuse phase of the ISM, with an extension to the grain dynamics in dense clouds at the density of n_H = 10⁴ cm⁻³. This density is our initial density before the collapse proceeds. For this reason, and also to facilitate comparisons of our results with other works of dust coagulation in cores, we chose to study the impact of the dynamics produced HD turbulence, not MHD turbulence. Still, our introduction of ambipolar diffusion is already a first step towards the entire inclusion of MHD aspects of grain dynamics, which we leave for a future work.

5.4 Grain growth: comparison with other works

In this section, we compare our results regarding the grain growth rate with other studies of dust coagulation in static and collapsing cores. Figure 8 has shown that hydrodynamical turbulence can, according to our calculations, make grain grow up to ~1−10 μm in a free-fall time. Flower et al. (2005) obtained similar results using another modified version of the Paris-Durham shock code and an approach based on a equivalent, single size, grain that grows in time. Other authors concentrating on grain growth in grain–grain processes towards the building-up of very large (mm) grains, found higher growth rates in their calculations (e.g. Weidenschilling & Ruzmaikina 1993; Ossenkopf 1993; Ormel et al. 2009). Coming after numerous detailed studies on the subject, we did not attempt in this paper to propose a realistic study of grain growth in cores. Our intent was to study the coupling between dust evolution and the ionisation and resistivities of the plasma, concentrating on the evolution of the smallest grains (a ≪ 0.1 μm) of the size distribution that a study of grain growth can reasonably ignore. Still, we can make the following comments and comparison with these works. Ossenkopf (1993) included in its modelling the follow-up of the grain porosity, also adding the effect of gravitation on the grain dynamics and coulombian effects on the grain–grain collision rate, but worked at constant density, and ignored fragmentation. Weidenschilling & Ruzmaikina (1993) considered fractal grains, included other dynamical processes such as the dust settling along the disk and a modelling of collisional fragmentation of aggregates, and considered both static and collapsing cloud, but like Ormel et al. (2009) ignored the evolution of grains smaller than 0.1 μm. Ormel et al. (2009), with a more physical description of the growth of aggregates through fragmentation, coagulation and compaction, started with a single size distribution of bare or ice-coated monomers of 0.1 μm. At constant density, they obtain little grain growth on a free-fall timescale. Over a longer timescale, higher growth rates are obtained as soon as grain enter into a regime of self-enhanced coagulation generated by the increasing grain porosity. Overall, these studies, which include a detailed description of the physics of evolution of dust aggregates, generally find higher growth rates than in our study, as expected (Ormel et al. 2009). Whether this regime of high porosity exists or not in the presence of ambipolar diffusion with a size distribution of monomers like here is an open question. In any case, this higher growth rate, which is due to the effect of aggregate porosity that we ignored here, is most probably not as essential to the calculation of the MHD resistivities of the plasma as the effect of fragmentation that we also ignored because fragmentation more than grain growth controls the number density of the smaller grains (see Sect. 5.1 for a discussion on the importance of dust fragmentation in the evaluation of the MHD resistivities of the collapsing cloud).

5.5 Can the magnetic flux be regulated by the abundance of small grains ?

The feedback of dust coagulation onto the evolution of the ambipolar diffusion velocity through the collapse (Eq. (17) and Fig. 11) raises the possibility of the self-regulation of the magnetic flux in the parent cloud. The ambipolar diffusion velocity scales with the square of the magnetic field intensity. A higher value for the magnetic field intensity in the parent cloud therefore tends to increase the value of V_AD, which increases the coagulation rate between small grains (coupled to the magnetic field) and large grains (coupled to the gas). The removal of small grains in turn decreases the coupling of the magnetic field with the gas, leading to a loss of magnetic flux. The ambipolar diffusion velocity then decreases, which stops the removal of small grains. In particular, the process described here could limit the magnetic intensity inside dense cores. This possibility, which is described here only qualitatively, could be investigated using MHD simulations by following the MHD dynamics and coagulation of a simplified grain size distribution composed of two fluids of grains (of 10 and 100 nm, for example) and building on the recipes for the grain charges and ionisation equilibrium in a dusty plasma from Ivlev et al. (2016).

6 Summary

In this article, extending the work by Marchand et al. (2016), we studied how the evolution of the dust size distribution may affect the MHD properties of a collapsing cloud. For this purpose, we have used the DUSTDaP code, a numerical tool that was designed to follow the evolution of the dust size distribution and its feedback on the evolution of the ionisation, chemical content and dynamics of interstellar shocks (Guillet et al. 2011). The code was amended to follow grain–grain coagulation in an element of volume undergoing isothermal compression in a free-fall collapse. For the simplicity of the analysis, we have for the moment ignored the possible production of small grains by fragmentation and craterisation in grain–grain collisions (Ormel et al. 2009). The variation of the cosmic ray ionisation rate with the column density (Padovani et al. 2013), an important factor influencing the ionisation and resistivities of the plasma, was also ignored.

We have considered two physical models for the grain size-dependent relative velocities responsible for grain–grain collisions and coagulation in the collapsing cloud: stochastic velocities triggered by hydrodynamical turbulence (following the recipe by Ormel & Cuzzi 2007), and systematic velocities induced by ambipolar diffusion (following our approach of grain dynamics in C-type shocks, Guillet et al. 2007). In the former model small grains are coupled to the gas, while they are coupled to the magnetic field and therefore decoupled from the gas and from larger grains in the latter.

Here is what we find: In a first phase of the isothermal collapse (n_H = 10⁴−10⁸ cm⁻³), grain–grain coagulation primarily removes small grains from the size distribution by sticking them onto larger grains. This depletion process, which happens faster and therefore sooner when ambipolar diffusion-induced velocities are self-consistently taken into account, does not significantly modifies the upper limit of the grain size distribution. In a second phase (n_H = 10⁸−10¹² cm⁻³), the mean grain size increases from ~0.1m to ~10 μm owing to thelarge, turbulent-induced, relative velocities between grains of comparable sizes while ambipolar diffusion has no more impact on the grain dynamics.

We have further studied the impact of grain–grain coagulation onto the evolution of the ionisation, conductivities and resistivities of the plasma, and therefore on the coupling of the magnetic field with the collapsing gas. The presence of small dust grains modifies the conductivities and resistivities of the plasma both directly and indirectly; directly through their own contributionto the conductivity, and indirectly through their control over the ionisation fraction and therefore over the conductivity of ions. If the total perpendicular conductivity is dominated by ions (as is typically the case in the parent cloud), removing small grains will decrease the Hall and ambipolar diffusion resistivities. If the total perpendicular conductivity is dominated by dust grains (as is the case during the collapse at densities lower than ~ 10⁷ cm⁻³ when starting with an MRN size distribution), removing small grains will on the contrary increase the MHD resistivities.

Starting with an MRN size distribution in the parent cloud, the amplitude increase in the Hall and ambipolar diffusion resistivities induced by grain–grain coagulation depends on the model of dust dynamics. It is stronger when the drift velocities induced by ambipolar diffusion are taken into account than when only turbulent velocities are considered. The ambipolar diffusion (respectively Hall) resistivity is comparable (respectively higher by one order of magnitude) in our model where dust grains coagulate during the cloud collapse than in a toy model without coagulation during the collapse where the removal of grains smaller than 0.1 μm would have already happened in the parent cloud (truncated-MRN size distribution, Zhao et al. 2016). It is therefore likely that the disappearance of small grains due to the differential velocity induced by ambipolar diffusion will have a significant impact on the formation of planet-forming disks. In particular it likely leads to the formation of somewhat larger disks, even when starting with an MRN size distribution in the parent cloud. Varying the uniform value assumed for the CR ionisation rate ζ in the 5 × 10⁻¹⁹−5 × 10⁻¹⁶ s⁻¹ range underthe simplifying assumption that the ambipolar diffusion velocity is constant in the parent cloud, we find that the Ohm resistivity (which is controlled by free electrons) is, as expected, very sensitive to ζ, but that the Hall and ambipolar diffusion resistivity are almost independent of it as long as the plasma total perpendicular conductivityis dominated by dust grains.

Our study confirms that the MHD resistivities of the collapsing core are very sensitive to the abundance of small grains, and therefore to the mechanisms and grain dynamics that drive their production and removal. The possible fragmentation of aggregates in mutual collision being ignored in our analysis, the modelling of the competition between the coagulation and fragmentation of dust grains appears to be the next important step towards a consistent modelling of the feedback of dust evolution onthe MHD dynamics of a collapsing core.

Acknowledgements

We thank the referee for her/his comments that helped improved the clarity of the manuscript. V.G thanks A. Ivlev and P. Lesaffre for stimulating discussions. This work was supported by the Programme National PCMI of CNRS/INSU.

Appendix A Coagulation rate in Larson’s collapse

In this section, we demonstrate that our model of a uniform compression of the cloud (Sect 2.2) is a good model to follow the dust size distribution induced by grain–grain coagulation through the cloud collapse up to a given density, even if it is a bad model to describe the density structure of the cloud at a given time.

To test our simple uniform compression model against a better representation of a cloud collapse, we compare the results we obtain with our model with those obtained for the isothermal collapse of a core sensitive to gravitational force and gas pressure (Larson 1969). We integrate Larson’s equation in a Lagrangian approach. The cloud density profile is presented in the top panel of Fig. A.1 at four epochs of the collapse. In the Larson scenario, we observe a uniform distribution of the gas density near the cloud centre, surrounded by an envelope where the density scales as R⁻² that extends to a distance close to the initial radius R_C of the cloud.On the contrary, in our simple free-fall model that ignores the effect of the pressure of the gas the density is uniform and the cloud dimension rapidly contracts.

Fig. A.1 Top: density profile of the collapsing cloud at four stages of its evolution for the Larson (Larson 1969) scenario (thick lines), and for our simple collapse model of Sect. 2.2. Bottom: ∫ nH d t for the Larson collapse for gas cells distributed at different radius before the cloud collapse. Results for models of uniform compression assuming nH ∝ R−3 (black thick line) and nH ∝ R−2 (black thin line) are over-plotted.

Our purpose in this article is to estimate the dust size distribution at any density during the collapse, not to describe the density distribution of the cloud. To check the reliability of our model for this purpose, we compare the total number of collisions experienced by a given grain particle in the Larson scenario with that obtained for our simple model.

The number of collisions experienced up to a time τ in the collapse by a given dust grain with other dust grains is primarily controlled by the evolution of the local gas density, which increases by orders of magnitude during the collapse. The scaling of the collisional rate with the relative velocity between impinging particles is of secondary importance for our purpose because these velocities depend on the model used for the grain dynamics and only mildly evolves through the collapse ($\Delta V\propto n_{\text{H}}^{-1/4}$ in the model by Ormel et al. 2009, see Fig. 2). To obtain a first order result, the total number of collisions experienced by a given particle during the collapse up to a time τ can be estimated through the proxy: $$N(\tau )\equiv {\displaystyle {\int }_0^{\tau }{n}_{\text{H}}}(t)\text{d}t.$$(A.1)

The bottom panel of Fig. A.1 compares the evolution of N(τ) in the Larson scenario and for our model of a uniform collapse. For the Larson scenario, the evolution of N(τ) with the local density is presented for six initial positions of the grain at 0.9, 0.8, 0.6, 0.4, 0.2 and 0.1 times the initial radius R_C of the cloud. We also present the evolution of N(τ) for our model of a uniform collapse, which by hypothesis does not depend on the initial position of the grain. The uniform compression scenario nicely reproduces the results obtained in the Larson scenario for all initial positions of the grain. This is clear for particles close to the centre (R∕R_C < 0.6), but also approximately for particles initially close to the boundary, within a factor of 1.2. These results are surprising but easily explained. For the Larson scenario, a gas cell almost stops its compression once reached by the rarefaction wave (top panel of Fig. A.1); time goes on, but the gas density and the dust size distribution do not evolve much anymore. The dust size distribution in a Larson collapse is therefore close to the one obtained for a uniform collapse, as long as one compares the results at the same gas density, and not at the same time.

Appendix B The coagulation algorithm

The evolution of the grain size distribution by grain–grain coagulation is driven by the Smoluchowski equation (Smoluchowski 1916; Mizuno et al. 1988): $$\begin{array}{lll}\frac{\text{d}\rho (m,t)}{\text{d}t}\hfill & =\hfill & -{\displaystyle {\int }_0^{\infty }}mK(m,m^{\prime })n(m,t)n(m^{\prime },t)\text{d}{m}^{\prime }\hfill \\ \hfill & \hfill & +\frac{1}{2}{\displaystyle {\int }_0^{m}m}K(m^{\prime }-m,m^{\prime })n(m^{\prime }-m,t)n(m^{\prime },t)\text{d}{m}^{\prime },\hfill \end{array}$$(B.1)

where m is the grain mass, n the density of grains, ρ = m × n the mass density of grains, and K(m, m′) the collision kernel (cm ⁻³ s⁻¹) between grains of mass m and m′.

B.1 Modelling of the size distribution

The grain size distribution is modelled by N_bins size bins, spread over the radius range [a₋ : a₊], or equivalently over the mass range [m₋ : m₊], assuming a specific density μ common to all grain cores. The lower and upper grain mass in bin k, m_k− and m_k+, respectively, follow a geometrical distribution in k of common ratio η (see Appendix A from Guillet et al. 2007): $$\eta ={\text{e}}^{-3\log (a_{+}/a_{-})/N_{\text{bins}}}<1,$$(B.2) $$\begin{array}{l}{m}_{k+}=m_{+}{\eta }^{k-1},\\ m_{k-}=m_{+}{\eta }^k.\end{array}$$

Bin 1 contains the largest grains and bin N contains the smallest grains, of the size distribution. We fix a₋ = a_N− = 5 nm, and a₊ = a₁₊ = 100 μm. Each bin k contains a mass density ρ_k of grain cores. At t = 0 bins where grain radius is larger than the MRN are empty. Silicate and carbon grains are mixed in a unique composite distribution of cores.

Following Mizuno et al. (1988), the mass density of grains in bins k − 1 and k + 1 is used to derive the spectral index β_k of the mass distribution within bin k: $$\frac{\text{d}n(m)}{\text{d}m}\propto m^{{\beta }_k},$$(B.5)

from which we can derive the average radius (a_k), cross-section (${\widehat{\sigma }}_k$) and mass (m_k) of grains in bin k. See Appendix A of Guillet et al. (2007) for more details on these calculations. The number of grains in bin k is simply n_k = ρ_k∕m_k.

B.2 Numerical implementation of the coagulation equation

To discretise the coagulation equation, we use the following procedure. We will consider collisions between grains of bin P (projectiles) and bin T (targets), P ≥ T implying that projectiles cannot be larger than targets. The coagulation rate between these two bins is $${\left(\frac{\text{d}n}{\text{d}t}\right)}_{\text{P},\text{T}}=n_{\text{T}}{n}_{\text{P}}{K}_{\text{P},\text{T}}=n_{\text{T}}{n}_{\text{P}}\pi {\left(a_{\text{T}}+a_{\text{P}}\right)}^2\Delta V_{\text{P},\text{T}},$$(B.6)

where ΔV_P,T is the rms relative velocity between the projectile and target grains.

Grain–grain coagulation leads to a transfer of mass from the projectile and target bins to other bins at the rate r. Taking into account that each bin k contains a size distribution of grains in the size range [m_k− : m_k+], the smallest coagulated grain, of mass m_P- + m_T-, must be transferred to the collector bin of index $$C=1+\text{int}\left(\frac{\log \left(\frac{m_{\text{P-}}+m_{\text{T-}}}{m_{+}}\right)}{\log \eta }\right),$$(B.7)

while the largest coagulated grains, of mass m_P+ + m_T+, must be transferred to collector bins C − 1 according to the same equation. There will be no transfer of mass from bins P and T to any other bins.

The fraction f of grain mass transferred to collector bin C − 1, and therefore 1 − f to collector C, is calculated as follows. We assume that all grains from the projectile bin are identical to the average grain of mass m_P. With respect to the size of the larger grains from the target bin, this is most of the time a reasonable assumption. If m_P + m_T- ≥ m_C+, all coagulated grains will be transferred to collector bin C − 1 (and therefore none to collector bin C): f = 1. If m_P + m_T+ ≤ m_C+, all coagulated grains will be transferred to collector bin C (and none to collector bin C − 1): f = 0. In the general case, the fraction f of coagulated grains transferred to collector bin C − 1 is $$f={\displaystyle {\int }_{m_{\text{P-}}}^{m_{\text{P+}}}}{\displaystyle {\int }_{m_{\text{C+}}-m_{\text{P}}}^{m_{\text{T+}}}(m+m^{\prime })}K(m,m^{\prime })m\text{{^}{\beta }_{\text{P}}\text{}}{m^{\prime }}^{{\beta }_{\text{T}}}\text{d}m\text{d}{m}^{\prime }{\displaystyle {\int }_{m_{\text{P-}}}^{m_{\text{P+}}}}{\displaystyle {\int }_{m_{\text{T-}}}^{m_{\text{T+}}}(m+m^{\prime })}K(m,m^{\prime })m\text{{^}{\beta }_{\text{P}}\text{}}{m^{\prime }}^{{\beta }_{\text{T}}}\text{d}m\text{d}{m}^{\prime }\text{}}.$$(B.8)

There exists analytical expressions for f whenever K is a function of powers of m and m′, which is usually the case.

The mass transfer rates from the projectile and target bins to the collector bins are therefore: $$\begin{array}{lll}\frac{\text{d}{\rho }_{\text{P}}}{\text{d}t}\hfill & =\hfill & -m_{\text{P}}{\left(\frac{\text{d}n}{\text{d}t}\right)}_{\text{P},\text{T}},\hfill \\ \frac{\text{d}{\rho }_{\text{T}}}{\text{d}t}\hfill & =\hfill & -m_{\text{T}}{\left(\frac{\text{d}n}{\text{d}t}\right)}_{\text{P},\text{T}},\hfill \\ \frac{\text{d}{\rho }_{\text{C}}}{\text{d}t}\hfill & =\hfill & \left(f-1\right)\left(m_{\text{P}}+m_{\text{T}}\right){\left(\frac{\text{d}n}{\text{d}t}\right)}_{\text{P},\text{T}},\hfill \\ \frac{\text{d}{\rho }_{\text{P}}}{\text{d}t}\hfill & =\hfill & f\left(m_{\text{P}}+m_{\text{T}}\right){\left(\frac{\text{d}n}{\text{d}t}\right)}_{\text{P},\text{T}}.\hfill \end{array}$$

The total grain mass is conserved. This procedure is reproduced for any pair of projectile bin/target bin where P ≥ T, that is when the projectile is smaller than the target, or of the same size as.

B.3 Tests of the coagulation algorithm

The Smoluchowski equation possesses self-similar solutions for the constant (K(m, m′) = 2), additive (K(m, m′) = m + m′), and multiplicative (K(m, m′) = m × m′) kernels (e.g. Menon & Pego 2003).

In log–log plots of the mass density, this self similarity appears clearly $$\begin{array}{lll}\frac{\text{d}\log {\rho }^C(m,t)}{\text{d}\log m}\hfill & =\hfill & {\left(\frac{m}{t}\right)}^2{\text{e}}^{-m/t},\hfill \\ \frac{\text{d}\log {\rho }^A(m,t)}{\text{d}\log m}\hfill & =\hfill & \sqrt{\frac{m{\text{e}}^{-2t}}{2\pi }}{\text{e}}^{-m{\text{e}}^{-2t}/2}.\hfill \end{array}$$

At t = 0, the self-similar solutions express $$\begin{array}{lll}\frac{\text{d}{\log }^C\rho (m,t=0)}{\text{d}m}\hfill & =\hfill & m{\text{e}}^{-m},\hfill \\ \frac{\text{d}{\log }^A\rho (m,t=0)}{\text{d}m}\hfill & =\hfill & \frac{{\text{e}}^{-m/2}}{\sqrt{2\pi m}},\hfill \end{array}$$

which can be integrated analytically to yield the mass density in each bin k: $${\rho }_k^C={\displaystyle {\int }_{m_{k-}}^{m_{k+}}m}{\text{e}}^{-m}\text{d}m=\left(1+m_{k-}\right){\text{e}}^{-m_{k-}}-\left(1+m_{k+}\right){\text{e}}^{-m_{k+}},$$(B.17)

and $${\rho }_k^A={\displaystyle {\int }_{m_{k-}}^{m_{k+}}\sqrt{\frac{{\text{e}}^{-m/2}}{2\pi m}}}\text{d}m=\text{erf}\left(\sqrt{\frac{m_{k+}}{2}}\right)-\text{erf}\left(\sqrt{\frac{m_{k-}}{2}}\right),$$(B.18)

where erf is the error function.

We test our coagulation algorithm by comparing our numerical results with analytical solutions for the constant and additive kernels (Fig. B.1). We use 35 bins, equally spaced in logarithmic scale on the range of adimensioned masses [4 × 10⁻⁶ : 4 × 10⁸] (η = 0.4)⁵. In both cases, our algorithm is able to reproduce analytical results over a 10⁶ -fold increase in the average grain mass (a factor of 100 in radius). The power-law of the smallest grains is well followed-up, which is particularly important in our study. Because of numerical errors, the exponential decay at large grain radius is overestimated in our simulation. This is however a small effect, which furthermore will have no impact on the ionisation and MHD properties of the gas which are dominated by the small grain population. B.4:

Fig. B.1 Comparison between the analytical (solid) and numerical (dashed) solutions to the Smoluchowski equation, at various adimensioned time steps, in the case of a constant (left) and additive (right) kernels, with 35 size bins equally spaced in logarithmic scale over a range [4 × 10−6 : 4 × 108] in mass (η = 0.4). Diamonds indicate the mean mass and mass density in each bin. Grain growth is linear with time for the constant kernel, and exponential with time for the additive kernel.

Fig. B.2 Check for the convergence of our coagulation/accretion algorithm at nH = 1012 cm−3, varying the number of logarithmic bins in the 5 nm–100 μm range.

B.4 Convergence tests

We test for the convergence of our coagulation algorithm using our standard model, but ignoring grain charge (which does not affect grain coagulation rates in our model). Figure B.2 shows that convergence is obtained as soon as η < 0.4, as for the self-similar solutions to the Smoluchowski equation (see Appendix. B.3). We will therefore use a total of 33 bins:13 bins for the MRN (0.005 μm < a < 0.25 μm) and 20 empty bins for 0.25 < a < 100 μm.

Appendix C Coagulation rates with a skewed Maxwellian velocity distribution

Following the approach by Flower et al. (2005), we consider the dynamics of particles having velocities v_x, v_y and v_z along $\widehat{x}$, ŷ and ẑ in the lab frame, respectively, each characterised by a normal distribution of zero expectation and same variance ΔV ²∕3. We observe the motion of those particles from a reference frame moving at a velocity μ >0 along the ẑ axis.

Defining $$\alpha \equiv \sqrt{\frac{3}{2}}\frac{1}{\Delta V},$$(C.1)

the probability distribution functions along for v_x, v_y, and v_z are $f(v_x)=\frac{\alpha }{\sqrt{\pi }}{\text{e}}^{-{\alpha }^2{v}_x^2}$, $f(v_y)=\frac{\alpha }{\sqrt{\pi }}{\text{e}}^{-{\alpha }^2{v}_y^2}$, and $f(v_z)=\frac{\alpha }{\sqrt{\pi }}{\text{e}}^{-{\alpha }^2{\left(v_z-\mu \right)}^2}$, respectively.

The probability distribution of v, the velocity vector measured in the frame of the observer, is: $$f(v)=f(v_x,v_y,v_z)={\left(\frac{\alpha }{\sqrt{\pi }}\right)}^3{\text{e}}^{-{\alpha }^2\left(v^2+{\mu }^2\right)}{\text{e}}^{-2{\alpha }^2{v}_z\mu }.$$(C.2)

If v is the velocity amplitude measured in the frame of the observer, and (θ, ϕ) the spherical coordinates of the velocity vector, then by definition v_z = vcosθ. To obtain the probability distribution of the velocity amplitude, f(v), we integrate f(v) over the spherical coordinates (θ, ϕ): $$\begin{array}{lll}f(v)\hfill & =\hfill & {\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }{v}^2}}f\left(v\right)\sin \theta \text{d}\theta \text{d}\varphi ,\hfill \\ \hfill & =\hfill & 2\pi {\left(\frac{\alpha }{\sqrt{\pi }}\right)}^3{\text{e}}^{-{\alpha }^2\left(v^2+{\mu }^2\right)}{\displaystyle {\int }_0^{\pi }{v}^2}{\text{e}}^{-2{\alpha }^{2}v\mu \cos \theta }\sin \theta \text{d}\theta ,\hfill \\ \hfill & =\hfill & 2\pi {\left(\frac{\alpha }{\sqrt{\pi }}\right)}^3{\text{e}}^{-{\alpha }^2\left(v^2+{\mu }^2\right)}{\left[\frac{v}{2{\alpha }^2\mu }{\text{e}}^{-2{\alpha }^{2}v\mu \cos \theta }\right]}_0^{\pi },\hfill \\ \hfill & =\hfill & \frac{1}{\sqrt{\pi }}\frac{\alpha }{\mu }v{\text{e}}^{-{\alpha }^2\left(v^2+{\mu }^2\right)}\left({\text{e}}^{2{\alpha }^{2}v\mu }-{\text{e}}^{-2{\alpha }^{2}v\mu }\right),\hfill \\ \hfill & =\hfill & \frac{1}{\sqrt{\pi }}\frac{\alpha }{\mu }v\left({\text{e}}^{-{\alpha }^2{\left(v-\mu \right)}^2}-{\text{e}}^{-{\alpha }^2{\left(v+\mu \right)}^2}\right).\hfill \end{array}$$

The mean collisional velocity among those collisions leading to grain–grain coagulation is $${\displaystyle {\int }_0^{\mathcal{V}}v}f(v)\text{d}v=\frac{1}{\sqrt{\pi }}\frac{\alpha }{\mu }{\displaystyle {\int }_0^{\mathcal{V}}{v}^2}\left({\text{e}}^{-{\alpha }^2{\left(v-\mu \right)}^2}-{\text{e}}^{-{\alpha }^2{\left(v+\mu \right)}^2}\right).$$(C.8)

We first calculate: $$\begin{array}{lll}{\displaystyle {\int }_0^{\mathcal{V}}{v}^2}{\text{e}}^{-{\alpha }^2{(v+\mu )}^2}\text{d}v\hfill & =\hfill & {\displaystyle {\int }_{\mu }^{\mathcal{V}+\mu }\left(v^2+{\mu }^2-2v\mu \right)}{\text{e}}^{-{\alpha }^2{v}^2}\text{d}v.\hfill \end{array}$$(C.9)

We have: $$\begin{array}{l}{\displaystyle \int {\text{e}}^{-{\alpha }^2{x}^2}}\text{d}x=\frac{\sqrt{\pi }}{2}\frac{\text{erf}\left(\alpha x\right)}{\alpha },\\ {\displaystyle \int x}{\text{e}}^{-{\alpha }^2{x}^2}\text{d}x=-\frac{1}{2{\alpha }^2}{\text{e}}^{-{\alpha }^2{x}^2},\\ {\displaystyle \int x^2}{\text{e}}^{-{\alpha }^2{x}^2}\text{d}x=\frac{1}{2{\alpha }^2}\left(\frac{\sqrt{\pi }}{2}\frac{\text{erf}\left(\alpha x\right)}{\alpha }-x{\text{e}}^{-{\alpha }^2{x}^2}\right).\end{array}$$

Then $$\begin{array}{lll}{\displaystyle {\int }_0^{\mathcal{V}}{v}^2}{\text{e}}^{-{\alpha }^2{(v+\mu )}^2}\text{d}v\hfill & =\hfill & \frac{1}{2{\alpha }^2}{\left[\left(2\mu -v\right){\text{e}}^{-{\alpha }^2{v}^2}+\left(2{\alpha }^2{\mu }^2+1\right)\frac{\sqrt{\pi }}{2}\frac{\text{erf}\left(\alpha v\right)}{\alpha }\right]}_{\mu }^{\mathcal{V}+\mu }\hfill \\ \hfill & =\hfill & \frac{1}{2{\alpha }^2}\left(\left(\mu -\mathcal{V}\right){\text{e}}^{-{\alpha }^2{(\mathcal{V}+\mu )}^2}-\mu {\text{e}}^{-{\alpha }^2{\mu }^2}+\left(2{\alpha }^2{\mu }^2+1\right)\frac{\sqrt{\pi }}{2\alpha }\left[\text{erf}\left(\alpha \mathcal{V}+\alpha \mu \right)-\text{erf}\left(\alpha \mu \right)\right]\right).\hfill \end{array}$$(C.13)

Similarly, $$\begin{array}{lll}{\displaystyle {\int }_0^{\mathcal{V}}{v}^2}{\text{e}}^{-{\alpha }^2{(v-\mu )}^2}\text{d}v\hfill & =\hfill & \frac{1}{2{\alpha }^2}\left(\left(-\mu -\mathcal{V}\right){\text{e}}^{-{\alpha }^2{(\mathcal{V}-\mu )}^2}+\mu {\text{e}}^{-{\alpha }^2{\mu }^2}+\left(2{\alpha }^2{\mu }^2+1\right)\frac{\sqrt{\pi }}{2{\alpha }^2}\left[\text{erf}\left(\alpha \mathcal{V}-\alpha \mu \right)+\text{erf}\left(\alpha \mu \right)\right]\right).\hfill \end{array}$$(C.14)

As a consequence, $$\begin{array}{l}{\displaystyle {\int }_0^{\mathcal{V}}v}f(v)\text{d}v=\frac{1}{\sqrt{\pi }}\frac{\alpha }{\mu }\frac{1}{2{\alpha }^2}\\ \times \left[-\left(\mu +\mathcal{V}\right){\text{e}}^{-{\alpha }^2{(\mathcal{V}-\mu )}^2}-\left(\mu -\mathcal{V}\right){\text{e}}^{-{\alpha }^2{(\mathcal{V}+\mu )}^2}+2\mu {\text{e}}^{-{\alpha }^2{\mu }^2}+\left(2{\alpha }^2{\mu }^2+1\right)\frac{\sqrt{\pi }}{2\alpha }\left(\text{erf}\left(\alpha \mathcal{V}-\alpha \mu \right)-\text{erf}\left(\alpha \mathcal{V}+\alpha \mu \right)+2\text{erf}\left(\alpha \mu \right)\right)\right]\\ =\frac{1}{\alpha \sqrt{\pi }}\left[{\text{e}}^{-{\alpha }^2{\mu }^2}-\frac{1}{2}\left({\text{e}}^{-{\alpha }^2{(\mathcal{V}-\mu )}^2}\left[1+\frac{\mathcal{V}}{\mu }\right]+{\text{e}}^{-{\alpha }^2{(\mathcal{V}+\mu )}^2}\left[1-\frac{\mathcal{V}}{\mu }\right]\right)\right)+\left(\mu +\frac{1}{2{\alpha }^2\mu }\right)\left(\text{erf}\left(\alpha \mu \right)-\frac{\text{erf}\left(\alpha \mathcal{V}+\alpha \mu \right)-\text{erf}\left(\alpha \mathcal{V}-\alpha \mu \right)}{2}\right]\\ =\frac{1}{\alpha }\left[\frac{1}{\sqrt{\pi }}\left({\text{e}}^{-{\xi }^2}-\frac{1}{2}\left({\text{e}}^{-{\left(\chi -\xi \right)}^2}\left[1+\frac{\chi }{\xi }\right]+{\text{e}}^{-{\left(\chi +\xi \right)}^2}\left[1-\frac{\chi }{\xi }\right]\right)\right)+\left(\xi +\frac{1}{2\xi }\right)h(\chi ,\xi )\right],\\ \text{where}\xi =\alpha \mu ,\chi =\alpha \mathcal{V},\text{and}h(\chi ,\xi )=\text{erf}\left(\xi \right)-\left(\text{erf}\left(\chi +\text{}\xi \right)-\text{erf}\left(\chi -\xi \right)\right)/2.\end{array}$$

Here are the two asymptotic cases:

No drift

When μ (and ξ) tend towards zero, h(χ, ξ)∕ξ tends towards erf′(0) −erf′(χ) (where erf′ is the first derivative of erf), that is towards $\frac{2}{\sqrt{\pi }}\left(1-{\text{e}}^{-{\chi }^2}\right)$. We get: $$\begin{array}{lll}{\displaystyle {\int }_0^{\mathcal{V}}v}f(v)\text{d}v\hfill & =\hfill & \frac{1}{\alpha }\left[\frac{1}{\sqrt{\pi }}\left(1-\left(1+2{\chi }^2\right){\text{e}}^{-{\chi }^2}\right)+\frac{1}{\sqrt{\pi }}\left(1-{\text{e}}^{-{\chi }^2}\right)\right]\hfill \\ \hfill & =\hfill & \frac{2}{\alpha \sqrt{\pi }}\left(1-\left(1+{\chi }^2\right){\text{e}}^{-{\chi }^2}\right)\hfill \\ \hfill & =\hfill & \sqrt{\frac{8}{3\pi }}\Delta V\left(1-\left(1+{\chi }^2\right){\text{e}}^{-{\chi }^2}\right),\hfill \end{array}$$(C.16)