EDP Sciences
Free Access
Issue
A&A
Volume 592, August 2016
Article Number A18
Number of page(s) 14
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201526780
Published online 04 July 2016

© 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 H2 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 H2 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 ~104 astronomical units (AU) with an average density of a few 103 cm-3, while the protostellar core measures only 10-3 AU at birth, with densities over 1020 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 (1)where B denotes the magnetic field, | | B | | its L2 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 where the parallel, perpendicular, and Hall conductivites are in turn expressed, respectively, as with , the cyclotron frequency, and (8)Here i stands for any charged particle of charge qi, with a charge particle number-density ni, m denotes the mass of a particle, and c is the speed of light. The factor aiHe 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). σcollwi is the rate constant for collisions between a particle i and H2 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 qi. 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 102 cm-3<n< 1026 cm-3 and 10 K <T< 105 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)(9)with n the total density and (10)and T0 = 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).

thumbnail Fig. 1

Solid line: temperature as a function of density for the equation of state. Dashed line: magnetic field prescription following Li et al. (2011).

Open with DEXTER

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 (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 Nbins of size bins with equal widths in log space. We define a minimum and maximum grain size amin and amax, respectively, with a number density of grains of radius between a and a + da(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 amin = 0.0181 μm and amax = 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,...,Nbins) where . The total number density of dust ng,tot is determined by constraining the total grain mass density in the size distribution to be (14)which yields (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 a0 (see Kunz & Mouschovias 2009). Because (16)this yields (17)where, for a dust to gas ratio of 1%, (ρ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 .

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)SiO4) and aluminium oxyde (Al2O3). 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)SiO4 14.4%, and Al2O3 0.6%, we obtain the three-step evolution curve displayed in Fig. 2.

thumbnail Fig. 2

Abundance of grains (relative to their value at T< 750 K) as a function of temperature. The most abundant species (C, MgFeSiO4 and Al2O3) evaporate during the three represented main destruction stages.

Open with DEXTER

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 (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 (H2, CO, O2, H2O, 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, 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+.

thumbnail Fig. 3

Evolution of the fractional abundances for the charged species with density for our fiducial case: 10 bins of grains, barotropic equation of state and a cosmic-ray ionisation rate ζCR = 10-17 s-1.

Open with DEXTER

Let αij represent the reactions of ionisation of species j into i(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(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 (21)Here γij = γji. We then solved the complete set of equations for each charged species (written in dimensionless form) (22)where N is the total number of species (both neutrals and charged particles), nH is the density of neutrals (here the density of hydrogen molecules), and xi denotes the fractional abundances of various particles, , and .

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 , and Taylor-expand the right-hand side of Eq. (22). This yields (23)where J is the Jacobian matrix, x = (x1,x2,...,xN), and δxi denotes the variation in abundance of the species i between the time steps n and n + 1. In matrix form, this reads (24)where I is the identity matrix, and F(x) = (F1(x),...,Fn(x)).

We limited by constraining the maximum allowed relative variation during one time step (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 singular1. We used the Singular Value Decomposition (SVD), described in Press et al. (2007) to solve system (24).

We imposed Fi(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.

thumbnail Fig. 4

Conductivities as a function of density for the fiducial case. Light blue: negative Hall conductivity, dark blue: positive Hall conductivity, red: parallel conductivity, green: perpendicular conductivity.

Open with DEXTER

Electrons and M+ are the dominant charged species at low densities. For nH > 108 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 nH > 109 cm-3, eventually becoming the main charge carriers alongside positively charged grains. At 1016 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 nH = 1018 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 nH < 1015 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 106 <nH < 1011 cm-3 for an MRN grain-size distribution, while it is slightly higher in our case. There is a peak in resistivities at nH ≈ 1018 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 , so that their resistivity scales as (because we consider ). Our ambipolar resistivity is close to theirs at low densities, but the two values diverge around nH = 1010 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 nH > 1010 cm-3. The model of Ohmic resistivity from Machida et al. (2007) also matches ours for nH ≤ 1015. 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 (nH ≤ 10[15−16] cm-3 for the Ohmic resistivity and nH ≤ 1010 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.

thumbnail Fig. 5

Resistivities as a function of density for the fiducial case. Light blue: negative Hall resistivity, dark blue: positive Hall resistivity, red: ambipolar diffusion resistivity, green: Ohm resistivity.

Open with DEXTER

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 1016−1018 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).

thumbnail Fig. 6

Resistivities as a function of density for the fiducial case, and a magnetic field three times higher and lower. Same colour coding as Fig. 5.

Open with DEXTER

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 nH > 1016 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 (> 1020 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, 1015 < nH < 1022 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 nH ≈ 1022 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).

thumbnail Fig. 7

Fractional abundances, using ten bins, without grain evaporation.

Open with DEXTER

thumbnail Fig. 8

Comparison of the resistivities with and without grain evaporation. Same colour coding as in Fig. 5. Dashed lines represent the absence of evaporation, whereas solid lines stand for the fiducial case.

Open with DEXTER

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 nH ≈ 1017 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.

thumbnail Fig. 9

Resistivities for different cosmic-ray ionisation rates, with the barotropic EOS. Solid lines: fiducial value ζCRref = 1 × 10-17 s-1. Dashed lines: . Dotted lines: . Same colour coding as Fig. 5.

Open with DEXTER

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 , where |||| is the L2 norm of the abundance vector x, xNbins is the abundance vector for the number of bins Nbins considered, and xref 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 1017cm-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 1012−1017cm-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.

thumbnail Fig. 10

Relative error in the species abundances () for different numbers of size bins (relatively to the ten bins results) for grains only (top panel), and for other species (bottom panel) in the density range before grain evaporation.

Open with DEXTER

thumbnail Fig. 11

Resistivities for different numbers of size bins. Solid line: ten bin fiducial case, dotted line: five bins, dash-dotted line: one bin. Same colour coding as Fig. 5.

Open with DEXTER

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 105 cm-3 and 1017 cm-3, there is a slight difference in the partitioning for the dominant effect. For example, around nH = 1014 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).

thumbnail Fig. 12

Ambipolar diffusion resistivity for the MRN λ = −3.5 (full line) and the Compiègne et al. (2011)λ = −2.8 (dashed line) distributions, with the same (our) barotropic EOS.

Open with DEXTER

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, (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.

thumbnail Fig. 13

Typical chemical equilibrium timescale (solid red) and dynamical timescale (dash green) as a function of density.

Open with DEXTER

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 species2. 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 H2 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.


1

Because of the neutral grains, although it is very close to singular even without them, because of very rare species.

2

It can also be used to derive out-of-equilibrium abundances of given species during the collapse, if necessary.

Acknowledgments

The research leading to these results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/20072013 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 with k the reaction rate. The units of k are cm3s-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+.

Table A.1

Chemical reactions and coefficient rates of the chemical network.

We mention that the two reactions H + e H + H + H and H + e H2 + H occur at the same rate. The ionisation rate for potassium, sodium and hydrogen are given by Pneuman & Mitchell (1965)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 ofPe = 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 qg and another species with a mass m and a charge qs is (A.4)for qsqg< 0, (A.5)for qsqg> 0 and (A.6)for qg = 0.

Table A.2

Abundances of neutrals relative to H.

For two grains of opposite charges q and q with radii a and a and a reduced mass , the collision rate is (A.7)and (A.8)between a charged and a neutral grain. 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 nH> 1012 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.

thumbnail Fig. B.1

Abundances of the chemical species with double-charge grains, without Na, the thermal ionisation of H and the thermionic emission.

Open with DEXTER

thumbnail Fig. B.2

Resistivity comparison for the simplified case. Dashed line: network with only singly charged grains, solid line: network with singly and doubly charged grains.

Open with DEXTER

thumbnail Fig. B.3

Highest reactions rates for grain/ion and grain/electron reactions (red curve), positively charged grain/negatively charged grain reactions (green curve) and neutral grain/charged grain reactions (blue curve).

Open with DEXTER

thumbnail Fig. B.4

Comparison of the grain abundances with (dashed lines) and without (solid lines) charge transfer between grains for singly charged grains (top panel) and singly and doubly charged grains (bottom panel).

Open with DEXTER

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 nH> 1012 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.

thumbnail Fig. B.5

Comparison of the abundances of grains computed with our code (dashed lines) and predicted from the ions and electron abundances (solid lines) for singly charged grains (top panel) and singly and doubly charged grains with one bin (bottom panel).

Open with DEXTER

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 1014<nH< 1017 cm-3, but the agreement is very good otherwise. Figure B.6 shows the relative error 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 nH> 1012 cm-3. Figure B.7 shows the corresponding resistivities. Both models are relatively similar except in the density ranges mentioned above. For 1014<nH< 1017 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).

thumbnail Fig. B.6

Relative error between the numerical abundances of grains and the predicted results of Fig. B.5 for singly and doubly charged grains.

Open with DEXTER

thumbnail Fig. B.7

Comparison of the resistivities using abundances computed with our code (dashed lines) and abundances predicted with the Draine & Sutin (1987) formulae (solid lines) for singly and doubly charged grains, without charge transfer between grains.

Open with DEXTER

Appendix C: Details of the numerical method

thumbnail Fig. C.1

Matrix visualisation of the system to be solved after linearisation in time.

Open with DEXTER

The Jacobian general expression (for charged particles alone, neutral grains are taken care of separately) reads: (C.1)The Jacobian expression for neutral grains in the bin α is: (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 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 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 a0. Let δ1 be the fraction of C and O in the gas phase, and δ2 the fraction of metals. Our calculations are shown in Fig. D.1, with δ1 = 0.2 and δ2 = 0.02 for the top panel, and δ1 = δ2 = 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 nH< 1010 cm-3, most of the grains have a negative charge and electrons and ions are the dominant species. For nH> 1013 cm-3, charged grains are more abundant than electrons and ions.

thumbnail Fig. D.1

Evolution of the abundances of various charged particles as a function of the number density of hydrogen atoms, at T = 10 K. Solid lines are our results, and circles are taken from Figs. 2 and 4 of Umebayashi & Nakano (1990). Top: δ1 = 0.2, δ2 = 0.02. Bottom: δ1 = δ2 = 0.

Open with DEXTER

All Tables

Table A.1

Chemical reactions and coefficient rates of the chemical network.

Table A.2

Abundances of neutrals relative to H.

All Figures

thumbnail Fig. 1

Solid line: temperature as a function of density for the equation of state. Dashed line: magnetic field prescription following Li et al. (2011).

Open with DEXTER
In the text
thumbnail Fig. 2

Abundance of grains (relative to their value at T< 750 K) as a function of temperature. The most abundant species (C, MgFeSiO4 and Al2O3) evaporate during the three represented main destruction stages.

Open with DEXTER
In the text
thumbnail Fig. 3

Evolution of the fractional abundances for the charged species with density for our fiducial case: 10 bins of grains, barotropic equation of state and a cosmic-ray ionisation rate ζCR = 10-17 s-1.

Open with DEXTER
In the text
thumbnail Fig. 4

Conductivities as a function of density for the fiducial case. Light blue: negative Hall conductivity, dark blue: positive Hall conductivity, red: parallel conductivity, green: perpendicular conductivity.

Open with DEXTER
In the text
thumbnail Fig. 5

Resistivities as a function of density for the fiducial case. Light blue: negative Hall resistivity, dark blue: positive Hall resistivity, red: ambipolar diffusion resistivity, green: Ohm resistivity.

Open with DEXTER
In the text
thumbnail Fig. 6

Resistivities as a function of density for the fiducial case, and a magnetic field three times higher and lower. Same colour coding as Fig. 5.

Open with DEXTER
In the text
thumbnail Fig. 7

Fractional abundances, using ten bins, without grain evaporation.

Open with DEXTER
In the text
thumbnail Fig. 8

Comparison of the resistivities with and without grain evaporation. Same colour coding as in Fig. 5. Dashed lines represent the absence of evaporation, whereas solid lines stand for the fiducial case.

Open with DEXTER
In the text
thumbnail Fig. 9

Resistivities for different cosmic-ray ionisation rates, with the barotropic EOS. Solid lines: fiducial value ζCRref = 1 × 10-17 s-1. Dashed lines: . Dotted lines: . Same colour coding as Fig. 5.

Open with DEXTER
In the text
thumbnail Fig. 10

Relative error in the species abundances () for different numbers of size bins (relatively to the ten bins results) for grains only (top panel), and for other species (bottom panel) in the density range before grain evaporation.

Open with DEXTER
In the text
thumbnail Fig. 11

Resistivities for different numbers of size bins. Solid line: ten bin fiducial case, dotted line: five bins, dash-dotted line: one bin. Same colour coding as Fig. 5.

Open with DEXTER
In the text
thumbnail Fig. 12

Ambipolar diffusion resistivity for the MRN λ = −3.5 (full line) and the Compiègne et al. (2011)λ = −2.8 (dashed line) distributions, with the same (our) barotropic EOS.

Open with DEXTER
In the text
thumbnail Fig. 13

Typical chemical equilibrium timescale (solid red) and dynamical timescale (dash green) as a function of density.

Open with DEXTER
In the text
thumbnail Fig. B.1

Abundances of the chemical species with double-charge grains, without Na, the thermal ionisation of H and the thermionic emission.

Open with DEXTER
In the text
thumbnail Fig. B.2

Resistivity comparison for the simplified case. Dashed line: network with only singly charged grains, solid line: network with singly and doubly charged grains.

Open with DEXTER
In the text
thumbnail Fig. B.3

Highest reactions rates for grain/ion and grain/electron reactions (red curve), positively charged grain/negatively charged grain reactions (green curve) and neutral grain/charged grain reactions (blue curve).

Open with DEXTER
In the text
thumbnail Fig. B.4

Comparison of the grain abundances with (dashed lines) and without (solid lines) charge transfer between grains for singly charged grains (top panel) and singly and doubly charged grains (bottom panel).

Open with DEXTER
In the text
thumbnail Fig. B.5

Comparison of the abundances of grains computed with our code (dashed lines) and predicted from the ions and electron abundances (solid lines) for singly charged grains (top panel) and singly and doubly charged grains with one bin (bottom panel).

Open with DEXTER
In the text
thumbnail Fig. B.6

Relative error between the numerical abundances of grains and the predicted results of Fig. B.5 for singly and doubly charged grains.

Open with DEXTER
In the text
thumbnail Fig. B.7

Comparison of the resistivities using abundances computed with our code (dashed lines) and abundances predicted with the Draine & Sutin (1987) formulae (solid lines) for singly and doubly charged grains, without charge transfer between grains.

Open with DEXTER
In the text
thumbnail Fig. C.1

Matrix visualisation of the system to be solved after linearisation in time.

Open with DEXTER
In the text
thumbnail Fig. D.1

Evolution of the abundances of various charged particles as a function of the number density of hydrogen atoms, at T = 10 K. Solid lines are our results, and circles are taken from Figs. 2 and 4 of Umebayashi & Nakano (1990). Top: δ1 = 0.2, δ2 = 0.02. Bottom: δ1 = δ2 = 0.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.