© ESO, 2014

1. Introduction

The study of the formation of circumstellar discs around protostars still presents considerable theoretical challenges. Because of the presence of magnetic fields in the parent molecular clouds and cores (Crutcher 2012), assuming the strict conservation of angular momentum during cloud collapse and star formation is not warranted. In fact, according to recent numerical and analytical studies, the main effect of a magnetic field entrained by a collapsing cloud is to brake any rotational motion, at least as long as the field remains frozen in the gas and the rotation axis of the cloud is close to the mean direction of the field (e.g. see Galli et al. 2006; Mellon & Li 2008; Hennebelle & Fromang 2008). However, discs around Classes I and II young stellar objects are commonly observed (e.g. Williams & Cieza 2011; Takakuwa et al. 2012), and there is also some evidence of discs around Class 0 objects (Tobin et al. 2012; Murillo et al. 2013).

Different mechanisms have been invoked as alleviating the problem of magnetic braking during cloud collapse: (i) non-ideal magnetohydrodynamic (MHD) effects (Shu et al. 2006; Dapp & Basu 2010; Krasnopolsky et al. 2011; Braiding & Wardle 2012a,b); (ii) misalignment between the main magnetic field direction and the rotation axis (Hennebelle & Ciardi 2009; Joos et al. 2012); (iii) turbulent diffusion of the magnetic field (Seifried et al. 2012; Santos-Lima et al. 2013; Joos et al. 2013); (iv) flux redistribution driven by the interchange instability (Krasnopolsky et al. 2012); and (v) depletion of the infalling envelope anchoring the magnetic field (Mellon & Li 2009; Machida et al. 2011).

Non-ideal MHD effects, namely ambipolar, Hall, and Ohmic diffusion, depend on the abundances of charged species and on their mass and charge. The ionisation fraction, in turn, is controlled by cosmic rays (CRs) in cloud regions of relatively high column density (visual extinction A_v ≳ 4, McKee 1989) where star formation takes place. The CR ionisation rate is usually assumed to be equal to a “standard” (constant) value of ζ^H₂ ≈ 10^-17 s^-1, often called the “Spitzer” value (Spitzer & Tomasko 1968). However, CRs interacting with H₂ in a molecular cloud lose energy by several processes, mainly by ionisation losses (see Padovani et al. 2009, hereafter PGG09). As a consequence, while low-energy CRs (E ≲ 100 MeV) are possibly prevented from entering a molecular cloud because of streaming instability (Cesarsky & Völk 1978), high-energy CRs are slowed down to energies that are relevant for ionisation (ionisation cross sections for protons and electrons colliding with H₂ peak at about 100 keV and 0.1 keV, respectively).

PGG09 show that ζ^H₂ can decrease by about two orders of magnitude from “diffuse” clouds of column density ~10²¹ cm^-2 to “dense” clouds and massive envelopes with column densities of ~10²⁴ cm^-2. A similar attenuation can in principle take place during the process of cloud collapse, resulting in a decrease in ζ^H₂ in the inner region of a core where the formation of a protostellar disc is expected to occur. A further attenuation is predicted by the toroidal field component generated by rotation that increases the particles’ path length and enhances the losses by magnetic mirroring (Padovani & Galli 2011, hereafter PG11). A reduced CR ionisation rate results in a more efficient ambipolar diffusion, which may help to alleviate the magnetic braking problem (Mellon & Li 2009). The aim of this paper is to evaluate how variations in ζ^H₂ can affect the resistivity of the gas, and, eventually, how CRs influence the dynamics of collapse and the formation of a circumstellar disc.

A full treatment of this problem, in which the propagation of CRs is computed self-consistently with the evolution of density and magnetic field, following at the same time the formation and destruction of chemical species, would be prohibitively time-consuming from a numerical point of view. Our approach is therefore a simpler one. First, (a) we sacrifice the self-consistency by taking snapshots at particular evolutionary times of magnetic field configurations and density distributions from ideal MHD simulations that do not include any resistivity of the gas. Second, (b) we propagate CRs in these configurations, and compute the spatial distribution of the CR ionisation rate. Then (c) we build a simplified chemical model (that can be used as a fast subroutine in any dynamical code) to approximately evaluate the chemical composition at each spatial position using as input the distribution of CR ionisation rates determined at the previous step. Finally, (d) we compute the microscopic resistivities (ambipolar, Hall, and Ohmic) and compare the time scale of magnetic field diffusion t_B to the dynamical time scale t_dyn at each point in the model to determine the region of magnetic decoupling, where t_B<t_dyn. In this region, the dynamics of the gas is unimpeded by magnetic forces, and whatever angular momentum is carried by fluid particles crossing its border, it will be conserved in its interior. In particular, centrifugally supported discs can form inside this region (but not outside). Clearly, the assumption of ideal MHD on which the simulations are based becomes invalid (by definition) inside the decoupling region. In this sense, our calculation is not fully self-consistent. However, our aims in this preliminary investigation are: first, to show whether the extent of the decoupling region can be sufficiently large to allow the formation of a realistic disc (at least ~10 AU in radius); second, to compare the size of this region obtained with the accurate treatment of CR attenuation with column density and magnetic field developed in our previous studies and with a constant value of ζ^H₂, as usually assumed.

The paper is organised as follows. In Sect. 2, we provide a description of the model of CR propagation as well as of the numerical simulations adopted. In Sect. 3 we describe the basic features of the chemical code that we employed to compute the ionisation fractions, exploring the effect of different grain size distributions on abundances. In Sect. 4 we calculate the diffusion coefficients, and, in Sect. 5, we use them to calculate the corresponding diffusion time scales. Finally in Sect. 6 we summarise our conclusions. A full description of the chemical code is given in Appendix A. The dependence of the diffusion coefficients on the grain size distribution is examined in detail in Appendix B. Finally, in Appendix C we describe the contribution to the total diffusion time due to the different diffusion processes.

2. Description of the models

PG11 showed that the magnetic fields of dense molecular cloud cores (modelled as equilibrium configurations) can influence the penetration and propagation of CRs from the intercloud medium to the core centre, affecting the spatial distribution of the ionisation rate. In fact, since charged particles spiral around field lines, CRs propagating in complex magnetic configuration “see” a higher column density than particles moving on straight trajectories, and therefore suffer larger energy losses. In addition, Padovani et al. (2013, hereafter PHG13) found that the mirroring effect becomes stronger when the toroidal field component is larger than about 40% of the total field, in the central 300–400 AU where density is higher than 10⁹ cm^-3. This makes the CR ionisation rate, ζ^H₂ (hereafter we refer specifically to the ionisation of H₂), to drop well below 10^-18 s^-1 down to about 10^-20 s^-1, roughly equal to the ionisation level arising from the decay of long- and short-lived radionuclides within protoplanetary discs (Umebayashi & Nakano 1981; Cleeves et al. 2013). PHG13 performed a numerical study of the propagation of CRs in collapsing clouds using several snapshots of numerical simulations, adopting the formalism developed in PG11. They found that the value of ζ^H₂ can be orders of magnitude lower than the standard “Spitzer” value in the inner region of a core where the formation of a protostellar disc is expected to take place. From their numerical results, they derived an approximate analytical expression to compute the CR ionisation rate taking both column density and magnetic field effects into account. This formula provides a fast and efficient way to calculate the distribution of ζ^H₂ that can be employed to estimate the ionisation fractions of charged particles in order to determine the diffusion coefficients.

In the following, we compute the distribution of ζ^H₂ for three different configurations obtained from numerical simulations performed with the AMR code RAMSES¹ (Teyssier 2002; Fromang et al. 2006) with the goal to compute the diffusion time scales in the collapse region. The first two cases are low-mass models from Joos et al. (2012), while the third one represents a high-mass model from Commerçon et al. (in prep.). Table 1 summarises the parameters. We assumed the dependence of ζ^H₂ on column density given by model ℳ in PHG13 (see their Fig. 1). This model is obtained by adopting the CR interstellar proton and electron spectra presented in Webber (1998) and Strong et al. (2000), respectively.

2.1. Low-mass models

Models L₁ and L₂ have been already used by PHG13 to calculate maps of the CR ionisation rate (see their Figs. 7 and 10) and they represent two extreme situations: one snapshot taken at the beginning of the simulation and one at the end (L₁ and L₂, respectively). These ideal MHD numerical simulations describe the collapse of a rotating core with initial mass of 1 M_⊙ and a density profile given by a modified power-law, $\mathit{n}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{0}}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}\hspace{0.17em}\mathit{,}$(1)\noindentwhere n₀ = 7.8 × 10⁶ cm^-3 and r₀ = 4.68 × 10^-3 pc according to observations (André et al. 2000; Belloche et al. 2002). The ratio of thermal-to-gravitational energy is about 0.25 and the ratio of rotational-to-gravitational energy ratio is about 0.03. Model L₁ is a parallel rotator, i.e. the rotation axis is aligned to the magnetic axis, defined as the average direction of the field (α_B,J = 0). No rotationally supported disc has formed at the time of this particular snapshot (t = 0.924 kyr after the formation of the first Larson’s core). In this case, ζ^H₂ decreases to values of 2−4 × 10^-18 s^-1 in the inner 100–200 AU radius (see Fig. 7 in PHG13). Model L₂ is a perpendicular rotator (α_B,J = π/ 2) and a late-time configuration (t = 10.756 kyr), showing a Keplerian disc perpendicular to the rotation axis. Here ζ^H₂ drops to very low values, down to 2 × 10^-21 s^-1 in the inner few tens of AU and it is lower than 10^-18 s^-1 in a region of about 200 AU radius (see Fig. 10 in PHG13). Clearly, at these low levels, CRs compete with short-lived radionuclides in determining the ionisation fraction (see e.g. Cleeves et al. 2013).

2.2. High-mass model

Model H has similar initial conditions as in Commerçon et al. (2011). It consists in a 100 M_⊙ uniform temperature (T = 10 K) dense core, with the same density profile shape as in models L₁ and L₂ (n₀ = 1.2 × 10⁷ cm^-3, r₀ = 1.87 × 10^-2 pc, and a factor of 10 in density contrast between the centre and the border of the sphere). The turbulent-to-gravitational energy ratio is about 0.2, corresponding to an initial Mach number of 7, and the thermal-to-gravitational energy ratio is about 0.01. Following Hennebelle et al. (2011) and Commerçon et al. (2011), we apply initial perturbations to the velocity field only to account for initial turbulence in the core, which is not driven at large scales after the start of the calculations. There is no global initial rotation in the model, namely the angular momentum is built from the initial velocity fluctuations. The initial magnetic field is aligned with the x-axis and its intensity is proportional to the column density through the cloud.

In the following we focus on the densest fragment formed in the calculations whose properties are summarised in Table 1 and we discuss the results only for the (z,x) plane, since it is close to the disc plane. As model L₂, model H shows ζ^H₂ ≃ 10^-21 s^-1 in the inner few tens of AU, but in the latter case the region with ζ^H₂ ≲ 10^-19 s^-1 is even larger, with a mean radius of 150 AU. At a radius of 300 AU ζ^H₂ increases only to up to 5 × 10^-18 s^-1, while at the same radius model L₂ has already reached CR ionisation rate values larger than 10^-17 s^-1. This can be explained by noting the larger extent of the region with high density in the high-mass case with respect to the low-mass cases (see Fig. 1).

Fig. 1 CR ionisation maps and iso-density contours (black solid lines) for the case H in Table 1. Left panels: entire computational domain while right panels: a zoom in the inner region. Upper and lower panels: two perpendicular planes both containing the density peak. Labels show log 10 [ n/ cm-3 ].

3. Abundances of charged species

To compute the ionisation fraction we adopt a “minimal” chemical network (a simplified version of more extensive networks, like e.g. in Umebayashi & Nakano 1990) that computes the steady-state abundance of H⁺, ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, a typical molecular ion mH⁺ (e.g. HCO⁺), a typical “metal” ion M⁺ (e.g. Mg⁺), electrons and dust grains as a function of the H₂ density, temperature and CR ionisation rate ζ^H₂ at each spatial position in our models. For each species i, the abundance is given by x(i) ≡ n(i) /n(H₂). Neutral hydrogen is assumed to be in the form of H₂, namely n(H) = 2n(H₂). The fraction of charged vs. neutral grains is computed in a simplified way: if the density of electrons is higher than the density of grains, all grains are assumed to carry one electron; otherwise, the fraction of negatively charged grains is determined by charge balance with the positively charged species, and the residual number of grains is assumed to be neutral. Positively charged grains as well as multiple charged grains are ignored for simplicity. A detailed description of the chemical network adopted and a summary of the reactions included is given in Appendix A.

As a benchmark test, we compared our results with those obtained with the publicly available code ASTROCHEM², that includes about a thousand of reactions, finding comparable results within a factor of 2 to 5. The use of our simplified code is justified by the fact that it allows the calculation of the ionisation fractions for each point of our models, in about 0.7 ms on average against about 20 s needed by ASTROCHEM. Thus, our minimal chemical network, combined with the fitting formula given by Eqs. (19)–(24) of PHG13 to compute ζ^H₂(N,B), is an effective tool to rapidly compute the fractional abundances and diffusion coefficients at each time step in non-ideal MHD simulations.

3.1. Effects of the grain size distribution

Grains play a decisive role in determining the degree of coupling between the gas and the magnetic field. In fact, the electrical resistivity of the gas depends on the abundance and the size distribution of charged grains: larger grains have a smaller Hall parameter (see Sect. 4) than smaller grains, and therefore are less coupled to the magnetic field. To cover all possible situations, we run our chemical model for three different grain size distribution. In particular, we fix the maximum size, a_max = 3 × 10^-5 cm (Nakano et al. 2002), while we vary the minimum grain size: (i) a_min = 10^-5 cm, representative of large grains formed by compression and coagulation during the collapse (Flower et al. 2005); (ii) a_min = 10^-6 cm, the minimum grain radius of a MRN size distribution (Mathis et al. 1977) that gives the same grain opacity found by Flower et al. (2005); and (iii) a_min = 10^-7 cm, a typical size for very small grains.

In Fig. 2 we compare the abundances computed with our minimal chemical network assuming a constant ζ^H₂ = 5 × 10^-17 s^-1 (hereafter “constant-ζ”) with those obtained from the spatially-resolved values of ζ^H₂ computed by PHG13 (hereafter “variable-ζ”). The comparison is shown only for the model L₂ (models L₁ and H give similar results). As a general remark, independently of a_min, the variable-ζ model gives larger abundances of charged species than the constant-ζ model at n ≲ 10⁶ cm^-3. Conversely, at higher densities, abundances from the variable-ζ model are well below those resulting from the constant-ζ model. This happens because the variable-ζ model is higher than the constant-ζ at n ≲ 10⁶ cm^-3, while it quickly decreases below ζ^H₂ = 5 × 10^-17 s^-1 at higher densities.

Fig. 2 Ionisation fractions for model L2 as a function of the volume density computed with a constant CR ionisation rate ζH2= 5 × 10-17 s-1 (constant-ζ, solid lines) and with a spatially resolved ζH2 (variable-ζ, dots). The three panels are for amin = 10-5 cm (top), amin = 10-6 cm (middle) and amin = 10-7 cm (bottom).

3.2. Dependence of chemical abundances on ζH2

The relation between densities of charged species and neutrals is usually expressed in the form $\mathit{n}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{\propto }\mathrm{(}{\mathit{\zeta }}^{{\mathrm{H}}_{\mathrm{2}}}{\mathrm{)}}^{{\mathit{k}}^{\mathrm{\prime }}}\mathit{,}$(2)with k′ ≈ 1 / 2 (e.g. see Ciolek & Mouschovias 1994; 1995). However, k′ can differ from 1 / 2 depending on the grain size (except for molecular ions). Figure 3 shows the dependence of the chemical species computed with our minimal model as a function of ζ^H₂. As shown by the figure, in the case of large grains the abundance of electrons and metal ions follows Eq. (2) with k′ = 1 / 2, while negatively charged grains are independent of ζ^H₂. For small grain size k′ increases towards 1 and 1 / 2 for electrons and negative grains, respectively, whereas metal ions become independent of ζ^H₂. It is worth noting that in the case of strong depletion, ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and H⁺ are the most abundant species and both $\mathit{x}\mathrm{\left(}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\right)}$ and x(H⁺) are proportional to ζ^H₂ irrespective of grain size and density.

Fig. 3 Ionisation fractions as a function of the CR ionisation rate at density n(H2) = 106 cm-3 assuming large grains (amin = 10-5 cm, solid lines) and small grains (amin = 10-7 cm, dotted lines). Dashed grey lines show the trend for xi ∝ ζH2 and xi ∝ (ζH2)1 / 2.

4. Diffusion coefficients

The electrical resistivity of a plasma is a measure of the ability of the magnetic field and the charges attached to it to move (diffuse) with respect to the neutrals and/or other charges. In the weakly ionised gas that characterise dense cores, the resistivity is dominated by ambipolar diffusion, a process by which charged particles develop a drift velocity with respect to the neutral component, and the Lorentz force acting on the charges is conveyed to the neutral gas through collisions.

For each charged species i in a sea of neutrals (molecular hydrogen), the parameter gauging the relative importance of the Lorentz and drag forces is the Hall parameter (e.g. Wardle & Ng 1999) defined by ${\mathit{\beta }}_{\mathit{i,}{\mathrm{H}}_{\mathrm{2}}}\mathrm{=}\left(\frac{{\mathit{Z}}_{\mathit{i}}\mathit{eB}}{{\mathit{m}}_{\mathit{i}}\mathit{c}}\right)\frac{{\mathit{m}}_{\mathit{i}}\mathrm{+}{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}}}{\mathit{\mu }{\mathit{m}}_{\mathrm{H}}\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}\mathrm{⟨}\mathit{\sigma v}{\mathrm{⟩}}_{\mathit{i,}{\mathrm{H}}_{\mathrm{2}}}}\mathit{,}$(3)where m_i and Z_ie are the mass and the charge of the species i, respectively, and μ = 2.36 is the molecular weight for the assumed fractional abundances of H₂ and He. The momentum transfer rate coefficients ⟨ σv ⟩ _i,H2 have been parameterised as a function of temperature and relative speed in Pinto & Galli (2008).

Fig. 4 Density contours (solid blue lines with labels indicating log 10 [ n/ cm-3 ]) superposed to the velocity field (red arrows). The shaded areas show regions dominated by ambipolar (white), Hall (grey), and Ohmic (black) diffusion. The resistivities are calculated for three different values of amin and for the three models L1, L2, and H.

Drifts of charged species with respect to neutrals determines different regimes for the magnetic diffusivity. The induction equation then becomes $\begin{array}{ccc}\frac{\mathit{\partial }{B}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{\times }\mathrm{(}{B}\mathrm{\times }{U}\mathrm{)}& \mathrm{=}& \mathrm{\nabla }\mathrm{\times }\{{\mathit{\eta }}_{\mathrm{O}}\mathrm{\nabla }\mathrm{\times }{B}\mathrm{+}{\mathit{\eta }}_{\mathrm{H}}\mathrm{(}\mathrm{\nabla }\mathrm{\times }{B}\mathrm{)}\mathrm{\times }\frac{{B}}{\mathit{B}}\\ & & \mathrm{+}{\mathit{\eta }}_{\mathrm{AD}}\left[\mathrm{(}\mathrm{\nabla }\mathrm{\times }{B}\mathrm{)}\mathrm{\times }\frac{{B}}{\mathit{B}}\right]\mathrm{\times }\frac{{B}}{\mathit{B}}\}\mathit{,}\end{array}$(4)where U is the fluid velocity and B the magnetic field vector. Ambipolar, Hall, and Ohmic resistivities (η_AD, η_H, and η_O, respectively) can be written as a function of the parallel (σ_∥), Pedersen (σ_P) and Hall (σ_H) conductivities (e.g. see Wardle 2007; Pinto et al. 2008) $\begin{array}{ccc}& & {\mathit{\eta }}_{\mathrm{AD}}\mathrm{=}\frac{{\mathit{c}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }}\left(\frac{{\mathit{\sigma }}_{\mathrm{P}}}{{\mathit{\sigma }}_{\mathrm{P}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{H}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{\parallel }}}\right)\mathit{,}\\ & & {\mathit{\eta }}_{\mathrm{H}}\mathrm{=}\frac{{\mathit{c}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }}\left(\frac{{\mathit{\sigma }}_{\mathrm{H}}}{{\mathit{\sigma }}_{\mathrm{P}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{H}}^{\mathrm{2}}}\right)\mathit{,}\\ & & {\mathit{\eta }}_{\mathrm{O}}\mathrm{=}\frac{{\mathit{c}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }{\mathit{\sigma }}_{\mathrm{\parallel }}}\mathit{,}\end{array}$which are defined by $\begin{array}{ccc}& & {\mathit{\sigma }}_{\mathrm{\parallel }}\mathrm{=}\frac{\mathit{ecn}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}}{\mathit{B}}\sum _{\mathit{i}}{\mathit{Z}}_{\mathit{i}}{\mathit{x}}_{\mathit{i}}{\mathit{\beta }}_{\mathit{i,}{\mathrm{H}}_{\mathrm{2}}}\mathit{,}\\ & & {\mathit{\sigma }}_{\mathrm{P}}\mathrm{=}\frac{\mathit{ecn}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}}{\mathit{B}}\sum _{\mathit{i}}\frac{{\mathit{Z}}_{\mathit{i}}{\mathit{x}}_{\mathit{i}}{\mathit{\beta }}_{\mathit{i,}{\mathrm{H}}_{\mathrm{2}}}}{\mathrm{1}\mathrm{+}{\mathit{\beta }}_{\mathit{i,}{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{2}}}\hspace{0.17em}\mathit{,}\\ & & {\mathit{\sigma }}_{\mathrm{H}}\mathrm{=}\frac{\mathit{ecn}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}}{\mathit{B}}\sum _{\mathit{i}}\frac{{\mathit{Z}}_{\mathit{i}}{\mathit{x}}_{\mathit{i}}}{\mathrm{1}\mathrm{+}{\mathit{\beta }}_{\mathit{i,}{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{2}}}\mathrm{·}\end{array}$In general, the ambipolar resistivity term controls diffusion in low density regions (n_H2 ≲ 10⁸−10⁹ cm^-3), whereas Hall diffusion dominates at intermediate densities (10⁸−10⁹ cm^-3 ≲ n_H2 ≲ 10¹¹ cm^-3) and Ohmic dissipation sets in at even higher densities (n_H2 ≳ 10¹¹ cm^-3), see Umebayashi & Nakano (1981). However, the extent to which a diffusion process dominates over the others hinges on several factors, one of which is the assumed grain size distribution. To show this effect, we computed the ionisation fractions using the values of ζ^H₂ for the three models (L₁, L₂, and H) and we compared the spatial distribution of the corresponding resistivities varying a_min. Background colours in Fig. 4 show the predominant diffusion mechanism. The distributions with a_min = 10^-5 and 10^-7 cm lead to similar results, independently on the distribution of ζ^H₂. This can be explained by looking at the trend of the different resistivity contributions as a function of the grain size and n(H₂) (see Appendix B). In fact, the ambipolar diffusion term does not increase monotonically with grain radius but, at densities of 10⁸−10⁹ cm^-3, it shows an absolute minimum around a_min = 10^-6 cm, while the Hall term continues to grow. At the highest densities the Ohmic term prevails, except for model L₁. The Hall term can play an important role down to densities of about 10⁸ cm^-3 for a_min = 10^-6 cm.

5. Diffusion time scales

The drift velocity of the magnetic field can be represented by the velocity of the charged species, which are frozen with field lines, with respect to neutrals. From the comparison of this velocity with the fluid velocity, it is possible to assess the degree of diffusion of the field and then to estimate the size of the region where gas and magnetic field are decoupled. In principle, a more direct estimate of the reduction of magnetic braking resulting from a decrease in the CR ionisation rate would come from a comparison between the magnetic braking time with the dynamical time of the flow. However, a simple estimate of the former is not straightforward, as it depends crucially on the field morphology, strength, and relative orientation of the average field direction and the cloud’s angular momentum. For this reason, we prefer to present a comparison, at each spatial position at a given time step, between the diffusion time of the field and the dynamical time of the flow as defined in the following.

The magnetic field drift velocity U_B can be written as a function of resistivities (Nakano et al. 2002), allowing to isolate the ambipolar (AD), Hall (H), and Ohmic (O) contributions, namely $\begin{array}{ccc}{{U}}_{\mathit{B}}\mathrm{=}{{U}}_{\mathrm{AD}}\mathrm{+}{{U}}_{\mathrm{H}}\mathrm{+}{{U}}_{\mathrm{O}}\mathit{,}& & \end{array}$(11)where

$\begin{array}{ccc}& & {{U}}_{\mathrm{AD}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }{\mathit{\eta }}_{\mathrm{AD}}}{\mathit{c}{\mathit{B}}^{\mathrm{2}}}{j}\mathrm{\times }{B}\mathit{,}\\ & & {{U}}_{\mathrm{H}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }{\mathit{\eta }}_{\mathrm{H}}}{\mathit{c}{\mathit{B}}^{\mathrm{3}}}\mathrm{(}{j}\mathrm{\times }{B}\mathrm{)}\mathrm{\times }{B}\mathit{,}\\ & & {{U}}_{\mathrm{O}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }{\mathit{\eta }}_{\mathrm{O}}}{\mathit{c}{\mathit{B}}^{\mathrm{2}}}{j}\mathrm{\times }{B}\hspace{0.17em}\end{array}$and ${{U}}_{\mathit{B}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathit{c}{\mathit{B}}^{\mathrm{2}}}\left[\left({\mathit{\eta }}_{\mathrm{AD}}\mathrm{+}{\mathit{\eta }}_{\mathrm{O}}\right){j}\mathrm{\times }{B}\mathrm{+}{\mathit{\eta }}_{\mathrm{H}}\left({j}\mathrm{\times }{B}\right)\mathrm{\times }\frac{{B}}{\mathit{B}}\right]\mathit{,}$(15)where j = (c/ 4π)∇ × B is the current density. Thus, the diffusion time of the magnetic field, t_B, can be written as a function of the time scales associated to the three diffusion processes, $\frac{\mathrm{1}}{{\mathit{t}}_{\mathit{B}}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{AD}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{H}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{O}}}\mathit{,}$(16)where t_k = R/U_k (k = AD,H,O) and R a typical length scale of the region (in the following at each point we take R equal to the distance from the density peak). The diffusion time of the magnetic field can then be compared to the time scale of evolution of the fluid (for example, Nakano et al. 2002 compare t_B to the free-fall time of a spherical homogeneous cloud). In this work we define the dynamical time scale of the cloud as t_dyn = R/U, where U is the fluid velocity, including both infall and rotation. In regions where t_B<t_dyn the magnetic field is partially decoupled and therefore has less influence on the gas dynamics while, if t_B>t_dyn, diffusion is not efficient enough and the magnetic field remains well coupled to the gas. To stress the importance of properly taking the propagation of CRs inside a molecular cloud into account, we also compare the dynamical time with the diffusion time computed in the constant-ζ case assuming ζ^H₂= 5 × 10^-17 s^-1 (hereafter t_{B,ζ − const}), and in the variable-ζ case, taking ζ^H₂ from PHG13 (hereafter t_{B,ζ − var}).

5.1. Low-mass models

Figure 5 shows contours of diffusion and dynamical time scales for the aligned rotator model L₁. Although the disc is not formed at this early stage (t = 0.824 kyr), for a_min = 10^-5 cm in the variable-ζ case the magnetic diffusion time becomes shorter than the dynamical time in a central region with a radius of about 20 AU. For smaller grains, t_{B,ζ − var} is always larger than t_dyn and no decoupling zone is formed. The diffusion time computed in the constant-ζ case is about one order of magnitude longer. Even at later times (t = 11.025 kyr), Joos et al. (2012) find no disc formation for this aligned rotator simulation. However, the low ionisation rate found by PHG13 in a small central region of density higher than ~10⁹ cm^-3 (see their Fig. 8) can presumably promote the formation of a centrifugally supported disc via enhanced Ohmic and Hall diffusion, at least in the case of large grains.

Fig. 5 Diffusion time contours (white solid lines) evaluated with constant ζH2= 5 × 10-17 s-1 (left column) and ζH2 from model L1 in Table 1 (right column) compared with dynamical time contours (black dashed lines) for three different values of the minimum grain size. Labels show log 10(t/ yr).

A more interesting case is shown by the model L₂ (Fig. 6), where a protostellar disc of radius ~200 AU has formed. While t_{B,ζ − const} is always larger than t_dyn, t_{B,ζ − var} is lower than t_dyn inside a decoupling zone whose size depends on a_min. Magnetic decoupling is favoured by large grains (a_min = 10^-5 cm) that are expected to form by coagulation of smaller grains during the collapse. In particular, the lower left panel shows that in the inner region with a radius of about 50 AU the gas experiences collapse while the magnetic field diffuses. The region of decoupling shrinks with decreasing a_min, but even with the smallest a_min = 10^-7 cm the decoupling occurs inside a region of about 20 AU of radius. This comparison definitely proves that CRs play a crucial role in determining the protostellar collapse time scale. In fact, the correct evaluation of the CR ionisation rate as a function of density and magnetic field allows the diffusion time to decrease up to three order of magnitudes. In Appendix C we also show separately the contribution to t_B resulting from ambipolar, Hall, and Ohmic diffusion, respectively.

Fig. 6 Diffusion time contours (white solid lines) evaluated with constant ζH2= 5 × 10-17 s-1 (left column) and ζH2 from model L2 in Table 1 (right column) compared with dynamical time contours (black dashed lines) for three different values of the minimum grain size. Labels show log 10(t/ yr).

5.2. High-mass model

The high-mass case (Fig. 7) shows a similar behaviour to model L₂, but in the case of a_min = 10^-5 cm, the decoupling between gas and field is allowed in a even larger region, of about 100 AU of radius where both t_{B,ζ − var} and t_{B,ζ − const} are lower than t_dyn. Decreasing the grain size, the region of decoupling becomes narrower and vanishes for a_min = 10^-7 cm. The lack of a decoupling region in this high-mass case is at variance with the low-mass case discussed in the previous Section. This depends on the field diffusing faster in the high-mass case owing to the turbulent nature of the flow, as already stressed by Hennebelle et al. (2011) and Joos et al. (2013). This diffusion is largely due to numerical effects associated to the finite spatial resolution of the simulation. This “numerical” field diffusion results in a reduced strength of the magnetic field with respect to the (non-turbulent) low-mass case, and this in turn makes the microscopic resistivity smaller (especially ambipolar diffusion, which is proportional to B²). In this case, the microscopic resistivity is not sufficient to produce a significant decoupling region. Of course, if turbulence actually affects the field diffusion, the microscopic resistivity is not relevant.

Fig. 7 Diffusion time contours (white solid lines) evaluated with constant ζH2= 5 × 10-17 s-1 (left column) and ζH2 from model H in Table 1 (right column) compared with dynamical time contours (black dashed lines) for three different values of the minimum grain size. Labels show log 10(t/ yr).

6. Conclusions

In this paper we investigated the role of CRs in the formation of protostellar discs, using the results of PHG13 to evaluate CR ionisation rate in collapsing low-mass clouds accounting for both column density and magnetic effects. We developed a simple chemical code based on a minimal set of reactions to calculate, at each position in space at any given time, the abundances of the charged species, the corresponding electric resistivities, and the times scales associated to the main magnetic diffusion processes (ambipolar, Hall, and Ohmic). We applied this formalism to selected snapshots of numerical MHD simulations of the collapse of low- and high-mass clouds. Comparing the magnetic diffusion and the dynamical time scale, we determined the extent of the area where the gas is dynamically decoupled from the field. Inside this region, field freezing is invalid, and magnetic braking is ineffective.

We performed our calculations in two cases: assuming a spatially uniform CR ionisation rate of ζ^H₂ = 5 × 10^-17 s^-1, and adopting the formalism of PHG13 to evaluate the attenuation of CRs in a magnetised cloud. We found that the ionisation fraction is significantly lower in the second case, and therefore the coupling with the magnetic field is weaker than usually assumed in the central region of a collapsing cloud. In particular:

• 1.

  For model L₂, a late-time snapshot of the collapse of a low-mass cloud with mean magnetic field perpendicular to the rotation axis, a decoupling zone of radius of 50–100 AU around the central protostar is found in the case of variable ζ^H₂, but not when ζ^H₂ is assumed constant. This size compares well with the size of protostellar discs. This stresses the importance of accounting for CR attenuation when computing ionisation fractions.

• 2.

  No decoupling zone is found in model L₁, an early time snapshot of the collapse of a low-mass cloud with aligned field and rotation axis. In this model, the time scale of field diffusion remains longer than the dynamical time scale everywhere, either for a uniform or attenuated CR ionisation rate. In this case, the relatively modest increase in the density and field strength by compression are insufficient to attenuate the CR flux to levels low enough to produce significant decoupling.

• 3.

  A decoupling zone of size ~100 AU is also found in the case of model H, a snapshot from a typical high-mass collapse simulation. However, its size becomes smaller for smaller grains, and disappears altogether for grains smaller than 10^-6 cm.

Although the models adopted do not represent a time sequence, they nevertheless suggest that a decrease in the ionisation and/or an increase in the resistivity occurs in the innermost region of a cloud some time after the onset of collapse, but not earlier. In fact, the conditions for a substantial increase in the magnetic diffusion time are that the field is considerably twisted (compare the late-time misaligned model L₂ to the early-time aligned model L₁) and that dust grains had time to grow by coagulation (compare the upper and lower panels of Figs. 5–7). It is tempting to speculate that large, 100 AU-size discs are only allowed to form at a later stage when the powerful magnetic brake on the infalling gas has been relieved by either (or a combination) of these effects.

Other results of our study are:

• 1.

  The dominant diffusion processes are ambipolar and Halldiffusion, with Ohmic resistivity becoming important only at thehighest densities reached in our models (n ≳ 10¹¹ cm^-3), where η_O is a factor of few larger than η_AD and η_H, (see black regions in Fig. 4). In general, ambipolar and Hall diffusion are comparable over the innermost few hundred AU around a forming protostar. Thus, Hall diffusion should not be neglected in non-ideal MHD collapse calculations. However, the region where a diffusion process dominates over the others depends sensitively on the grain size distribution. In fact, while Hall and Ohmic resistivities have a monotonic dependence on the grain size (all else being equal), ambipolar diffusion has a minimum for a mean grain radius of about 1−3 × 10^-6 cm especially at high densities.

• 2.

  In general, the size of the decoupling zone decreases for smaller grain size. The maximum grain size assumed here, (a_min = 10^-5 cm), is likely a realistic value for the condition expected in disc-forming regions, where larger grains are predicted to form by compression and coagulation of smaller grains. Both Models L₂ and H, referring to low- and high-mass cases, respectively, predict the formation of a protostellar disc, but for model H the decoupling region disappears for smaller grains (a_min = 10^-7 cm), since the magnetic field becomes weaker.

• 3.

  The dependence of the fractional abundances of charged species on the CR ionisation rate may differ from the usually assumed $\sqrt{{\mathit{\zeta }}^{{\mathrm{H}}_{\mathrm{2}}}}$ dependence. Our minimal chemical model shows that depending on the grain size distribution, the abundance of several species can scale as ζ^H₂ or become independent on ζ^H₂, and in general the dependence is not the same for all charged species. Thus, care should be taken when the ionisation fraction is scaled with the CR flux.

To summarise, we demonstrated that a correct treatment of CR propagation can explain the occurrence of a decoupling region between gas and magnetic field that in turn affects the disc formation. We emphasise, however, that our calculations are not strictly self-consistent, because we computed microscopic resistivities using the density and magnetic field strength obtained from ideal MHD simulations. If diffusion processes were included self-consistently in the numerical simulation itself, the line twisting would be presumably reduced. This would attenuate the effect of cosmic-ray mirroring, and the decrease in ζ^H₂ would not be so strong as found in this paper. In this sense, our study should be considered a proof of concept showing how a correct evaluation of ζ^H₂ can affect the protostellar disc formation.

Online material

Appendix A: Chemical model

In this Appendix we describe the “minimal” chemical model introduced in Sect. 3 to compute the ionisation fraction. The charged species considered are H⁺, ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, molecular ions mH⁺ (e.g. HCO⁺), “metal” ions M⁺ (e.g. Mg⁺), electrons e, and negatively charged dust grains, whereas the neutral species are H₂, heavy molecules m (e.g. CO), metal atoms M, (e.g. Mg) and neutral grains. Charged and neutral grains are collectively indicated as g. We indicate with x(i) the abundance of each species i with respect to H₂. The abundance of the neutral species is fixed. In particular, we assume x(m) ≃ 6 × 10^-4 and x(M) ≃ 4 × 10^-8. All rate coefficients are estimated at T = 10 K.

Appendix A.1: Minimal chemical network

Protons are produced by CR ionisation of H₂ at a rate ϵζ^H₂, with ϵ ≃ 0.05 (Shah & Gilbody 1982). They are mainly destroyed by charge transfer (CT) with molecules (at a rate β ≃ 10^-9 cm³ s^-1) and by recombination on grains (at a rate α_gr, see Eq. (A.8)) $\mathit{ϵ}{\mathit{\zeta }}^{{\mathrm{H}}_{\mathrm{2}}}\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}\mathrm{\left[}\mathit{\beta n}\mathrm{\right(}\mathit{m}\mathrm{)}\mathrm{+}{\mathit{\alpha }}_{\mathrm{gr}}\mathit{n}\mathrm{(}\mathit{g}\mathrm{\left)}\mathrm{\right]}\mathit{n}\mathrm{\left(}{\mathrm{H}}^{\mathrm{+}}\mathrm{\right)}\mathit{.}$(A.1)The formation of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ is driven by CR ionisation of H₂ at a rate (1 − ϵ)ζ^H₂, while destruction is due to CT with heavy molecules, dissociative recombination (DR, at a rate α_dr ≃ 10^-6 cm³ s^-1), and recombination on grains $\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}{\mathit{\zeta }}^{{\mathrm{H}}_{\mathrm{2}}}\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}\mathrm{\left[}\mathit{\beta n}\mathrm{\right(}\mathit{m}\mathrm{)}\mathrm{+}{\mathit{\alpha }}_{\mathrm{dr}}\mathit{n}\mathrm{(}\mathit{e}\mathrm{)}\mathrm{+}{\mathit{\alpha }}_{\mathrm{gr}}\mathit{n}\mathrm{(}\mathit{g}\mathrm{\left)}\mathrm{\right]}\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\right)}\mathit{.}$(A.2)The formation of molecular ions mH⁺ occurs by CT of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and heavy molecules, while destruction occurs by DR and recombination on grains $\mathit{\beta n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\right)}\mathit{n}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}\mathrm{\left[}{\mathit{\alpha }}_{\mathrm{dr}}\mathit{n}\mathrm{\right(}\mathit{e}\mathrm{)}\mathrm{+}{\mathit{\alpha }}_{\mathrm{gr}}\mathit{n}\mathrm{(}\mathit{g}\mathrm{\left)}\mathrm{\right]}\mathit{n}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{.}$(A.3)Metal ions are formed by CT with ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and mH⁺, and destroyed by recombination with free electrons and on grains $\mathit{\beta n}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{\left[}\mathit{n}\mathrm{\right(}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{)}\mathrm{+}\mathit{n}\mathrm{(}\mathit{m}{\mathrm{H}}^{\mathrm{+}}\mathrm{\left)}\mathrm{\right]}\mathrm{=}\mathrm{\left[}{\mathit{\alpha }}_{\mathrm{rec}}\mathit{n}\mathrm{\right(}\mathit{e}\mathrm{)}\mathrm{+}{\mathit{\alpha }}_{\mathrm{gr}}\mathit{n}\mathrm{(}\mathit{g}\mathrm{\left)}\mathrm{\right]}\mathit{n}\mathrm{\left(}{\mathit{M}}^{\mathrm{+}}\mathrm{\right)}\mathit{.}$(A.4)Note that CT with metal atoms can be neglected with respect to DR if x(e) ≫ (β/α_dr)x(M) ≃ 10^-3x(M).

Dust grains are assumed to be negatively charged (charge −1) or neutral. The total number density of grains is obtained from the MRN size distribution (Mathis et al. 1977) $\frac{\mathrm{d}\mathit{n}\mathrm{\left(}\mathit{g}\mathrm{\right)}}{\mathrm{d}\mathit{a}}\mathrm{=}\mathit{C}{\mathit{a}}^{-3.5}\mathit{,}$(A.5)between a minimum (a_min) and a maximum (a_max) grain radius (see also Sect. 3.2). The normalisation constant C is obtained by imposing that the mass density of grains is equal to a fraction q = 0.01 of the gas density. We assume that grains are spherical and have density ρ_g = 2 g cm^-3 (Flower et al. 2005). For a_max ≫ a_min we obtain $\mathit{C}\mathrm{=}\frac{\mathrm{3}\mathit{q}{\mathit{m}}_{\mathrm{H}}}{\mathrm{4}\mathit{\pi }{\mathit{\rho }}_{\mathrm{g}}{\mathit{a}}_{\max }^{\mathrm{0.5}}}\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}\mathit{.}$(A.6)Under these assumptions, the number density of grains is $\mathit{n}\mathrm{\left(}\mathit{g}\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}\mathit{q}{\mathit{m}}_{\mathrm{H}}}{\mathrm{10}\mathit{\pi }{\mathit{\rho }}_{\mathrm{g}}{\mathit{a}}_{\max }^{\mathrm{0.5}}{\mathit{a}}_{\min }^{\mathrm{2.5}}}\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}\mathit{,}$(A.7)and is strongly dependent on a_min.

The coefficient of recombination of positive ions on negatively charged grains was computed by Draine & Sutin (1987) assuming the MRN size distribution, $\begin{array}{ccc}{\mathit{\alpha }}_{\mathrm{gr}}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\frac{\mathrm{10}}{\mathrm{3}}{\mathit{e}}^{\mathrm{2}}{\mathit{a}}_{\min }{\left(\frac{\mathrm{8}\mathit{\pi }}{{\mathit{m}}_{\mathrm{i}}\mathit{kT}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\left[\mathrm{1}\mathrm{+}\frac{\mathrm{3}\mathit{kT}{\mathit{a}}_{\min }}{\mathrm{2}{\mathit{e}}^{\mathrm{2}}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\psi }\mathrm{)}\right]\\ & \mathrm{=}& \hspace{0.17em}\mathrm{1.6}\mathrm{\times }{\mathrm{10}}^{-7}\left(\frac{{\mathit{a}}_{\min }}{\AA }\right){\left(\frac{\mathit{T}}{\mathrm{10}\mathrm{K}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\\ & & \hspace{0.17em}\mathrm{\times }\hspace{0.17em}\left[\mathrm{1}\mathrm{+}\mathrm{3.6}\mathrm{\times }{\mathrm{10}}^{-4}\left(\frac{{\mathit{a}}_{\min }}{\AA }\right)\left(\frac{\mathit{T}}{\mathrm{10}\mathrm{K}}\right)\left(\frac{\mathrm{1}\mathrm{-}\mathit{\psi }}{\mathrm{4}}\right)\right]\mathit{,}\end{array}$(A.8)where ψ is a numerical coefficient equal to −2.5 for an e–H⁺ plasma and −3.8 for a heavy-ion plasma. For simplicity, we adopt the same value of α_gr for all positively charged ions, assuming a typical ion mass m_i = 25m_H. The actual value of α_gr is larger by a factor 5 and 3 for H⁺ and H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$, respectively.

Appendix A.2: Evaluation of the ionisation fractions

The equation of charge neutrality is $\mathit{n}\mathrm{\left(}\mathit{g}\mathrm{\right)}\mathrm{+}\mathit{n}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathrm{=}\mathit{n}\mathrm{\left(}{\mathrm{H}}^{\mathrm{+}}\mathrm{\right)}\mathrm{+}\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\right)}\mathrm{+}\mathit{n}\mathrm{\left(}\mathit{m}{\mathrm{H}}^{\mathrm{+}}\mathrm{\right)}\mathrm{+}\mathit{n}\mathrm{\left(}{\mathit{M}}^{\mathrm{+}}\mathrm{\right)}\mathit{.}$(A.9)Equations (A.1)–(A.4) then become

$\begin{array}{ccc}& & \mathit{x}\mathrm{\left(}{\mathrm{H}}^{\mathrm{+}}\mathrm{\right)}\mathrm{=}\frac{\mathit{ϵA}}{\mathit{q}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{+}\mathit{rx}\mathrm{\left(}\mathit{g}\mathrm{\right)}}\mathit{,}\\ & & \mathit{x}\mathrm{\left(}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\mathit{A}}{\mathit{q}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{+}\mathit{rx}\mathrm{\left(}\mathit{g}\mathrm{\right)}\mathrm{+}\mathit{x}\mathrm{\left(}\mathit{e}\mathrm{\right)}}\mathit{,}\\ & & \mathit{x}\mathrm{\left(}\mathit{m}{\mathrm{H}}^{\mathrm{+}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\mathit{A}}{\mathrm{\left[}\mathit{rx}\mathrm{\right(}\mathit{g}\mathrm{)}\mathrm{+}\mathit{x}\mathrm{(}\mathit{e}\mathrm{\left)}\mathrm{\right]}\mathrm{\left[}\mathit{q}\mathrm{\right(}\mathit{m}\mathrm{)}\mathrm{+}\mathit{rx}\mathrm{(}\mathit{g}\mathrm{)}\mathrm{+}\mathit{x}\mathrm{(}\mathit{e}\mathrm{\left)}\mathrm{\right]}}\mathit{q}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{,}\\ & & \mathit{x}\mathrm{\left(}{\mathit{M}}^{\mathrm{+}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\mathit{A}}{\mathrm{\left[}\mathit{rx}\mathrm{\right(}\mathit{g}\mathrm{)}\mathrm{+}\mathit{sx}\mathrm{(}\mathit{e}\mathrm{\left)}\mathrm{\right]}\mathrm{\left[}\mathit{rx}\mathrm{\right(}\mathit{g}\mathrm{)}\mathrm{+}\mathit{x}\mathrm{(}\mathit{e}\mathrm{\left)}\mathrm{\right]}}\mathit{q}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathit{,}\end{array}$where we have defined $\begin{array}{ccc}& & \mathit{A}\mathrm{=}\frac{{\mathit{\zeta }}^{{\mathrm{H}}_{\mathrm{2}}}}{\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}{\mathit{\alpha }}_{\mathrm{dr}}}\mathit{,q}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}\frac{\mathit{\beta x}\mathrm{\left(}\mathit{m}\mathrm{\right)}}{{\mathit{\alpha }}_{\mathrm{dr}}}\mathit{,q}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{=}\frac{\mathit{\beta x}\mathrm{\left(}\mathit{M}\mathrm{\right)}}{{\mathit{\alpha }}_{\mathrm{dr}}}\mathit{,}\\ & & \mathit{r}\mathrm{=}\frac{{\mathit{\alpha }}_{\mathrm{gr}}}{{\mathit{\alpha }}_{\mathrm{dr}}}\mathit{,s}\mathrm{=}\frac{{\mathit{\alpha }}_{\mathrm{rec}}}{{\mathit{\alpha }}_{\mathrm{dr}}}\mathrm{·}\end{array}$(A.14)Thus, the equation of charge neutrality $\mathit{x}\mathrm{\left(}\mathit{g}\mathrm{\right)}\mathrm{+}\mathit{x}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathrm{=}\mathit{x}\mathrm{\left(}{\mathrm{H}}^{\mathrm{+}}\mathrm{\right)}\mathrm{+}\mathit{x}\mathrm{\left(}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\right)}\mathrm{+}\mathit{x}\mathrm{\left(}\mathit{m}{\mathrm{H}}^{\mathrm{+}}\mathrm{\right)}\mathrm{+}\mathit{x}\mathrm{\left(}{\mathit{M}}^{\mathrm{+}}\mathrm{\right)}\mathit{,}$(A.15)becomes $\begin{array}{ccc}\mathit{x}\mathrm{\left(}\mathit{g}\mathrm{\right)}\mathrm{+}\mathit{x}\mathrm{\left(}\mathit{e}\mathrm{\right)}\hspace{0.17em}& \mathrm{=}& \hspace{0.17em}\frac{\mathit{ϵA}}{\mathit{q}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{+}\mathit{rx}\mathrm{\left(}\mathit{g}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\mathit{A}}{\mathit{rx}\mathrm{\left(}\mathit{g}\mathrm{\right)}\mathrm{+}\mathit{x}\mathrm{\left(}\mathit{e}\mathrm{\right)}}\\ & & \end{array}$(A.16)We solve Eq. (A.16) as a cubic equation for x(e) assuming that all grains have charge −1. When x(e) becomes negative, we set x(e) = 0 and we solve Eq. (A.16) for x(g).

Appendix B: Dependence of diffusion coefficients on grain radius

Fig. B.1 Ambipolar (upper panel), Hall (middle panel), and Ohmic (lower panel) diffusion coefficients as a function of the mean grain radius computed at different molecular hydrogen densities: 105 (thick solid line), 106 (thick dashed line), 107 (thick dash-dotted line), 108 (dotted line), 109 (thin solid line), 1010 (thin dashed line), 1011 cm-3 (thin dash-dotted line), and 1012 cm-3 (thin dotted line).

In this Appendix we compute the resistivities using Eqs. (5)–(7) varying the minimum radius of the grain size distribution and fixing the H₂ density to evaluate the sensitivity of the diffusive terms to the grain size. For the magnetic field strength we assume | B | = (n/ cm^-3)^0.47μG (Crutcher 1999), in order to be independent on specific models. The MRN size distribution given by Eq. (A.5) implies that the mean value of the square of the grain radius weighted on the grain distribution (namely the quantity that enters the equation for the momentum transfer rate coefficients), defined by $\mathrm{⟨}{\mathit{a}}_{\mathrm{gr}}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\frac{{\mathrm{\int }}_{{\mathit{a}}_{\min }}^{{\mathit{a}}_{\max }}{\mathit{a}}^{\mathrm{2}}\mathit{C}{\mathit{a}}^{-3.5}\mathrm{d}\mathit{a}}{{\mathrm{\int }}_{{\mathit{a}}_{\min }}^{{\mathit{a}}_{\max }}\mathit{C}{\mathit{a}}^{-3.5}\mathrm{d}\mathit{a}}\mathrm{=}\mathrm{5}\frac{{\mathit{a}}_{\max }^{-0.5}\mathrm{-}{\mathit{a}}_{\min }^{-0.5}}{{\mathit{a}}_{\max }^{-2.5}\mathrm{-}{\mathit{a}}_{\min }^{-2.5}}\mathit{,}$(B.1)

is close to ${\mathit{a}}_{\min }^{\mathrm{2}}$. We vary the minimum grain radius between 10^-7 to 10^-5 cm, fixing a_max = 3 × 10⁵ cm. This corresponds to values of $\mathrm{⟨}{\mathit{a}}_{\mathrm{gr}}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ ranging from 2.2 × 10^-7 to 1.5 × 10^-5 cm. Figure B.1 shows that the coefficient of ambipolar diffusion is not monotonic with the grain radius but presents an absolute minimum at larger radii with increasing H₂ densities. On the other hand, the Hall term increases with the grain radius for any value of n(H₂), starting to decrease only at very large grain size (a_min ≳ 10^-5 cm). Finally, the Ohmic resistivity becomes important at high densities (n ≳ 10¹¹ cm^-3) and increases monotonically with grain size.

Appendix C: Contributions to the diffusion time of the magnetic field

Fig. C.1 Ambipolar, Hall, and Ohmic contribution to the diffusion time (red dotted lines) compared with the dynamical time (black solid lines). Labels show log 10(t/ yr).

For the sake of completeness, in Fig. C.1 we show separately the diffusion time for ambipolar, Hall, and Ohmic diffusion computed with ζ^H₂ from model L₂. As expected, at the high densities reached in the disc region, ambipolar diffusion time is more than one order of magnitude larger than the dynamical time scale. On the contrary, Hall and Ohmic diffusion times are comparable and lower than the dynamical time scale in a region whose radius decreases from about 50 AU to 25 AU with decreasing minimum grain size. This is another way to say that gas-magnetic field decoupling is due to Hall and Ohmic diffusion at densities higher than n ~ 10¹⁰ cm^-3.

Acknowledgments

M.P. and P.H. acknowledge the financial support of the Agence National pour la Recherche (ANR) through the COSMIS project. M.P. and D.G. also acknowledge the support of the CNRS-INAF PICS project “Pulsar wind nebulae, supernova remnants and the origin of cosmic rays”. This work has been carried out thanks to the support of the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the “Investissements d’Avenir” French government programme managed by the ANR. B.C. acknowledges the financial support of the ANR through the “Feedback ISM” project.

References

All Tables

All Figures

	Fig. 1 CR ionisation maps and iso-density contours (black solid lines) for the case H in Table 1. Left panels: entire computational domain while right panels: a zoom in the inner region. Upper and lower panels: two perpendicular planes both containing the density peak. Labels show log 10 [ n/ cm-3 ].
In the text

	Fig. 2 Ionisation fractions for model L2 as a function of the volume density computed with a constant CR ionisation rate ζH2= 5 × 10-17 s-1 (constant-ζ, solid lines) and with a spatially resolved ζH2 (variable-ζ, dots). The three panels are for amin = 10-5 cm (top), amin = 10-6 cm (middle) and amin = 10-7 cm (bottom).
In the text

	Fig. 3 Ionisation fractions as a function of the CR ionisation rate at density n(H2) = 106 cm-3 assuming large grains (amin = 10-5 cm, solid lines) and small grains (amin = 10-7 cm, dotted lines). Dashed grey lines show the trend for xi ∝ ζH2 and xi ∝ (ζH2)1 / 2.
In the text

	Fig. 4 Density contours (solid blue lines with labels indicating log 10 [ n/ cm-3 ]) superposed to the velocity field (red arrows). The shaded areas show regions dominated by ambipolar (white), Hall (grey), and Ohmic (black) diffusion. The resistivities are calculated for three different values of amin and for the three models L1, L2, and H.
In the text

	Fig. 5 Diffusion time contours (white solid lines) evaluated with constant ζH2= 5 × 10-17 s-1 (left column) and ζH2 from model L1 in Table 1 (right column) compared with dynamical time contours (black dashed lines) for three different values of the minimum grain size. Labels show log 10(t/ yr).
In the text

	Fig. 6 Diffusion time contours (white solid lines) evaluated with constant ζH2= 5 × 10-17 s-1 (left column) and ζH2 from model L2 in Table 1 (right column) compared with dynamical time contours (black dashed lines) for three different values of the minimum grain size. Labels show log 10(t/ yr).
In the text

	Fig. 7 Diffusion time contours (white solid lines) evaluated with constant ζH2= 5 × 10-17 s-1 (left column) and ζH2 from model H in Table 1 (right column) compared with dynamical time contours (black dashed lines) for three different values of the minimum grain size. Labels show log 10(t/ yr).
In the text

Fig. B.1 Ambipolar (upper panel), Hall (middle panel), and Ohmic (lower panel) diffusion coefficients as a function of the mean grain radius computed at different molecular hydrogen densities: 105 (thick solid line), 106 (thick dashed line), 107 (thick dash-dotted line), 108 (dotted line), 109 (thin solid line), 1010 (thin dashed line), 1011 cm-3 (thin dash-dotted line), and 1012 cm-3 (thin dotted line).

In the text

	Fig. C.1 Ambipolar, Hall, and Ohmic contribution to the diffusion time (red dotted lines) compared with the dynamical time (black solid lines). Labels show log 10(t/ yr).
In the text