© ESO, 2016

1. Introduction

It is well established today that stars form within turbulent complexes that range from tens to hundreds of parsecs, that harbour thousands of solar masses of cold gas: the molecular clouds (e.g. McKee & Ostriker 2007,for a review). Turbulence within these clouds generates over-dense regions where the balance between the thermal and non-thermal (outward) pressure and the gravitational (inward) force is disrupted, causing the gas to collapse. The mass distribution of these dense cores shows many similarities with the stellar mass function (see e.g. Könyves et al. 2010), making them the most obvious stellar progenitors. The details of the formation process are still actively debated however. The (initially optically thin) core is believed to first contract isothermally as the compression heating is lost through radiation, until the density is high enough to render the cooling ineffective. This leads to the formation of a hydrostatic body, known as the first Larson core (Larson 1969), which accretes material from the surrounding envelope. The sustained increase in mass, density, and temperature eventually triggers the dissociation of H₂ molecules above 2000 K. This leads to the second phase of collapse because of the endothermic nature of the dissociation process. The collapse ceases when most or all of the H₂ molecules have been dissociated, at which point a second much more dense and compact hydrostatic core (Larson’s second core) is formed at the centre (Larson 1969; Masunaga & Inutsuka 2000; Vaytet et al. 2013). The temperature inside the second core continues to rise until the nuclear reactions are ignited: the young star is born.

This problem entails many physical processes over a very wide range of spatial scales, making any modelling of the entire process challenging; the cloud core has an initial radius of ~10⁴ astronomical units (AU) with an average density of a few 10³ cm^-3, while the protostellar core measures only 10^-3 AU at birth, with densities over 10²⁰ cm^-3. Describing and understanding the physical processes at work involves an intricate interplay between large-scale environmental factors, which regulate the supply of mass, angular momentum, and magnetic flux, and small-scale processes, which control the evolution and dynamics in protostellar systems. An accurate description of a large number of complex physical mechanisms, involving the magnetic field, gravity, radiative transfer, time-dependent chemistry, and dust physics, is necessary to derive realistic models of the global star formation process.

Recently, various studies have been devoted to the role magnetic fields play in collapsing systems and the effect on the transport of angular momentum. The first studies were based on ideal magnetohydrodynamics (MHD), implying numerical instead of physical diffusivity. The results concerning the disk formation are thus of dubious validity (Allen et al. 2003; Matsumoto & Tomisaka 2004; Galli et al. 2006; Price & Bate 2007; Hennebelle & Teyssier 2008b; Hennebelle & Fromang 2008a; Commerçon et al. 2010; Masson et al. 2016). The reconnection and diffusion of the field are essential components of the magnetic braking process. Using the framework of non-ideal MHD allows describing the diffusivity and the magnetic properties of the charged fluid accurately (Machida et al. 2006; Duffin & Pudritz 2008; Mellon & Li 2009; Machida & Matsumoto 2011). This requires the accurate calculations of magnetic resistivities, however. This remains challenging as they depend on many factors such as temperature, density, chemical abundances and magnetisation.

The interstellar gas from which the cores (and thus stars) form is essentially composed of hydrogen (~74% by mass), helium (~25%), and heavier elements (~1%) such as carbon, oxygen and heavy metals (Na, Fe, etc.), whose respective abundances are determined by an elaborate set of chemical reactions. The complexity is even increased by of the formation of dust grains, which arise from the aggregation of several molecules and can reach micrometer sizes (Kunz & Mouschovias 2009). Not only do grains react with other elements, they also carry electric charges and play the role of catalysts for other chemical reactions. The grain ionisation rates, or even the chemical reaction rates themselves, remain poorly constrained in environments typical of prestellar core collapse. The relative abundances and the degree of ionisation of the grains and molecules define the magnetic resistivities that regulate the dissipation of the magnetic flux through ambipolar and Ohmic diffusion and the Hall effect (Krasnopolsky et al. 2010; Machida et al. 2006; Li et al. 2011). These dissipation processes have a fundamental role in the formation and evolution of the first and second Larson cores and their associated disks and outflows. In particular, the interplay between flux-freezing and condensation of the global angular momentum strongly depends on whether ambipolar diffusion, Hall effect, or Ohmic diffusion is the dominant process. To compute the accurate resisitivities under the conditions typical of collapsing molecular clouds, we have developed a relevant reduced chemical network. This allows us to test the effect of various parameters, such as the equilibrium abundances of chemical species, the evolution of grains of different sizes, and the cosmic-ray ionisation rate on the various non-ideal MHD diffusion coefficients. This network significantly extends these and other previous prescriptions (Umebayashi & Nakano 1990; Nishi et al. 1991; Wardle & Ng 1999; Nakano et al. 2002; Wardle 2007; Kunz & Mouschovias 2009; Ilgner & Nelson 2006; Bai 2011) by including new pieces of physics that are necessary to precisely describe the chemical evolution of the interstellar gas.

The paper is organised as follows. In Sect. 2, we outline the main processes and physical ingredients relevant to star formation. In Sect. 3, we present our chemical network and the numerical method used to solve the reaction equations. The resulting resistivities, the effect of each ingredient on the chemical equilibrium of the gas, and the consequences for star-forming systems are discussed in Sect. 4, while Sect. 5 is devoted to the conclusion.

2. Physics of collapsing prestellar cores

2.1. Magnetic resistivities

The complete induction equation in non-ideal MHD reads $\begin{array}{ccc}\frac{\mathit{\partial }{B}}{\mathit{\partial t}}\mathrm{-}\mathrm{\nabla }\mathrm{\times }\mathrm{(}{v}\mathrm{\times }{B}\mathrm{)}\mathrm{=}& & \mathrm{-}\frac{{\mathit{c}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }}\mathrm{\nabla }\mathrm{\times }[{\mathit{\eta }}_{\mathrm{\Omega }}\mathrm{(}\mathrm{\nabla }\mathrm{\times }{B}\mathrm{)}\\ & & \mathrm{+}{\mathit{\eta }}_{\mathrm{H}}\left\{\mathrm{(}\mathrm{\nabla }\mathrm{\times }{B}\mathrm{)}\mathrm{\times }\frac{{B}}{\mathrm{\left|}\mathrm{\right|}{B}\mathrm{\left|}\mathrm{\right|}}\right\}\\ & & \end{array}$(1)where B denotes the magnetic field, | | B | | its L₂ norm, and v the fluid velocity. η_Ω, η_H and η_AD represent the Ohmic, Hall, and ambipolar resistivities, respectively. These dissipative terms account for collisions between neutral and charged species. They strongly depend on the chemical equilibrium and thus on the thermodynamic conditions. In contrast, ideal MHD considers the fluid as a mixture of perfectly coupled charged fluids, which corresponds to η_Ω = η_H = η_AD = 0, ensuring flux-freezing of the magnetic field.

The resistivities are defined in terms of the conductivities of the gas-dust mixture as $\begin{array}{ccc}{\mathit{\eta }}_{\mathrm{\Omega }}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{\parallel }}}\mathit{,}\\ {\mathit{\eta }}_{\mathrm{H}}& \mathrm{=}& \frac{{\mathit{\sigma }}_{\mathrm{H}}}{{\mathit{\sigma }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{H}}^{\mathrm{2}}}\mathit{,}\\ {\mathit{\eta }}_{\mathrm{AD}}& \mathrm{=}& \frac{{\mathit{\sigma }}_{\mathrm{\perp }}}{{\mathit{\sigma }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{H}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{\parallel }}}\mathit{,}\end{array}$where the parallel, perpendicular, and Hall conductivites are in turn expressed, respectively, as $\begin{array}{ccc}{\mathit{\sigma }}_{\mathrm{\parallel }}& \mathrm{=}& \sum _{\mathit{i}}{\mathit{\sigma }}_{\mathit{i}}\mathit{,}\\ {\mathit{\sigma }}_{\mathrm{\perp }}& \mathrm{=}& \sum _{\mathit{i}}\frac{{\mathit{\sigma }}_{\mathit{i}}}{\mathrm{1}\mathrm{+}\mathrm{(}{\mathit{\omega }}_{\mathit{i}}{\mathit{\tau }}_{\mathit{i}\mathrm{n}}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}\\ {\mathit{\sigma }}_{\mathrm{H}}& \mathrm{=}& \mathrm{-}\sum _{\mathit{i}}\frac{{\mathit{\sigma }}_{\mathit{i}}{\mathit{\omega }}_{\mathit{i}}{\mathit{\tau }}_{\mathit{i}\mathrm{n}}}{\mathrm{1}\mathrm{+}\mathrm{(}{\mathit{\omega }}_{\mathit{i}}{\mathit{\tau }}_{\mathit{i}\mathrm{n}}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}\end{array}$with ${\mathit{\sigma }}_{\mathit{i}}\mathrm{=}\frac{{\mathit{n}}_{\mathit{i}}{\mathit{q}}_{\mathit{i}}^{\mathrm{2}}{\mathit{\tau }}_{\mathit{i}\mathrm{n}}}{{\mathit{m}}_{\mathit{i}}}$, ${\mathit{\omega }}_{\mathit{i}}\mathrm{=}\frac{{\mathit{q}}_{\mathit{i}}\mathit{B}}{{\mathit{m}}_{\mathit{i}}\mathit{c}}$ the cyclotron frequency, and ${\mathit{\tau }}_{\mathit{i}\mathrm{n}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{a}}_{\mathit{i}\mathrm{He}}}\frac{{\mathit{m}}_{\mathit{i}}\mathrm{+}{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}}}{{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}}}\frac{\mathrm{1}}{{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{⟨}{\mathit{\sigma }}_{\mathrm{coll}}\mathit{w}{\mathrm{⟩}}_{\mathit{i}}}\mathrm{·}$(8)Here i stands for any charged particle of charge q_i, with a charge particle number-density n_i, m denotes the mass of a particle, and c is the speed of light. The factor a_iHe 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). ⟨ σ_collw ⟩ _i is the rate constant for collisions between a particle i and H₂ molecules (Desch & Mouschovias 2001; Pinto & Galli 2008). The values are taken from Pinto & Galli (2008). We note that these rates were calculated in a three-fluid formalism, and we use them in a multifluid context, but their dependance with temperature is a better approximation than the Langevin model Pinto & Galli (2008). ω_i can be negative, depending on the charge q_i. This is important when calculating the conductivities, as the Hall conductivity σ_H can be either positive or negative. This implies that the Hall resistivity can be negative, since it is of the same sign as σ_H.

2.2. Density and temperature range

Our chemical network was used to compute conductivities for an entire density-temperature range spanning 10² cm^-3<n< 10²⁶ cm^-3 and 10 K <T< 10⁵ K. This covers typical (parent) molecular cloud conditions (low temperature and density) as well as the interior of first and second Larson cores, where densities and temperatures have increased by many orders of magnitude, as a result of the strong gravitational compression of the gas. However, for clarity as well as for comparison with previous studies, the majority of the results presented in this work are presented only as a function of density. To mimic the natural temperature evolution as a function of density in a collapsing dense core (as opposed to using a single constant temperature), we used the piecewise barotropic equation of state (EOS) from Machida et al. (2006)$\mathit{T}\mathrm{=}{\mathit{T}}_{\mathrm{0}}\sqrt{\mathrm{1}\mathrm{+}{\left(\frac{\mathit{n}}{{\mathit{n}}_{\mathrm{1}}}\right)}^{\mathrm{2}{\mathit{g}}_{\mathrm{1}}}}{\left(\mathrm{1}\mathrm{+}\left(\frac{\mathit{n}}{{\mathit{n}}_{\mathrm{2}}}\right)\right)}^{{\mathit{g}}_{\mathrm{2}}}{\left(\mathrm{1}\mathrm{+}\left(\frac{\mathit{n}}{{\mathit{n}}_{\mathrm{3}}}\right)\right)}^{{\mathit{g}}_{\mathrm{3}}}\mathit{,}$(9)with n the total density and $\begin{array}{c}\\ {\mathit{n}}_{\mathrm{1}}& \mathrm{=}& {\mathrm{10}}^{\mathrm{11}}{\mathrm{cm}}^{-3}\mathrm{;}& & {\mathit{g}}_{\mathrm{1}}& \mathrm{=}& \mathrm{0.4}\mathrm{;}\\ {\mathit{n}}_{\mathrm{2}}& \mathrm{=}& {\mathrm{10}}^{\mathrm{16}}{\mathrm{cm}}^{-3}\mathrm{;}& & {\mathit{g}}_{\mathrm{2}}& \mathrm{=}& \mathrm{-}\mathrm{0.3}\mathrm{;}\\ {\mathit{n}}_{\mathrm{3}}& \mathrm{=}& {\mathrm{10}}^{\mathrm{21}}{\mathrm{cm}}^{-3}\mathrm{;}& & {\mathit{g}}_{\mathrm{3}}& \mathrm{=}& \mathrm{0.56667}\mathit{,}\end{array}$(10)and T₀ = 10 K. This EOS successively characterises the isothermal phase of the collapse, the adiabatic phase during the first Larson core evolution, the second collapse and the second Larson core evolution. It is important to use such a typical EOS to present our results because the temperature increase during the collapse strongly influences on the chemistry of the grains (see Sect. 2.4.2) and on the ionisation of potassium (Sect. 3.1). The EOS is represented in Fig. 1 with our magnetic field prescription (see next paragraph).

2.3. Magnetic field

Computing the resistivities also requires a knowledge of the magnetic field intensity. To present our magnetic resistivities as a function of density alone, we assumed that the magnetic intensity scales as $\mathit{B}\mathrm{\left(}{\mathit{n}}_{\mathrm{H}}\mathrm{\right)}\mathrm{=}\mathrm{1.43}\mathrm{\times }{\mathrm{10}}^{-7}\sqrt{{\mathit{n}}_{\mathrm{H}}}$ (Li et al. 2011), which corresponds to magnetic flux conservation (flux-freezing approximation). It is represented in Fig. 1. Although the scaling in real collapse calculations is more complicated (see e.g. Masson et al. 2016), the magnetic field here is only used to illustrate the behaviour of the resistivities. It has no influence on the chemistry. In detailed prestellar collapse calculations, the magnetic field evolution is properly calculated and the resistivities are computed consistently during the collapse (see Masson et al. 2016).

2.4. Dust grain model

Grains can be the main charge carriers and thus need to be accurately described, both in size and number density, because these two quantities determine the surface area available for chemical reactions. Furthermore, grain evaporation is a process of prime importance in the present context because it occurs at temperatures close to the second collapse, after the first core formation.

2.4.1. Grain size

The grain reaction rates and the conductivity of the dust-gas mixture depend on the grain cross section. In our calculations, we included a power-law grain-size distribution by considering a finite number N_bins of size bins with equal widths in log space. We define a minimum and maximum grain size a_min and a_max, respectively, with a number density of grains of radius between a and a + da$\mathrm{d}{\mathit{n}}_{\mathrm{g}\mathit{,}\mathrm{tot}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\mathrm{𝒞}\hspace{0.17em}{\mathit{a}}^{\mathit{\lambda }}\hspace{0.17em}\mathrm{d}\mathit{a,}$(11)where the subscript g denotes the grains and is a normalisation constant. Unless otherwise stated, we used the standard MRN distribution with λ = −3.5 (Mathis et al. 1977) throughout this work. For the sake of generality, however, we write here the equations for any power-law index λ. Following Kunz & Mouschovias (2009), we chose for the minimum and maximum radii a_min = 0.0181 μm and a_max = 0.9049 μm. Each size bin is defined by a lower and upper radius that give the number density and size for the αth bin (α = 1,2,...,N_bins) $\begin{array}{ccc}{\mathit{n}}_{\mathrm{g}\mathit{,\alpha }}& \mathrm{=}& {\mathit{n}}_{\mathrm{g}\mathit{,}\mathrm{tot}}{\mathit{\xi }}^{\frac{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{+}\mathrm{1}\mathrm{)}\mathit{\alpha }}{{\mathit{N}}_{\mathrm{bins}}}}\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\xi }}^{\frac{\mathit{\lambda }\mathrm{+}\mathrm{1}}{{\mathit{N}}_{\mathrm{bins}}}}}{{\mathit{\xi }}^{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{-}\mathrm{1}}\right)\mathit{,}\\ {\mathit{a}}_{\mathit{\alpha }}& \mathrm{=}& {\mathit{a}}_{\min }{\mathit{\xi }}^{\frac{\mathrm{-}\mathit{\alpha }}{{\mathit{N}}_{\mathrm{bins}}}}{\left[\left(\frac{\mathit{\lambda }\mathrm{+}\mathrm{1}}{\mathit{\lambda }\mathrm{+}\mathrm{3}}\right)\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\xi }}^{\frac{\mathit{\lambda }\mathrm{+}\mathrm{3}}{{\mathit{N}}_{\mathrm{bins}}}}}{\mathrm{1}\mathrm{-}{\mathit{\xi }}^{\frac{\mathit{\lambda }\mathrm{+}\mathrm{1}}{{\mathit{N}}_{\mathrm{bins}}}}}\right)\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathit{,}\end{array}$where $\mathit{\xi }\mathrm{=}\frac{{\mathit{a}}_{\min }}{{\mathit{a}}_{\max }}$. The total number density of dust n_g,tot is determined by constraining the total grain mass density in the size distribution to be ${\mathit{\rho }}_{\mathrm{g}\mathit{,}\mathrm{tot}}\mathrm{=}{\mathrm{\int }}_{{\mathit{a}}_{\min }}^{{\mathit{a}}_{\max }}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{\rho }}_{\mathit{s}}{\mathit{a}}^{\mathrm{3}}\mathrm{𝒞}{\mathit{a}}^{\mathit{\lambda }}\mathrm{d}\mathit{a,}$(14)which yields ${\mathit{n}}_{\mathrm{g}\mathit{,}\mathrm{tot}}\mathrm{=}\left(\frac{{\mathit{\rho }}_{\mathrm{g}\mathit{,}\mathrm{tot}}}{\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{\rho }}_{\mathit{s}}{\mathit{a}}_{\min }^{\mathrm{3}}}\right)\left(\frac{\mathit{\lambda }\mathrm{+}\mathrm{4}}{\mathit{\lambda }\mathrm{+}\mathrm{1}}\right)\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\xi }}^{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{+}\mathrm{1}\mathrm{)}}}{{\mathit{\xi }}^{\mathit{\lambda }\mathrm{+}\mathrm{4}}\mathrm{-}\mathrm{1}}\right){\mathit{\xi }}^{\mathit{\lambda }\mathrm{+}\mathrm{4}}\mathit{.}$(15)Finally, we constrained the factor by choosing the total grain density ρ_g,tot to obtain the same total surface area as a fiducial uniform distribution with a grain size a₀ (see Kunz & Mouschovias 2009). Because ${\mathit{n}}_{\mathrm{g}\mathit{,}\mathrm{tot}}\mathrm{=}\frac{{\mathit{a}}_{\max }^{\mathit{\lambda }\mathrm{+}\mathrm{1}}\mathrm{𝒞}}{\mathit{\lambda }\mathrm{+}\mathrm{1}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\xi }}^{\mathit{\lambda }\mathrm{+}\mathrm{1}}\mathrm{)}\mathit{,}$(16)this yields ${\mathit{\rho }}_{\mathrm{g}\mathit{,}\mathrm{tot}}\mathrm{=}{\mathit{\rho }}_{\mathrm{g}\mathit{,}\mathrm{tot}}^{\mathrm{fiducial}}\left(\frac{{\mathit{a}}_{\min }}{{\mathit{a}}_{\mathrm{0}}}\right)\left(\frac{\mathit{\lambda }\mathrm{+}\mathrm{3}}{\mathit{\lambda }\mathrm{+}\mathrm{4}}\right)\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\xi }}^{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{+}\mathrm{4}\mathrm{)}}}{\mathrm{1}\mathrm{-}{\mathit{\xi }}^{\mathrm{-}\mathrm{(}\mathit{\lambda }\mathrm{+}\mathrm{3}\mathrm{)}}}\right)\mathit{,}$(17)where, for a dust to gas ratio of 1%, ${\mathit{\rho }}_{\mathrm{g}\mathit{,}\mathrm{tot}}^{\mathrm{fiducial}}\mathrm{=}\mathrm{0.01}{\mathit{\rho }}_{\mathrm{n}\mathit{,}\mathrm{tot}}$ (ρ_n,tot represents the density of neutral gas). In the MRN distribution, small grains vastly outnumber the large grains, contributing dominantly to the total surface area. The total grain density is then ${\mathit{\rho }}_{\mathrm{g}\mathit{,}\mathrm{tot}}^{\mathrm{\left(}\mathrm{MRN}\mathrm{\right)}}\mathrm{=}\mathrm{0.0341}{\mathit{\rho }}_{\mathrm{n}\mathit{,}\mathrm{tot}}$.

2.4.2. Grains at high temperature

Grain evaporation.

The three main grain constituents are carbon (essentially amorphous), silicates (here represented by the molecule (MgFe)SiO₄) and aluminium oxyde (Al₂O₃). For the sake of simplicity, we assumed that each grain is composed of only one of these three materials, instead of considering a layered structure, as suggested by various studies (e.g. Semenov et al. 2003) The precise evolution of the grain population is a complex issue. Grains evaporate during the first core contraction, before the second collapse. Lenzuni et al. (1995) proposed two main processes of grain destruction: thermal evaporation (destruction directly through thermal vibration), and chemisputtering (reactions between dust and gas). The authors found that carbon evaporates between 750 K and 1100 K, silicates between 1200 K and 1300 K, and aluminium oxides between 1600 K and 1700 K. We here assumed that for each material, the quantity of evaporated grains grows linearly with temperature inside the above ranges until total depletion of this type of grain. Based on the relative fractional abundances of each material calculated from Table 2 in Lenzuni et al. (1995) C 85%, (MgFe)SiO₄ 14.4%, and Al₂O₃ 0.6%, we obtain the three-step evolution curve displayed in Fig. 2.

Thermionic emission.

Thermal agitation on grains induces spontaneous emission of electrons adsorbed on the grain surface (Desch & Turner 2015). Richardson’s law gives the rate of emission $\mathit{\phi }\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{a}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{R}}\frac{\mathrm{4}\mathit{\pi }{\mathit{m}}_{\mathrm{e}}\mathrm{(}{\mathit{k}}_{\mathrm{B}}\mathit{T}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{h}}^{\mathrm{2}}}\exp \left(\frac{\mathit{W}\mathrm{+}\frac{\mathit{Z}{\mathrm{e}}^{\mathrm{2}}}{\mathit{a}}}{{\mathit{k}}_{\mathrm{B}}\mathit{T}}\right)\mathit{,}$(18)with λ_R = 0.5, W = 5eV, Z the grain charge and a its radius.

2.4.3. Grain charges

Grains can hold several electric charges (Draine & Sutin 1987). However, multiply charged grains are weakly abundant in comparison to the singly charged or neutral grains (Sano et al. 2000). Including these grains means handling abundances that range over 20 orders of magnitude, which is numerically difficult to achieve accurately. Furthermore, current analytical models do not correctly describe the charge distribution when grain-grain reactions dominate (see Appendix B for more details). In the present paper, we considered grains holding only one electric charge but we rook the grain-grain reactions into account. We acknowledge that multiply charged grains may change our results, and we will address this issue in future work.

3. Chemistry

3.1. Chemical network

We considered the following elements and their ionised counterparts: H, He, C, O and heavy metallic elements such as Na, Mg, Al, Ca, Fe, Ni, and Si. In conditions typical of molecular clouds and cold neutral medium, H, C and O are primarily found in their molecular forms (H₂, CO, O₂, H₂O, OH). We assumed this to be still the case after the second collapse, but we kept in mind that we lack a precise description of the evolution of these molecules at T> 2000 K. The charged particles taken into account are electrons e⁻, H⁺, He⁺, C⁺, H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$, molecular ions m⁺, and metallic ions M⁺. We considered grains of various sizes (see Sect. 2.4.1), either neutral or with an electric charge ± e. Potassium, sodium, and hydrogen are major contributors to ionised species in number density at high temperature, because of their low ionisation energy. These ionisation reactions (described in Pneuman & Mitchell 1965, see Appendix A) become relevant at T> 1500 K and are included in our network. Since Na and K are alkali metal, we assumed that their reactions and associated rate coefficients are the same as for the other metallic ions M⁺.

Let α_ij represent the reactions of ionisation of species j into i$\begin{array}{ccc}\mathit{j}\mathrm{\to }\mathit{i}\mathrm{+}{\mathrm{e}}^{\mathrm{-}}\mathit{,}& & \end{array}$(19)with the ionisation rate of hydrogen molecules ζ. In our context, UV and radionucleides contributions to ionisation rate are negligible compared to cosmic rays that can deeply penetrate dense cores, and we therefore write ζ = ζ_CR. Let β_ijk represent the reactions between j and k to form i$\begin{array}{ccc}\mathit{j}\mathrm{+}\mathit{k}\mathrm{\to }\mathit{i}\mathrm{+}\mathit{,}& & \end{array}$(20)and β_{∗ ij} (where ∗ denotes any other species that might be present) the reactions between i and j to form another species. We note that β_ijk = β_ikj. We also defined γ_ij to represent the reactions between i and j to form another species, such as β_{∗ ij}, but γ_ij specifically characterises the destruction of i and j rather than the creation of a given species $\begin{array}{ccc}\mathit{i}\mathrm{+}\mathit{j}\mathrm{\to }\mathrm{\ast }\mathit{.}& & \end{array}$(21)Here γ_ij = γ_ji. We then solved the complete set of equations for each charged species (written in dimensionless form) $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \\ \end{array}& & \end{array}$(22)where N is the total number of species (both neutrals and charged particles), n_H is the density of neutrals (here the density of hydrogen molecules), and x_i denotes the fractional abundances of various particles, ${\mathit{x}}_{\mathit{i}}\mathrm{=}\frac{{\mathit{n}}_{\mathit{i}}}{{\mathit{n}}_{\mathrm{H}}}$, and $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathrm{d}\mathit{t}\end{array}\mathrm{=}\mathrm{d}\mathit{t\zeta }$.

We considered that neutral abundances are constant, with values taken from Umebayashi & Nakano (1990), and we solved for the eight above-mentioned cations, plus electrons and grains. The reaction rates were taken from the UMIST database (McElroy et al. 2013) for gas species and Kunz & Mouschovias (2009) for the interactions with and between grains. More details on the chemical network (the considered reactions, the initial abundances, etc.) are given in Appendix A.

3.2. Numerical method

Our resolution method is a semi-implicit scheme. It is unconditionally stable and permits either accurate following of the temporal evolution of the solution (using a stringent constraint on the time step) or acceleration of the equilibrium abundances calculations (with a larger time step). In the latter case, we lose precision on the temporal evolution of the network, but the convergence toward equilibrium is unconditional, which is our main interest.

We write $\mathit{F}\mathrm{(}{x}{\mathrm{)}}_{\mathit{i}}\mathrm{=}{\sum }_{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{N}}[{\mathit{\alpha }}_{\mathit{ij}}{\mathit{x}}_{\mathit{j}}\mathrm{+}\frac{{\mathit{n}}_{\mathrm{H}}}{\mathrm{2}\mathit{\zeta }}{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{\beta }}_{\mathit{ijk}}{\mathit{x}}_{\mathit{j}}{\mathit{x}}_{\mathit{k}}\mathrm{-}\frac{{\mathit{n}}_{\mathrm{H}}}{\mathit{\zeta }}{\mathit{\gamma }}_{\mathit{ij}}{\mathit{x}}_{\mathit{j}}{\mathit{x}}_{\mathit{i}}]$, and Taylor-expand the right-hand side of Eq. (22). This yields $\begin{array}{ccc}{\mathit{x}}_{\mathit{i}}^{\mathit{n}\mathrm{+}\mathrm{1}}& \mathrm{=}& {\mathit{x}}_{\mathit{i}}^{\mathit{n}}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ \mathrm{d}\mathit{t}\end{array}{\mathit{F}}_{\mathit{i}}\mathrm{\left(}{{x}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\right)}\\ & \mathrm{=}& {\mathit{x}}_{\mathit{i}}^{\mathit{n}}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ \mathrm{d}\mathit{t}\end{array}\left[{\mathit{F}}_{\mathit{i}}\mathrm{\left(}{{x}}^{\mathit{n}}\mathrm{\right)}\mathrm{+}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{N}}{J}_{\mathit{ij}}^{\mathit{n}}\mathit{\delta }{\mathit{x}}_{\mathit{j}}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}\mathrm{\right|}\mathrm{\left|}\mathit{\delta }{x}\mathrm{\right|}\mathrm{\left|}\mathrm{\right)}\right]\mathit{,}\end{array}$(23)where J is the Jacobian matrix, x = (x₁,x₂,...,x_N), and δx_i denotes the variation in abundance of the species i between the time steps n and n + 1. In matrix form, this reads $\begin{array}{ccc}\mathrm{(}I\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathrm{d}\mathit{t}\end{array}{J}^{\mathit{n}}\mathrm{)}\hspace{0.17em}\mathit{\delta }{x}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ \mathrm{d}\mathit{t}\end{array}\hspace{0.17em}{F}\mathrm{(}{{x}}^{{n}}\mathrm{)}\mathit{,}& & \end{array}$(24)where I is the identity matrix, and F(x) = (F₁(x),...,F_n(x)).

We limited $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathrm{d}\mathit{t}\end{array}$ by constraining the maximum allowed relative variation during one time step $\begin{array}{ccc}\frac{\mathrm{\Delta }{x}}{{x}}\mathit{<}\mathit{ϵ,}& & \end{array}$(25)with the control parameter ϵ = 10^-2. This permits equilibrium to be rapidly reached and with good precision. The matrix on the left-hand side of Eq. (24) is singular, or numerically very close to singular¹. We used the Singular Value Decomposition (SVD), described in Press et al. (2007) to solve system (24).

We imposed F_i(x) = 0 for species whose abundances are more than eight orders of magnitude smaller than the most abundant ones, because their variations have a negligible influence on the chemical evolution, and this avoids unnecessary small time steps. For T > 1700 K, we also set up the fractional abundance of grains to 10^-30 instead of to 0 to avoid numerical problems.

4. Results

4.1. Fiducial case

Our fiducial case includes ten bins of grain sizes and a standard cosmic-ray ionisation rate of 10^-17 s^-1 (as in Umebayashi & Nakano 1990). The evolution of the various abundances of charged species as function of density during the global collapse of a prestellar core is presented in Fig. 3.

Electrons and M⁺ are the dominant charged species at low densities. For n_H > 10⁸ cm^-3, grains take over while the abundances of all other species decrease by several orders of magnitude. Negatively charged grains are at first more abundant until the neutral grain prevalence for n_H > 10⁹ cm^-3, eventually becoming the main charge carriers alongside positively charged grains. At 10¹⁶ cm^-3, the thermionic emission of grains becomes relevant and releases many electrons. Meanwhile, grain evaporation proceeds through the three stages of destruction that are clearly visible until n_H = 10¹⁸ cm^-3. This affects the abundances of other species, especially M⁺, K⁺, Na⁺ and e⁻. Immediately after the complete destruction of the grains, the thermal ionisation of K, Na and H become important, and their ionised counterparts become the dominant species along with M⁺ and e⁻. Eventually, all neutral K, Na and M atoms disappear, leaving H⁺ as the most abundant ionised species. Grain evaporation takes place at the end of the first core phase, and the thermal ionisations essentially take place during the second collapse. The corresponding conductivities and resistivities are plotted in Figs. 4 and 5. Ohmic and ambipolar resistivities are positive, but Hall conductivities and resistivities have negative values at low densities (light blue curves), before becoming positive (dark blue curves). For n_H < 10¹⁵ cm^-3, these figures are qualitatively comparable to Fig. 7 of Kunz & Mouschovias (2010). We recovered the result of Wardle & Ng (1999) concerning the Hall conductivity, becoming comparable to the Pedersen conductivity for 10⁶ <n_H < 10¹¹ cm^-3 for an MRN grain-size distribution, while it is slightly higher in our case. There is a peak in resistivities at n_H ≈ 10¹⁸ cm^-3 that is due to grain evaporation, where all three resistivities have similar values (but are still dominated by the Ohmic diffusion). The peak does not extend over a wide range of densities because the number density of charged species increases again as soon as thermal ionisations begin, which drastically decreases the resistivities. After the peak, resistivity is dominated by the Ohmic and Hall contributions, which remain comparable, until Ohmic resistivity eventually prevails. Figure 5 also clearly highlights the differences, in the evolution of the various resistivities, between the present calculations and the results of previous studies. For the ambipolar resistivity, Duffin & Pudritz (2008) used the simple analytical expression ${\mathit{\eta }}_{\mathrm{AD}}\mathrm{\propto }\frac{{\mathit{B}}^{\mathrm{2}}}{{\mathit{n}}_{\mathrm{H}}\sqrt{{\mathit{n}}_{\mathrm{H}}}}$, so that their resistivity scales as $\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{H}}}$ (because we consider $\mathit{B}\mathrm{\propto }\sqrt{{\mathit{n}}_{\mathrm{H}}}$). Our ambipolar resistivity is close to theirs at low densities, but the two values diverge around n_H = 10¹⁰ cm^-3. This is due to their assumption that ions are perfectly coupled to the magnetic field, and that the ion-neutral collision time is much shorter than the other characteristic physical times of the system. The authors acknowledge that their model fails for dense regions, with n_H > 10¹⁰ cm^-3. The model of Ohmic resistivity from Machida et al. (2007) also matches ours for n_H ≤ 10¹⁵. The chemical network they considered (Nakano et al. 2002), however, does not directly include potassium in this form, neither the present updated reaction rates nor grain evaporation. At low density (n_H ≤ 10^[15−16] cm^-3 for the Ohmic resistivity and n_H ≤ 10¹⁰ cm^-3 for ambipolar resistivity), these two commonly used models are then in agreement with our results, but the additional physics included in our work significantly improves the situation for the conditions inside the first and second cores.

4.2. Magnetic field variations

The magnetic field only influences the resistivities and the conductivities and has no effect on the chemistry. Figure 6 shows the resistivities for the fiducial case and for magnetic fields three times lower and higher. Hall and ambipolar resistivities are slightly shifted, while the Ohmic resistivity is of course not affected because it does not depend on B. Varying the magnetic field has the strongest effect on the first collapse density range. The density at which the Hall resistivity changes sign during the first core contraction is also shifted, and the ratio between the three resistivities is strongly affected, especially in the 10¹⁶−10¹⁸ cm^-3 density range. Because the resistivities were computed on the fly with the magnetic field of the simulation, this variation will influence the gravitational collapse. As the Hall effect is strongly sign dependent, the formation of some structures, like the first Larson core or the protoplanetary disk, may be delayed or made impossible (Tsukamoto et al. 2015; Wurster et al. 2016).

4.3. Effect of grain evaporation

Figure 7 presents the same calculations as for the fiducial case, but without grain evaporation. There are substantial changes for n_H > 10¹⁶ cm^-3, at the onset of grain destruction. The abundances of both electrons and metallic ions increase as a result of the thermionic emission of the grains and the thermal ionisations. They are the dominant species along with negatively charged grains, which are more prone to form because of the large number of free electrons released by K, Na and H atoms. Consequently, the abundances of neutral and positive grains both drop, while the abundances of other species remain very small. Figure 8 shows the effect that removing grain evaporation has on the resistivities. At high densities (> 10²⁰ cm^-3) the resistivities still decrease because the number of ionised particles increases in the medium. In the density range typical of the second collapse and the second Larson core, 10¹⁵ < n_H < 10²² cm^-3, the Hall and the Ohmic resistivites show significant differences compared to the fiducial case. The Hall resistivity is first lower and then higher than the reference case, while the Ohmic resistivity is tremendously higher in this density range. At n_H ≈ 10²² cm^-3, the excessive abundance of electrons in the medium cancels out the difference between the two cases. An accurate description of grain evaporation is therefore mandatory to study the possible transformation of the first core into a disk around the second Larson core (as described in Machida et al. 2006).

4.4. Cosmic-ray ionisation rates

Cosmic-rays (hereafter CR) propagation along field lines is affected by two competing effects: magnetic focusing, which increases the ionisation rate, and magnetic mirroring, which prevents CRs from reaching deep parts of the cloud. Padovani & Galli (2011), Padovani et al. (2013, 2014a) found that mirroring always dominates focusing for a field topology obtained by ideal-MHD simulations of collapsing rotating clouds and that the ionisation rate could vary by up to a factor 50, depending on the mass-to-flux ratio and the magnetic flux tube considered. In addition, CR are partly absorbed by the dense medium, which lowers the CR ionisation rate (Padovani et al. 2009).

To quantify the effect of these uncertainties on the CR ionisation rate in a protostellar environment with a complex magnetic field topology (twisted field lines, misalignement, turbulence, etc.), we computed abundances and resistivities for two CR values in addition to our standard case, namely ζ_CR = 5 × 10^-18 s^-1 and ζ_CR = 1 × 10^-18 s^-1. The results are shown Fig. 9 (we do not show the abundances because they are similar to the fiducial case). The Hall and ambipolar diffusion contributions are alternatingly increase and decrease. In comparison, our fiducial value ζ_CR = 1 × 10^-17 s^-1 yields about an Ohmic dissipation that is about an order of magnitude smaller before the resistivities decrease at n_H ≈ 10¹⁷ cm^-3. Minor deviations are visible in the second collapse density range, but the overwhelming contribution of electrons and hydrogen compensates for most of the variations.

4.5. Grain-size distribution

Most of the grain surface is due to the smallest grains, therefore a large enough number of grain bins is necessary to properly evaluate of the final conductivity. Similar as Kunz & Mouschovias (2009), we found that five size bins seem to be sufficient to reach a relative precision of the order of 1% for every considered species. In this part, our reference case includes ten size bins. The error on the abundances was calculated as $\frac{\mathrm{\left|}\mathrm{\right|}{{x}}_{{\mathit{N}}_{\mathrm{bins}}}\mathrm{-}{{x}}_{\mathrm{ref}}\mathrm{\left|}\mathrm{\right|}}{\mathrm{\left|}\mathrm{\right|}{{x}}_{\mathrm{ref}}\mathrm{\left|}\mathrm{\right|}}$, where |||| is the L₂ norm of the abundance vector x, x_Nbins is the abundance vector for the number of bins N_bins considered, and x_ref is the abundance vector for the ten bin case. The relative error on the least abundant gas molecules becomes fairly high when grains start to evaporate at 10¹⁷cm^-3, but this is inconsequential since these species hardly contribute to the resistivity at this stage. Five bin abundances clearly yield a smaller error than the ten bins case, which shows that our fiducial calculations have a good precision. The resistivities for different numbers of bins (one, five, and ten) are plotted in Fig. 11. When only one bin is considered, the resistivities are shifted up and down by about one order of magnitude and are strongly overestimated in the 10¹²−10¹⁷cm^-3 density range. However, the five and ten bins cases are extremely close to each other. A large enough number of bins is thus necessary to properly describe this highly dynamic phase of the prestellar core evolution.

Even though the MRN grain-size distribution is reasonable when considering interstellar dust, it becomes of questionable validity for denser media such as molecular clouds or Larson cores, for which a precise grain-size distribution is still lacking. We examined the effect of this uncertainty by conducting calculations with another distribution: λ = −2.8 (instead of λ = −3.5 for the MRN distribution), as suggested by Compiègne et al. (2011) for amorphous carbon grains between radii of~ 5 and few hundred nanometers. The result for the resistivities is shown in Fig. 12. Despite the similarity of the comportment of the resistivities between 10⁵ cm^-3 and 10¹⁷ cm^-3, there is a slight difference in the partitioning for the dominant effect. For example, around n_H = 10¹⁴ cm^-3, Ohmic diffusion takes the lead upon the ambipolar diffusion for a density ten times higher than the MRN distribution. This density range is typical of protostellar disks, which means that the respective evolution of these latter might be quite different because the size distributions may vary from one to another, with effects difficult to control (grain coalescence, premature destruction, etc.). Several authors proposed their own size distribution. Compiègne et al. (2011) also assumed a lognormal distribution for the smallest grains (PAH and small amorphous carbon) and a power law with an exponential cut-off for larger amorphous carbon grains and amorphous silicate grains to reproduce observed emission and extinction spectra (see their Fig. 2).

4.6. Non-equilibrium chemistry

Until now, all the species abundances were calculated at chemical equilibrium. In real situations, the chemical reaction timescale could be greater than the dynamical (collapse) one. In that case, the environment conditions (density, temperature, etc.) can change significantly before chemical equilibrium is reached. We explored this possibility by comparing the timescale required for our chemical network to reach equilibrium with the typical free-fall time for a self-gravitating spherical cloud, ${\mathit{t}}_{\mathrm{ff}}\mathrm{=}\sqrt{\frac{\mathrm{3}\mathit{\pi }}{\mathrm{32}\mathit{G}\mathit{\rho ̅}}}$ (with the mean density of the cloud). The results are portrayed in Fig. 13. Clearly, the free-fall time is far longer than the chemical equilibrium time at all densities. Furthermore, the above free-fall time estimate is likely to be underestimated since real clouds are additionally supported by thermal, turbulent and magnetic pressures. Therefore, assuming chemical equilibrium seems to be valid in the context of collapsing prestellar cores, with the limitation that we did not consider the history of the chemical species during the collapse. Particles may be transported from one specific environment to another (e.g. from the core region to the outflow then to the disk). The only way to properly account for this complicated pattern is to calculate non-equilibrium chemical reactions on the fly during the dynamical collapse. Additionally the flow of the fluid changes the dynamics and statistics of the collisions, which results in chemical transformations and should therefore also be taken into account. This task remains computationally prohibitive for now.

5. Conclusion

We have developed a solver that computes a detailed network of the main chemical reactions relevant to the first and second collapse of prestellar cores. The network is based on the work of Umebayashi & Nakano (1990) but extends this study significantly by updating the reaction rates and including the effects of dust evaporation, thermal ionisation of potassium, and exploring various cosmic-ray ionisation rates. We also used a distribution of grain sizes, the MRN distribution, and explored the effect of the number of size-bins on the results. The solver yields the chemical equilibrium abundances of the various neutral and ionised species². The various abundances were used to calculate the non-ideal MHD resistivities, namely Ohmic, Hall and ambipolar resistivities, during the collapse, using a barotropic EOS to reproduce the typical density-temperature conditions. The resistivities determine the dynamics of the first and second collapse, and thus the properties of the first and second Larson cores.

Above a temperature T ~ 500 K , the effects of thermionic emission, grain evaporation and thermal ionisation become preponderant. An accurate description of these processes is mandatory to properly characterise the collapse because they occur during the first core contraction and influence the initial conditions of the second collapse. Dust destruction has a double effect on the collapse. First, it modifies the opacity of the medium (see e.g. Lenzuni et al. 1995), which regulates the radiative cooling of the system. Evaporation increases the efficiency of this process, and the collapse accelerates as the gas is less thermally supported. Second, it affects the various resistivities of the non-ideal MHD terms with direct effect on magnetic field topology and a diminution of the magnetic braking.

We did not include the molecular dissociation of elements in the chemical network, in particular the dissociation of H₂ at 2000 K, which leads to the second collapse. However, we do not expect the resistivities to be significantly affected by this process because H⁺ ions and electrons are the dominant charged species at these temperatures. As mentioned in the text, our grain model is simplified because we assumed that each grains was composed of only one material. A more realistic structure containing several layers of different materials that evaporate one after the other will be included in future work.

In addition to the general outcome of these calculations and their effect on prestellar core evolution, we wish to highlight the following points

• At least five bins for the grain-size distribution are necessary toreliably determine of the resistivities. As discussed inSect. 4.5, a precise knowledge of the grain-sizedistribution would certainly improve the reliability of the results.

• As shown in Sect. 4.4, resistivities change by up to two orders of magnitude when ionisation rates vary by a factor 10. This highlights the importance of shielding against cosmic rays in collapsing cores (Padovani et al. 2014b).

• In the temperature range 750 K ≲ T ≲ 1700 K, dust grains evaporate. This evaporation has tremendous consequences on the various chemical abundances and thus on the resistivities, since grains are the main contributors to the resistivities at these temperatures.

• Around 1500 K and above, at which thermal ionisation of metallic ions and hydrogen occur, grains have been entirely destroyed and H⁺ and electrons become the main charge carriers, which causes the resistivities to drop even further.

• Our chemical integration time is always shorter than the free-fall time. We can therefore assume equilibrium chemistry, which is less demanding than non-equilibrium chemistry, especially during hydrodynamics simulations.

This solver allows us to compute a large multi-dimensional multi-species equilibrium abundance table for a wide range in temperatures, densities and ionisation rates. This table can be used during simulations of the first and second collapse of prestellar cores, allowing a consistent dynamical-chemical description of this process. The table can be downloaded at https://bitbucket.org/pmarchan/chemistry and at the CDS.

---

Acknowledgments

The research leading to these results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007−2013 Grant Agreement No. 247060). B.C. gratefully acknowledges support from the French ANR Retour Postdoc program (ANR-11-PDOC-0031). N.V. is supported by the European Commission through the Horizon 2020 Marie Skłodowska-Curie Actions Individual Fellowship 2014 programme (Grant Agreement No. 659706). We acknowledge financial support from the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU, France. We thank T. Grassi for useful discussions during the writing of this paper. Finally, we would also like to thank the referee for very insightful comments that have greatly helped to improve the completeness of this work.

References

Appendix A: Details of the chemical network

Table A.1 gives a list of all the reactions implemented in the code, taken from Umebayashi & Nakano (1990). The reaction rates are taken from the UMIST Database (McElroy et al. 2013). The , , and coefficients are defined as $\begin{array}{ccc}\mathit{k}\mathrm{=}\mathit{\alpha ̅}{\left(\frac{\mathit{T}}{\mathrm{300}}\right)}^{\mathit{\beta ̅}}{\mathrm{e}}^{\mathrm{-}\frac{\mathit{\gamma ̅}}{\mathit{T}}}\mathit{,a}\mathrm{=}\mathrm{1}& & \end{array}$with k the reaction rate. The units of k are cm³s^-1 for reactions between species and s^-1 for ionisation and photo-reactions with cosmic rays. M stands for metals such as Mg, Al, Ca and Fe. Following Umebayashi & Nakano (1990), we used for all these elements the typical coefficient rates of Mg, which is the most abundant species. m stands for molecules that can be ionised and are represented by HCO⁺.

We mention that the two reactions H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ + e⁻→ H + H + H and H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ + e⁻→ H₂ + H occur at the same rate. The ionisation rate for potassium, sodium and hydrogen are given by Pneuman & Mitchell (1965)$\begin{array}{ccc}\frac{\mathrm{d}{\mathit{n}}_{{\mathrm{K}}^{\mathrm{+}}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{6.5}\mathrm{\times }{\mathrm{10}}^{-15}{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}{\mathit{T}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\times }{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{5.1}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{\mathrm{4}}\hspace{0.17em}\mathrm{K}}{\mathit{T}}}{\mathrm{cm}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}& & \\ \frac{\mathrm{d}{\mathit{n}}_{\mathrm{N}{\mathrm{a}}^{\mathrm{+}}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{1.4}\mathrm{\times }{\mathrm{10}}^{-15}{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}{\mathit{T}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\times }{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{6.0}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{\mathrm{4}}\hspace{0.17em}\mathrm{K}}{\mathit{T}}}{\mathrm{cm}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}& & \\ \frac{\mathrm{d}{\mathit{n}}_{{\mathrm{H}}^{\mathrm{+}}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{2.0}\mathrm{\times }{\mathrm{10}}^{-10}{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}{\mathit{T}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\times }{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{15.8}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{\mathrm{4}}\hspace{0.17em}\mathrm{K}}{\mathit{T}}}{\mathrm{cm}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{.}& & \end{array}$We considered recombination reactions at the surface of the grains. We used the collision rates of Draine & Sutin (1987) and the interactions described in Umebayashi & Nakano (1990), Ilgner & Nelson (2006), and Kunz & Mouschovias (2009):

• Neutral grains: when hit by an electron, the elec-tron sticks onto the grain with a probability ofP_e = 60%, while an ion always sticks.

• Negatively charged grains: when hit by an ion, the ion recombines.

• Positively charged grains: when hit by an electron, the grain becomes neutral.

• Charged grains: if two grains with opposite charges collide, one electric charge is transferred.

• Neutral grains: if hit by a charged grain, singly charged is transferred to the neutral grain.

The mean collision rate ⟨σv⟩ between a grain with a radius a and a charge q_g and another species with a mass m and a charge q_s is $⟨\mathit{\sigma v}⟩\mathrm{=}\mathit{\pi }{\mathit{a}}^{\mathrm{2}}{\left(\frac{\mathrm{8}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\pi m}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{q}}_{\mathrm{s}}{\mathit{q}}_{\mathrm{g}}}{\mathit{a}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\right)\left(\mathrm{1}\mathrm{+}{\left(\frac{\mathrm{2}}{\frac{\mathit{a}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{{\mathrm{e}}^{\mathrm{2}}}\mathrm{-}\mathrm{2}\frac{{\mathit{q}}_{\mathit{g}}}{{\mathit{q}}_{\mathit{s}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)$(A.4)for q_sq_g< 0, $\begin{array}{ccc}⟨\mathit{\sigma v}⟩\mathrm{=}\mathit{\pi }{\mathit{a}}^{\mathrm{2}}{\left(\frac{\mathrm{8}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\pi m}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathrm{1}\mathrm{+}{\left(\frac{\mathrm{4}\mathit{a}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{{\mathrm{e}}^{\mathrm{2}}}\mathrm{+}\mathrm{3}\frac{{\mathit{q}}_{\mathrm{g}}}{{\mathit{q}}_{\mathrm{s}}}\right)}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{2}}& & \\ \mathrm{\times }\exp \left(\mathrm{-}\frac{{\mathit{q}}_{\mathrm{g}}{\mathrm{e}}^{\mathrm{2}}}{\mathrm{\left(}{\mathit{q}}_{\mathrm{s}}\mathit{a}{\mathit{k}}_{\mathrm{b}}\mathit{T}\mathrm{\right)}\left(\mathrm{1}\mathrm{+}{\left(\frac{{\mathit{q}}_{\mathrm{s}}}{{\mathit{q}}_{\mathrm{g}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)}\right)& & \end{array}$(A.5)for q_sq_g> 0 and $⟨\mathit{\sigma v}⟩\mathrm{=}\mathit{\pi }{\mathit{a}}^{\mathrm{2}}{\left(\frac{\mathrm{8}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\pi m}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}\mathrm{+}{\left(\frac{\mathit{\pi }{\mathrm{e}}^{\mathrm{2}}}{\mathrm{2}\mathit{a}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\right)\mathrm{·}$(A.6)for q_g = 0.

For two grains of opposite charges q and q′ with radii a and a′ and a reduced mass ${\mathit{\mu }}_{\mathrm{g}}\mathrm{=}\frac{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}{\mathit{m}\mathrm{+}{\mathit{m}}^{\mathrm{\prime }}}$, the collision rate is $⟨\mathit{\sigma v}⟩\mathrm{=}\mathit{\pi }\mathrm{(}\mathit{a}\mathrm{+}{\mathit{a}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{2}}{\left(\frac{\mathrm{8}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\pi }{\mathit{\mu }}_{\mathrm{g}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\frac{\mathit{q}{\mathit{q}}^{\mathrm{\prime }}}{\mathrm{(}\mathit{a}\mathrm{+}{\mathit{a}}^{\mathrm{\prime }}\mathrm{)}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\right)\mathit{,}$(A.7)and $⟨\mathit{\sigma v}⟩\mathrm{=}\mathit{\pi }\mathrm{(}\mathit{a}\mathrm{+}{\mathit{a}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{2}}{\left(\frac{\mathrm{8}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\pi }{\mathit{\mu }}_{\mathrm{g}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}\mathrm{+}{\left(\frac{\mathit{\pi }{\mathrm{e}}^{\mathrm{2}}}{\mathrm{2}\mathrm{(}\mathit{a}\mathrm{+}{\mathit{a}}^{\mathrm{\prime }}\mathrm{)}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\right){\mathit{P}}_{\mathit{a}}\mathit{,}$(A.8)between a charged and a neutral grain. ${\mathit{P}}_{\mathit{a}}\mathrm{=}\frac{{\mathit{a}}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{{}^{\mathrm{\prime }\mathrm{2}}}}$ denotes the probability of a charge transfer to the neutral grain of radius a′.

The abundances of the neutrals are also taken from Umebayashi & Nakano (1990; Kunz & Mouschovias 2009, for potassium). Because of its high thermal ionisation rate, K is the only neutral species whose fractional abundance varies in our network.

Appendix B: Grain charges

Multiply charged grains are part of a chemical network (Nishi et al. 1991). Figure B.1 shows the abundances of species in a simplified case (only one bin of size, without Na, the thermionic emission and thermal ionisation of H) but with grains holding two electric charges. Although they are generally less abundant than single-charge grains, they effect the charge distribution, especially for n_H> 10¹² cm^-3. The main effect is that one-charge grains, the main charge carriers, seem to form more neutral grains and are less numerous than in the fiducial case. This of course affects the resistivities, as shown in Fig. B.2.

We were only able to produce such a result for one bin of size with the current method of resolution. The abundance gap between the two-charge holding grains and the remaining the species spans over 20 orders of magnitude (especially at low densities in this case). This gap widens for several number of bins because small grains tend to hold fewer charges, and large grains are less abundant. Dealing with so many orders of magnitude is problematic because numeric round-off errors may artificially change the contribution of the less abundant species. Several authors have used analytical expressions of charge distribution to avoid this numerical difficulty (Draine & Sutin 1987; Okuzumi 2009; Fujii et al. 2011). These models, however, do not take the charge transfer between grains into account. Figure B.3 shows the highest reaction rates of grain-ion, grain-electron and grain-grain reactions. Reaction rates between grains are largely dominant for n_H> 10¹² cm^-3, and not taking them into account leads to large discrepancies in the grains abundances in this density range (see Fig. B.4), with in turn affects the resistivities.

On the other hand, using the analytical expressions of Draine & Sutin (1987) with the ion and electron abundances calculated in our simplified case with only one charge and without considering grain-grain reactions enables us to predict the grain abundances for one and two charges with good precision. Figure B.5 depicts these predicted abundances of grains, compared to those computed with our code, for one and two charges. The abundances of the doubly charged grains are slightly overestimated for 10¹⁴<n_H< 10¹⁷ cm^-3, but the agreement is very good otherwise. Figure B.6 shows the relative error $\frac{{{x}}_{\mathrm{predict}}\mathrm{-}{{x}}_{\mathrm{num}}}{{{x}}_{\mathrm{num}}}$ of the model prediction for the numerical run with two charges allowed. All abundances agree within roughly 10%, except for g⁺⁺ at low densities and g⁻ and g⁺⁺ for n_H> 10¹² cm^-3. Figure B.7 shows the corresponding resistivities. Both models are relatively similar except in the density ranges mentioned above. For 10¹⁴<n_H< 10¹⁷ cm^-3, the model resistivities are slightly underestimated, but the differences with the numerical results remain similar. At low densities, the errors on the abundances result in a sign change of the Hall resistivity because even a small change allows the positive contribution of η_H to overcome the negative contribution. However, the Hall effect is not expected to play a role at these early stages of protostellar collapse. Therefore, even though we were unable to verify the exactness of the prediction for grains holding a larger number of charges (because of numerical difficulties), it seems reasonable to say that our method gives satisfying results at least for a small number of charges, which is the relevant case because grains holding more than three charges are not expected to be abundant enough to significantly modify the resistivities (Sano et al. 2000).

Appendix C: Details of the numerical method

The Jacobian general expression (for charged particles alone, neutral grains are taken care of separately) reads: $\begin{array}{ccc}{J}_{\mathit{ij}}& \mathrm{=}& \frac{\mathit{\partial }{\mathit{F}}_{\mathit{i}}}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}\\ & \mathrm{=}& {\mathit{\alpha }}_{\mathit{ij}}\mathrm{+}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{\beta }}_{\mathit{ijk}}{\mathit{x}}_{\mathit{k}}\mathrm{-}{\mathit{\gamma }}_{\mathit{ij}}{\mathit{x}}_{\mathit{i}}\mathrm{-}{\mathit{\delta }}_{\mathit{ij}}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{\gamma }}_{\mathit{ik}}{\mathit{x}}_{\mathit{k}}\mathit{.}\end{array}$(C.1)The Jacobian expression for neutral grains in the bin α is: $\begin{array}{ccc}{J}_{{\mathit{g}}_{\mathrm{0}}^{\mathit{\alpha }}\mathit{j}}\mathrm{=}\mathrm{-}{J}_{{\mathit{g}}_{\mathrm{+}}^{\mathit{\alpha }}\mathit{j}}\mathrm{-}{J}_{{\mathit{g}}_{\mathrm{-}}^{\mathit{\alpha }}\mathit{j}}\mathit{.}& & \end{array}$(C.2)This makes the Jacobian matrix J singular, which is a great problem when solving the system using Newton-Raphson iterations. For the semi-implicit method, however, a simple solution consists of reducing the calculated species to charged particles alone and then updating with the neutral grains at the end of the iteration.

A schematic diagram of the system is shown in Fig. C.1, with $\begin{array}{ccc}{\mathit{A}}_{\mathit{ij}}& \mathrm{=}& {\mathit{\delta }}_{\mathit{ij}}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathrm{d}\mathit{t}\end{array}\left[{\mathit{\alpha }}_{\mathit{ij}}\mathrm{-}\frac{{\mathit{n}}_{\mathrm{H}}}{\mathit{\zeta }}{\mathit{\gamma }}_{\mathit{ij}}{\mathit{x}}_{\mathit{i}}^{\mathit{n}}\mathrm{+}\frac{{\mathit{n}}_{\mathrm{H}}}{\mathit{\zeta }}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{(}{\mathit{\beta }}_{\mathit{ijk}}\mathrm{-}{\mathit{\gamma }}_{\mathit{ik}}{\mathit{\delta }}_{\mathit{ij}}\mathrm{)}{\mathit{x}}_{\mathit{k}}^{\mathit{n}}\right]\mathit{,}\\ {\mathit{B}}_{\mathit{i}}& \mathrm{=}& \frac{{\mathit{n}}_{\mathrm{H}}}{\mathit{\zeta }}\left[\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{\gamma }}_{\mathit{ij}}{\mathit{x}}_{\mathit{j}}^{\mathit{n}}{\mathit{x}}_{\mathit{i}}^{\mathit{n}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{N}}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{\beta }}_{\mathit{ijk}}{\mathit{x}}_{\mathit{j}}^{\mathit{n}}{\mathit{x}}_{\mathit{k}}^{\mathit{n}}\right]\begin{array}{c}\mathrm{˜}\\ \mathrm{d}\mathit{t}\end{array}\mathrm{+}{\mathit{x}}_{\mathit{i}}^{\mathit{n}}\mathit{.}\end{array}$While the system nearly reaches equilibrium, the abundances of the least abundant species may continue to vary significantly in relative values for a long time. For this reason, we chose to stop the calculations when the time-step $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathrm{d}\mathit{t}\end{array}$ reached a final constant value. The least abundant species may continue to evolve, but this is inconsequential for evaluating the equilibrium abundances of the dominant species.

Appendix D: Code validation

In order to test the code with the most simple parameter set, we compared it with the results of Umebayashi & Nakano (1990), which did not include potassium and grain evaporation. We considered the basic case of one bin of grains, with a size a₀. Let δ₁ be the fraction of C and O in the gas phase, and δ₂ the fraction of metals. Our calculations are shown in Fig. D.1, with δ₁ = 0.2 and δ₂ = 0.02 for the top panel, and δ₁ = δ₂ = 0 for the bottom panel. It represents the evolution of the relative abundances with the density at T = 10 K, compared with data from Umebayashi & Nakano (1990). It shows that the two agree very well. For n_H< 10¹⁰ cm^-3, most of the grains have a negative charge and electrons and ions are the dominant species. For n_H> 10¹³ cm^-3, charged grains are more abundant than electrons and ions.

All Tables

All Figures