Open Access
Issue
A&A
Volume 666, October 2022
Article Number A27
Number of page(s) 7
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202142551
Published online 30 September 2022

© P. Marchand et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In Marchand et al. (2021, hereafter Paper I), we presented a novel method to calculate the coagulation and ionization of grains at a low computational cost. The only ionization source in that model was cosmic rays, which are dominant for isolated dense cores in the interstellar medium. At high temperatures, however, thermal ionization plays a major role in determining the ionization equilibrium of the gas-grain mixture, with a significant impact on the nonideal magnetohydrodynamics (MHD) resistivities.

The first stage of the protostellar collapse is isothermal at ∼10 K until density reaches ∼10−13 g cm−3. At this point, the dust-gas mixture becomes opaque to its own thermal radiation and the temperature rises as a core forms and contracts slowly (Larson 1969). At 2000 K, the dissociation of H2 molecules absorbs energy, allowing a rapid second collapse that leads to the creation of the protostar when all H2 is depleted. The thermal ionization of hydrogen occurs during the second collapse, and all hydrogen becomes ionized early in the protostar’s life, provoking a drop in resistivities to virtually zero (Marchand et al. 2016). Accounting for the thermal ionization of hydrogen is therefore critical to accurately describe the transition between the first core and the protostar, and thus from the nonideal to ideal MHD regime.

Another example of how this could be applicable pertains to chondrule formation. Chondrules are molten grains found in meteorites, whose formation requires rapid heating to ∼2000 K (Ebel et al. 2012). While there is no consensus on this topic, it has been proposed that their creation takes place in magnetic current sheets in protoplanetary disks, which may reach temperatures > 1500 K (Joung et al. 2004; McNally et al. 2014). Those high temperatures would trigger the thermal ionization of K and Na, leading to a sharp decrease in resistivities. Subsequently, the downward gradient in resistivities may create an instability that would allow the magnetic field to pile up in the current sheet (Hubbard et al. 2012). That phenomenon may be the origin of thunderclaps and extremely localized heating of the grains and gas (McNally et al. 2013), which are necessary conditions for chrondrule formation. However, readers should refer to Desch & Turner (2015), who argue that insufficient alkali metals would evaporate from grains to allow for this instability to act.

In this paper, we focus on extending the ionization model of Paper I by including the thermal ionization of one gas-phase species. In Sect. 2 we analytically derive the grain charge (Sect. 2.1) and thermal ionization equilibrium (Sect. 2.2), whose numerical implementations are described in Sect. 3. In Sect. 4 we discuss applications, and in Sect. 5 we present our conclusions.

2. Analytical method

2.1. Grain charge

Let us consider the two ionic species i and s of number density ni and ns. Species i corresponds to all the ions that are exclusively ionized by cosmic rays in the same manner as in Paper I, for which we assume an average atomic mass μi = mi/mH. Species s, however, undergoes both cosmic-ray and thermal ionization. We define ns, 0 as the total abundance of species s (both neutral and charged), so that ns, 0 is an upper bound for ns. We also consider an arbitrary size distribution of dust grains.

Grain charges Ze, with e being the electron charge, fluctuate stochastically due to the collection of electrons and the recombination of ions on their surfaces. The grain charge equilibrium (Eq. (4.3) of Draine & Sutin 1987, and Eq. (24) of Paper I) can be written as

f ( Z , τ k ) ( J i ( Z , τ k ) + J s ( Z , τ k ) ) = f ( Z + 1 , τ k ) J e ( Z + 1 , τ k ) , $$ \begin{aligned} f(Z,\tau _k)(J_i(Z,\tau _k)+J_s(Z,\tau _k))=f(Z+1,\tau _k)J_\mathrm{e} (Z+1,\tau _k), \end{aligned} $$(1)

where

τ k = a k k B T e 2 $$ \begin{aligned} \tau _k= \frac{a_k k_\mathrm{B} T}{e^2} \end{aligned} $$(2)

is the reduced temperature of a grain of radius ak and temperature T, while kB is the Boltzmann constant. Furthermore, f(Z, τk) is the distribution function of grain charges for a grain of reduced temperature τk, and Ji, Js, and Je are the fluxes of species i, s, and electrons onto the grains, respectively. The fluxes are (Draine & Sutin 1987)

J j ( Z , τ k ) = n j s j v j π a k 2 J ( Z e / q j , τ k ) , $$ \begin{aligned} J_j(Z,\tau _k) = n_j s_j v_j \pi a_k^2 \tilde{J}(Ze/q_j,\tau _k), \end{aligned} $$(3)

with nj being the abundance of species j (for j = i, s, or e), sj being the sticking probability of species j on grains, mj and qj being the mass and charge of species j, and vj = (8kBT/πmj)1/2 being the thermal speed of j. The polarization factor of grains for species j is J $ \tilde{J} $, which depends on the relative signs of Z and qj.

Similarly to the case without thermal ionization, presented in Paper I, we need to solve Eq. (1) for small grains and low temperatures, for which τk ≪ 1 and f(Z) is only significant for Z = −1, 0 and 1; that is, these grains hold a maximum of one charge. Therefore

f ( 1 , τ k ) + f ( 0 , τ k ) + f ( 1 , τ k ) = 1 . $$ \begin{aligned} f(-1,\tau _k)+f(0,\tau _k)+f(1,\tau _k)=1. \end{aligned} $$(4)

Equation (1) can be rewritten for Z = 0 and Z = −1 as follows:

f ( 1 , τ k ) = f ( 0 , τ k ) J e ( 0 , τ k ) J i ( 1 , τ k ) + J s ( 1 , τ k ) , $$ \begin{aligned} f(-1,\tau _k)&=f(0,\tau _k)\frac{J_\mathrm{e} (0,\tau _k)}{J_\mathrm{i} (-1,\tau _k)+J_\mathrm{s} (-1,\tau _k)}, \end{aligned} $$(5)

f ( 1 , τ k ) = f ( 0 , τ k ) J i ( 0 , τ k ) + J s ( 0 , τ k ) J e ( 1 , τ k ) . $$ \begin{aligned} f(1,\tau _k)&=f(0,\tau _k)\frac{J_\mathrm{i} (0,\tau _k)+J_\mathrm{s} (0,\tau _k)}{J_\mathrm{e} (1,\tau _k)}. \end{aligned} $$(6)

As in Paper I, we have

J ( 0 , τ k ) ( π / 2 τ k ) 1 2 , $$ \begin{aligned} \tilde{J}(0,\tau _k)&\approx \left(\pi / 2\tau _k\right)^{\frac{1}{2}},\end{aligned} $$(7)

J ( 1 , τ k ) ( 2 / τ k ) , $$ \begin{aligned} \tilde{J}(-1,\tau _k)&\approx (2 / \tau _k), \end{aligned} $$(8)

where J ( 0 , τ k ) $ \tilde{J}(0,\tau_k) $ appears in Je(0, τk), Ji(0, τk), and Js(0, τk), while J ( 1 , τ k ) $ \tilde{J}(-1,\tau_k) $ appears in Je(1, τk), Ji(−1, τk), and Js(−1, τk). We can then solve the equation system (4)–(6). With si = ss = 1, we obtain

f ( 0 , τ k ) = 1 1 + 1 α k [ ϵ Θ + 1 ϵ Θ ] , $$ \begin{aligned} f(0,\tau _k)&= \frac{1}{ 1 + \frac{1}{\alpha _k} \left[\epsilon \Theta +\frac{1}{\epsilon \Theta } \right]},\end{aligned} $$(9)

f ( 1 , τ k ) = ϵ 2 Θ 2 1 + α k ϵ Θ + ϵ 2 Θ 2 , $$ \begin{aligned} f(-1,\tau _k)&= \frac{\epsilon ^2 \Theta ^2}{1 + \alpha _k \epsilon \Theta + \epsilon ^2 \Theta ^2}, \end{aligned} $$(10)

f ( 1 , τ k ) = 1 1 + α k ϵ Θ + ϵ 2 Θ 2 , $$ \begin{aligned} f(1,\tau _k)&= \frac{1}{1 + \alpha _k \epsilon \Theta + \epsilon ^2 \Theta ^2,} \end{aligned} $$(11)

where Θ = s e ( μ i m H / m e ) 1 2 $ \Theta=s_{\mathrm{e}} ({{{\upmu}}_\text{i}}{m_\text{H}}/ {m_\text{e}})^{\frac12} $, q is = ( μ i / μ s ) 1 2 $ {q_\text{is}}= ({{{\upmu}}_\text{i}}/ {{{\upmu}}_\text{s}})^\frac12 $, and

ϵ = n e n i + q is n s · $$ \begin{aligned} \epsilon =\frac{n_\mathrm{e} }{n_\mathrm{i} +q_\mathrm{is} n_\mathrm{s} }\cdot \end{aligned} $$(12)

The main difference with Paper I is the appearance of the term qis. While ni + ns is the total abundance of ions, ni + qisns is an effective abundance that reflects the relative flux of ions onto grains. The average charge of grains

Z k = Z = 1 Z = 1 Z f ( Z , τ k ) $$ \begin{aligned} Z_k = \sum _{Z=-1}^{Z=1} Zf(Z,\tau _k) \end{aligned} $$(13)

and the grain-ion recombination enhancement factor

J ( τ k ) = Z = 1 Z = 1 J ( Z e / q i , τ k ) f ( Z , τ k ) $$ \begin{aligned} \langle \tilde{J}(\tau _k) \rangle = \sum _{Z=-1}^{Z=1} \tilde{J}(Ze/q_\mathrm{i} ,\tau _k)f(Z,\tau _k) \end{aligned} $$(14)

thus are given by the same expressions as in Paper I,

Z k = 1 ϵ 2 Θ 2 1 + α k ϵ Θ + ϵ 2 Θ 2 , $$ \begin{aligned} Z_k&= \frac{1-\epsilon ^2 \Theta ^2}{1+\alpha _k \epsilon \Theta +\epsilon ^2 \Theta ^2},\end{aligned} $$(15)

J ( τ k ) = 2 τ k ( ϵ 2 Θ 2 + ϵ Θ ) ϵ 2 Θ 2 + α k ϵ Θ + 1 , $$ \begin{aligned} \langle \tilde{J}(\tau _k) \rangle&= \frac{\frac{2}{\tau _k} (\epsilon ^2 \Theta ^2 + \epsilon \Theta )}{\epsilon ^2 \Theta ^2 + \alpha _k \epsilon \Theta +1}, \end{aligned} $$(16)

where we neglected the recombination of ions on positively charged grains ( J ( 1 , τ k ) $ \tilde{J}(1,\tau_k) $).

For the larger grains (τk ≫ 1), the same kind of change needs to be made to the Spitzer equation (Spitzer 1949; Draine & Sutin 1987) that governs the grain’s electric potential ψ. For ψ < 0, eψ represents the repulsion of the flux of electrons by the negatively charged grains, while 1 − ψ characterizes the attraction of the flux of ions. The flux equilibrium can thus be written as

( n e v e s e ) e ψ = ( n i v i s i + n s v s s s ) ( 1 ψ ) . $$ \begin{aligned} (n_\mathrm{e} v_\mathrm{e} s_\mathrm{e} )e^\psi = (n_\mathrm{i} v_\mathrm{i} s_\mathrm{i} + n_\mathrm{s} v_\mathrm{s} s_\mathrm{s} )(1-\psi ). \end{aligned} $$(17)

Introducing the same notations as above, we can write

ϵ = 1 ψ Θ e ψ , $$ \begin{aligned} \epsilon = \frac{1-\psi }{\Theta e^\psi }, \end{aligned} $$(18)

which is the same equation as the one-ion model of Paper I (Eq. (34)) with the modified expression for ϵ. The average charge of large grains and the grain-ion recombination enhancement factor yield the same expressions as in Paper I,

Z k = ψ τ k , $$ \begin{aligned} Z_k&= \psi \tau _k,\end{aligned} $$(19)

J ( τ k ) = ( 1 ψ ) . $$ \begin{aligned} \langle \tilde{J}(\tau _k) \rangle&= (1-\psi ). \end{aligned} $$(20)

The average charge and recombination enhancement factor for a mix of small and large grains is assumed to be the sum of the contributions from both Eqs. (15),(19) and Eqs. (16),(20) (Draine & Sutin 1987).

2.2. Ionization equilibrium

We always assume charge neutrality,

n i + n s n e + n k Z k = 0 . $$ \begin{aligned} n_\mathrm{i} +n_\mathrm{s} - n_\mathrm{e} + \sum n_k Z_k = 0. \end{aligned} $$(21)

This allows us to find the ionization equilibrium for species i. In Paper I, we considered the balance between the creation of species i by cosmic-ray ionization, and the destruction of species i by recombination with electrons and with grains. In this two-ion model, we also need to consider the charge exchange reactions between species i and s. The one-ion model hides and summarizes all the charge transfer reactions between gas-phase species in the choice of μi. Here, we have to explicitly account for the creation of species i by the destruction of species s, and vice versa. The ionization equilibrium is then

ζ ( n H n s , 0 n i ) + k s , i ( n H n s , 0 n i ) n s = σ v ie n e n i + n i v i n k π a k 2 J k + k i , s ( n s , 0 n s ) n i , $$ \begin{aligned}&\zeta (n_\mathrm{H} -n_{\mathrm{s} ,0}-n_\mathrm{i} ) + k_\mathrm{s,i} (n_\mathrm{H} -n_{\mathrm{s} ,0}-n_\mathrm{i} )n_\mathrm{s} \nonumber \\&= \langle \sigma v \rangle _\mathrm{ie} n_\mathrm{e} n_\mathrm{i} + n_\mathrm{i} v_\mathrm{i} \sum n_k \pi a_k^2 J_k + k_\mathrm{i,s} (n_{\mathrm{s} ,0}-n_\mathrm{s} ) n_\mathrm{i} , \end{aligned} $$(22)

where ζ is the cosmic-ray ionization rate, ks,i and ki,s are the chemical reaction rates of species s → i and i → s, respectively, and we consider the recombination rate of ions i with electrons to be σ v ie = 2 × 10 7 ( T / 300 ) 1 2 $ {\langle\sigma v \rangle_\text{ie}}= 2\times 10^{-7} (T/300)^{\frac{1}{2}} $ cm3 s−1, based on the recombination rate of HCO+ taken from the UMIST database (McElroy et al. 2013). The terms of the form ka, bnanc are the transformation of species c to species b, through the chemical reaction with species a. Hence the term (nH − ns, 0 − ni) represents the total abundance of neutral species that can be ionized into ion i, as (ns, 0 − ns) is the total abundance of neutral species s that can be ionized to ion s.

The ionization equilibrium for species s is similar, with the addition of a thermal ionization term (Pneuman & Mitchell 1965, also see Sect. 4)

( d n s d t ) à ¾ thermal à ¾ ionization = β s n H T 1 / 2 e T 0 , s / T ( n s , 0 n s ) , $$ \begin{aligned} \left(\frac{\mathrm{d}n_\mathrm{s} }{\mathrm{d}t}\right)_{ \begin{array}{l} þ{\scriptstyle \mathrm{thermal}}\\ þ{\scriptstyle \mathrm{ionization}} \end{array}} = \beta _\mathrm{s} n_\mathrm{H} T^{1/2} {e}^{-T_{0,\mathrm{s} }/T} (n_{\mathrm{s} ,0}- n_\mathrm{s} ), \end{aligned} $$(23)

with the values of β and T0 depending on the species. Table 1 summarizes the values of constants specific to species s for the cases of sodium, potassium, and hydrogen. The ionization equilibrium equation is then

ζ ( n s , 0 n s ) + ( d n s d t ) à ¾ thermal à ¾ ionization + k i , s ( n s , 0 n s ) n i = σ v se n e n s + n s v s n k π a k 2 J k + k s , i ( n H n s , 0 n i ) n s . $$ \begin{aligned}&\zeta (n_{\mathrm{s} ,0}-n_\mathrm{s} ) + \left(\frac{\mathrm{d}n_\mathrm{s} }{\mathrm{d}t}\right)_{\begin{array}{l} þ{\scriptstyle \mathrm{thermal}}\\ þ{\scriptstyle \mathrm{ionization}} \end{array}} + k_\mathrm{i,s} (n_{\mathrm{s} ,0}-n_\mathrm{s} ) n_\mathrm{i} \nonumber \\&= \langle \sigma v \rangle _\mathrm{se} n_\mathrm{e} n_\mathrm{s} + n_\mathrm{s} v_\mathrm{s} \sum n_k \pi a_k^2 J_{k}+ k_\mathrm{s,i} (n_\mathrm{H} -n_{\mathrm{s} ,0}-n_\mathrm{i} )n_\mathrm{s} . \end{aligned} $$(24)

Table 1.

Constants depending on the choice of species s for sodium, potassium, and hydrogen.

Those equations are valid if the Saha equation is valid as well, meaning that there should be a large number of particles within a Debye length of each other. The validity condition is then

4 3 π n H λ D 3 1 , $$ \begin{aligned} \frac{4}{3}\pi n_\mathrm{H} \lambda _\mathrm{D} ^3 \gg 1, \end{aligned} $$(25)

with

λ D = k B T 4 π n i e 2 · $$ \begin{aligned} \lambda _\mathrm{D} =\sqrt{\frac{k_\mathrm{B} T}{4\pi n_\mathrm{i} e^2}}\cdot \end{aligned} $$(26)

3. Numerical implementation and tests

Equations (18), (21), (22), and (24) need to be solved for ψ, ϵ, ni, and ns. In this section, we discuss the solution for this system with four equations and four unknowns. The system could be reduced to three equations, as Eq. (18) is an explicit expression of ψ as a function of ϵ. That would, however, significantly increase the analytical and numerical complexity of the calculation, and it is unclear whether this would lead to better performances or not.

3.1. Numerical convergence

Although the system of equations is valid for a wide range of physically valid parameters, we need to be cautious to ensure numerical convergence toward the solution, especially at high density and temperature. At high density, ψ converges toward zero by a negative value and becomes very small in an absolute value. At high temperature, ns overwhelmingly dominates ni due to the thermal ionization. Therefore, the numerical implementation has to be robust for the cases |ψ|≪1 and ns/ni ≫ 1.

For this purpose, the four equations must be normalized so that they can be written in the form 1 + x = 0 to avoid sums of very large or very small numbers. It is therefore necessary to include species s in the normalization of the equations to avoid convergence issues at large temperatures when ns grows much larger than ni. We therefore normalized Eq. (21) by the total number of ions ni + ns, we used both the cosmic-ray ionization rate and the chemical reaction rate s → i to normalize Eq. (22), and we included the thermal ionization term in the normalization of Eq. (24).

Another issue arises from the average grain charge (15). Before the thermal ionization starts to be relevant and ni ≫ ns, ϵ converges toward 1/Θ as density increases (see Fig. 2 of Paper I). When the difference between ϵ and 1/Θ becomes close to machine precision, the term 1 − (ϵΘ)2 reaches a lower bound1 which prevents the convergence of the charge neutrality (Eq. (21)), as the grain charge fails to decrease. A solution to avoid this issue is to replace ϵ by ϵ′+ϵ0, with ϵ0 = 1/Θ. This substitution has to be made in all the equations. The new variable to find is ϵ′ which converges toward zero instead of 1/Θ. The term 1 − (ϵΘ)2 in Eq. (15) is then mathematically equal to and can be replaced by −2ϵ′Θ − (ϵ′Θ)2, which avoids the lower-bound issue.

3.2. The normalized system of equations

The numerical solution of this system requires its normalization, as discussed in Sect. 3.1. We define the function F(ψ, ϵ′,ni, ns) = (f1, f2, f3, f4), with

f 1 = 1 ψ ( ϵ + ϵ 0 ) Θ e ψ 1 = 0 , $$ \begin{aligned}&f_1 = \frac{1-\psi }{(\epsilon^\prime +\epsilon _0) \Theta e^\psi } -1 = 0,\end{aligned} $$(27)

f 2 = 1 ( ϵ + ϵ 0 ) n i + q is n s n i + n s + 1 n i + n s n k Z k = 0 , $$ \begin{aligned}&f_2 = 1 - (\epsilon^\prime +\epsilon _0)\frac{n_\mathrm{i} +q_\mathrm{is} n_\mathrm{s} }{n_\mathrm{i} +n_\mathrm{s} } +\frac{1}{n_\mathrm{i} +n_\mathrm{s} } \sum n_k Z_k = 0,\end{aligned} $$(28)

f 3 = 1 n s , 0 + n i n H σ v ie ( ϵ + ϵ 0 ) n i ( n i + q is n s ) ( ζ + k s , i n s ) n H n i v i n k π a k 2 J ( τ k ) ( ζ + k s , i n s ) n H k i , s n i ( n s , 0 n s ) ( ζ + k s , i n s ) n H = 0 , $$ \begin{aligned}&f_3 = 1- \frac{n_{\mathrm{s} ,0}+n_\mathrm{i} }{n_\mathrm{H} }- \frac{\langle \sigma v \rangle _\mathrm{ie} (\epsilon^\prime +\epsilon _0) n_\mathrm{i} (n_\mathrm{i} +q_\mathrm{is} n_\mathrm{s} )}{(\zeta +k_\mathrm{s,i} n_\mathrm{s} ) n_\mathrm{H} } \nonumber \\&\quad \quad - \frac{n_\mathrm{i} v_\mathrm{i} \sum n_k \pi a_k^2 \langle \tilde{J}(\tau _k) \rangle }{(\zeta +k_\mathrm{s,i} n_\mathrm{s} )n_\mathrm{H} } -\frac{k_\mathrm{i,s} n_\mathrm{i} (n_{\mathrm{s} ,0}-n_\mathrm{s} )}{(\zeta +k_\mathrm{s,i} n_\mathrm{s} )n_\mathrm{H} } = 0,\end{aligned} $$(29)

f 4 = ( 1 n s n s , 0 ) ( 1 + k i , s n i ζ + β T 1 2 e T 0 T n H ) σ v se ( ϵ + ϵ 0 ) n s ( n i + q is n s ) ( ζ + β T 1 2 e T 0 T n H ) n s , 0 n s v s n k π a k 2 J ( τ k ) ( ζ + β T 1 2 e T 0 T n H ) n s , 0 k s , i n s ( n H n s , 0 n i ) ( ζ + β T 1 2 e T 0 T n H ) n s , 0 = 0 , $$ \begin{aligned}&f_4 = \left(1-\frac{n_\mathrm{s} }{n_{\mathrm{s} ,0}}\right)\left(1 + \frac{k_\mathrm{i,s} n_\mathrm{i} }{\zeta +\beta T^\frac{1}{2}e^{-\frac{T_0}{T}}n_\mathrm{H} }\right) \nonumber \\&\quad \quad -\frac{\langle \sigma v \rangle _\mathrm{se} (\epsilon^\prime +\epsilon _0) n_\mathrm{s} (n_\mathrm{i} +q_\mathrm{is} n_\mathrm{s} )}{(\zeta +\beta T^\frac{1}{2}e^{-\frac{T_0}{T}}n_\mathrm{H} )n_{\mathrm{s} ,0}} \nonumber \\&\quad \quad - \frac{n_\mathrm{s} v_\mathrm{s} \sum n_k \pi a_k^2 \langle \tilde{J}(\tau _k) \rangle }{(\zeta +\beta T^\frac{1}{2}e^{-\frac{T_0}{T}}n_\mathrm{H} )n_{\mathrm{s} ,0}} \nonumber \\&\quad \quad - \frac{k_\mathrm{s,i} n_\mathrm{s} (n_\mathrm{H} - n_{\mathrm{s} ,0}-n_\mathrm{i} )}{(\zeta +\beta T^\frac{1}{2}e^{-\frac{T_0}{T}}n_\mathrm{H} )n_{\mathrm{s} ,0}}= 0, \end{aligned} $$(30)

where

Z k = ψ τ k + 2 ϵ Θ ( ϵ Θ ) 2 1 + α k ( ϵ + ϵ 0 ) Θ + ( ϵ + ϵ 0 ) 2 Θ 2 , $$ \begin{aligned} Z_k&= \psi \tau _k + \frac{-2\epsilon^\prime \Theta -(\epsilon^\prime \Theta )^2}{1+\alpha _k (\epsilon^\prime +\epsilon _0) \Theta +(\epsilon^\prime +\epsilon _0)^2 \Theta ^2},\end{aligned} $$(31)

J ( τ k ) = ( 1 ψ ) + 2 τ k ( ( ϵ + ϵ 0 ) 2 Θ 2 + ( ϵ + ϵ 0 ) Θ ) ( ϵ + ϵ 0 ) 2 Θ 2 + α k ( ϵ + ϵ 0 ) Θ + 1 , $$ \begin{aligned} \langle \tilde{J}(\tau _k) \rangle&= (1-\psi ) + \frac{\frac{2}{\tau _k} ((\epsilon^\prime +\epsilon _0)^2 \Theta ^2 + (\epsilon^\prime +\epsilon _0) \Theta )}{(\epsilon^\prime +\epsilon _0)^2 \Theta ^2 + \alpha _k (\epsilon^\prime +\epsilon _0) \Theta +1},\end{aligned} $$(32)

ϵ = n e n s = ϵ + ϵ 0 , $$ \begin{aligned} \epsilon&= \frac{n_\mathrm{e} }{n_\mathrm{s} }=\epsilon^\prime + \epsilon _0,\end{aligned} $$(33)

ϵ 0 = 1 Θ , $$ \begin{aligned} \epsilon _0&= \frac{1}{\Theta }, \end{aligned} $$(34)

and

v i , s = ( 8 k B T π μ i , s m H ) 1 / 2 . $$ \begin{aligned} v_\mathrm{i,s} = \left(\frac{8k_\mathrm{B} T}{\pi \upmu _\mathrm{i,s} m_\mathrm{H} }\right)^{1/2}. \end{aligned} $$(35)

Here, vi, s and μi, s stand for either vi and μi or vs and μs.

3.3. Solving the system

Similarly to Paper I, we used a Newton–Raphson method to solve the equation system. Let X = (ψ, ϵ, ni, ns). Starting from an educated guess X0, we iterated

X n + 1 = X n J ( X n ) 1 F ( X n ) , $$ \begin{aligned} \mathbf X _{n+1} = \mathbf X _n - \mathbb{J} (\mathbf X _n)^{-1} \mathbf F (\mathbf X _n), \end{aligned} $$(36)

until ||F(Xn)||< δ ≪ 1. The matrix 𝕁 is the Jacobian of the system defined by 𝕁i, j = ∂fi/∂Xj. The full analytic components of the Jacobian matrix are given in Appendix A.

For reliable convergence, this iterative solution for the system of equations is best started from as close an estimate as possible. For low density (nH < 107 cm−3) and low temperature (typically 10 K), a good starting point is

ψ = ψ 0 , $$ \begin{aligned} \psi&= \psi _0,\end{aligned} $$(37)

ϵ = 0.9999 , $$ \begin{aligned} \epsilon&= 0.9999,\end{aligned} $$(38)

n i = 10 7 n H , $$ \begin{aligned} n_\mathrm{i}&= 10^{-7} n_\mathrm{H} ,\end{aligned} $$(39)

n s = 10 4 n i , $$ \begin{aligned} n_\mathrm{s}&= 10^{-4} n_\mathrm{i} , \end{aligned} $$(40)

where ψ0 is the solution of

1 ψ Θ e ψ = 1 . $$ \begin{aligned} \frac{1-\psi }{\Theta e^{\psi }}=1. \end{aligned} $$(41)

For larger densities and temperatures, we recommend solving the system for a gradual increase in those quantities, using the previous solution as a first estimate. In particular, ns increases very quickly with density and temperature once the thermal ionization starts, and large leaps may lead to convergence failure.

3.4. Tests

We solved the normalized system of Eqs. (27)–(30). For testing purposes, we used typical parameters of star-forming environments to compare with the existing literature, but the reader should keep in mind that a wide range of physically sound parameters is possible. The density spans nH = 104–1025 cm−3, starting from the lowest density and increasing nH gradually. We assumed the same barotropic equation of state as in Marchand et al. (2016) to emulate the rise in temperature during a protostellar collapse

T = T 0 ( 1 + [ n H n 1 ] 0.8 ) 1 2 ( 1 + [ n H n 2 ] ) 0.3 ( 1 + [ n H n 3 ] ) 1.7 3 , $$ \begin{aligned} T = T_0 \left(1+\left[\frac{n_\mathrm{H} }{n_1}\right]^{0.8}\right)^{\frac{1}{2}}\left(1+\left[\frac{n_\mathrm{H} }{n_2}\right]\right)^{-0.3}\left(1+\left[\frac{n_\mathrm{H} }{n_3}\right]\right)^{-\frac{1.7}{3}}, \end{aligned} $$(42)

with T0 = 10 K, n1 = 1011 cm−3, n2 = 1016 cm−3, and n3 = 1021 cm−3. We assumed a nonevolving Mathis, Rumpl, Nordsieck (MRN) grain-size distribution (Mathis et al. 1977), with a slope of −3.5 between minimum grain size amin = 5 nm and maximum grain size amax = 250 nm, sampled by 26 bins. The grain bulk density is ρ = 2.9 g cm−3 and the dust-to-gas mass ratio is 1%, which is typical of the insterstellar medium (Bohlin et al. 1978). We set ζ = 5 × 10−17 s−1 (Padovani et al. 2013), se = 0.5 (Umebayashi & Nakano 1990), and μi = 25 (close to the molecular mass of Mg, Fe, or HCO+, Marchand et al. 2016).

In this test, the temperature exceeds several 103 K, above which all grains should be quickly destroyed by evaporation or sputtering (Lenzuni et al. 1995). This is, however, not an issue since at such high temperatures, the contribution of the (computed) charge of grains is negligible compared to that from ions and electrons. It is also possible to combine our method with any grain destruction model. Marchand et al. (2016) assumed that the grain evaporation would occur between ∼800 K and ∼1600 K in the density range nH ≈ 1016 − 1019 cm−3, which coincides with the beginning of the thermal ionization of K and Na. The species s is assumed to be K, Na, or H, using the values of Table 1. Desch & Turner (2015) show that K and Na may originally be confined to grains and have to be evaporated before being available for thermal ionization. Our method is compatible with models accounting for that process since ns, 0 can be freely modified at any time. For simplicity, however, we assume here that K and Na are already present in the gas phase in quantities given by ns, 0 in the table. Here, we present the evolution of ni/nH, ns/nH, ne/nH and the average grain charge for several bins. The results are displayed in Figs. 13, respectively.

thumbnail Fig. 1.

Test of our method for protostellar collapse conditions. Evolution of the average charge of grains (left axis) for several grain sizes (color lines), and the fractional abundance (right axis) of cosmic-ray ionized ions ni/nH (solid line), thermally ionized K+ ions ns/nH (dotted line), and electrons ne/nH (dashed line).

thumbnail Fig. 2.

Same as Fig. 1, but for Na. The only ions present are the cosmic-ray ionized ones i and Na (ion s), which is both cosmic-ray and thermally ionized.

thumbnail Fig. 3.

Same as Fig. 1, but for H. The only ions presents are the cosmic-ray ionized ones i and H (ions s), which is both cosmic-ray and thermally ionized.

These figures can be compared to Fig. 7 of Marchand et al. (2016). Although the abundances of Na+ and K+ seem overestimated compared to ne, and that of H+ seems underestimated (with a negligible impact on the resistivities), our model reproduces the evolution of abundances of this more detailed calculation at the key points fairly well: K and Na start their ionization around nH = 1018 cm−3 and T = 1600 K, and they saturate around nH = 1022 cm−3 and T = 1.5 × 104 K at the maximum fractional abundance of their respective species. Furethermore, H has a similar behavior, starting its thermal ionization at nH = 1020 cm−3 and T = 2650 K, turning virtually all neutrals into ions by nH = 1023 cm−3 and T = 7 × 104 K. The main difference with Fig. 7 of Marchand et al. (2016) is the earlier rise of electron density at nH ≈ 1016 cm−3 due to the thermionic emission of grains included in the complete chemical calculation. This emission is associated with a drop in neutral grain density and rise in positively and negatively charged grains.

In all three cases, the large input of electrons into the gas significantly increases the grain charges as well. We note that at the hydrogen ionization fraction nH+/nH > 0.1 reached at nH > 1023 cm−3, the Saha equation is no longer valid, so our model becomes imprecise. However, the resistivities are so low in this regime that ideal MHD is a valid approximation; thus this is never an issue in practice.

Figures 46 display the associated nonideal MHD resistivities using the formulae provided by Marchand et al. (2016). For display purposes only, we prescribe a magnetic field (Li et al. 2011)

B = 1.43 × 10 7 G n H 1 / 2 , $$ \begin{aligned} B=1.43\times 10^{-7} \mathrm{G} \ n_\mathrm{H} ^{1/2}, \end{aligned} $$(43)

thumbnail Fig. 4.

Evolution of the Ohmic (yellow), Hall (blue), and ambipolar (dark red) resistivities for the thermal ionization of K. The dashed line represents the Hall resistivity in negative values. The thin dotted lines are the resistivities in the absence of thermal ionization.

thumbnail Fig. 5.

Same as Fig. 4, but for Na.

thumbnail Fig. 6.

Same as Fig. 4, but for H.

which corresponds to the critical magnetic field strength of a spherical cloud. The magnetic field strength and the resistivities are therefore overestimated compared to protostellar collapse simulations with nonideal MHD for nH ≳ 1012 cm−3.

The nonideal MHD terms significantly affect the MHD evolution of protostellar collapse and protoplanetary disks for resistivities larger than ≈1018 cm2 s−1. Nonideal MHD terms are then important at all densities before thermal ionization starts (except the Ohmic diffusion at a low density). At this point, the abundance of charged species in the gas significantly increases, leading to a sharp decrease in all resistivities. The overall behavior is consistent with previous works (Kunz & Mouschovias 2010; Marchand et al. 2016; Wurster et al. 2016; Koga et al. 2019).

4. Discussion

The method presented here is applicable in a wide variety of environments, and it is particularly suited for modeling protostar formation and protoplanetary disks. At later stages of the star formation process, or in the presence of nearby massive stars, photoionization by UV or X-rays may become relevant (Getman & Feigelson 2021). In this case, their ionization rate can simply be added to the cosmic-ray ionization rate ζ in the system of Eqs. (27)–(30).

We describe our method to calculate the resistivities as fast in comparison to solving a full chemical network. We have implemented the algorithm in the 3D MHD RAMSES code (Teyssier 2002). The code previously calculated the resistivities by interpolating on the precalculated table of Marchand et al. (2016). Without the Hall effect, our thermal ionization algorithm is faster than reading the chemical table, as the calculation needs to be performed only once per cell per time-step. This is different with the Hall effect, which requires, in addition, the self-consistent calculation of resistivities on cell edges (Marchand et al. 2018). In this case, the code runs at similar speeds for both methods (it is important to note that this may vary with different implementations). However, the method presented in this paper is much more flexible than a precalculated table because the physical conditions, the chemical composition, and the grain size-distribution can be changed at any point for a self-consistent calculation.

In Table 1, we provide the thermal ionization coefficients for K, Na, and H from Pneuman & Mitchell (1965), which is a theoretical work, to match the rates used in Marchand et al. (2016). In the 1960s and 1970s, there were many discussions about the ionization rates of alkali metals. Flame experiments to measure those rates (Hollander et al. 1963; Ashton & Hayhurst 1973) resulted in much larger cross sections than theoretically predicted (Aller 1961; Hollenbach & Salpeter 1969, see Schofield 1965; Shui 1977 for review). Desch & Turner (2015) argue that the experimental value of Ashton & Hayhurst (1973) should be preferred to the theoretical value of Pneuman & Mitchell (1965). This is debatable as the conditions in flames may be difficult to control and different from astrophysical plasmas (pressure, chemical composition), where Na and K ionize by colliding with H2. The larger coefficient rates would suggest thermal ionization of K and Na at a lower temperature, typically T = 1200 K instead of T = 1600 K, which could be of importance in protoplanetary disks. For reference, we provide the experimental values of Ashton & Hayhurst (1973) for K and Na:

β K = 1.0 × 10 8 cm 3 s 1 K 1 2 , $$ \begin{aligned}&\beta _\mathrm{K} =1.0 \times 10^{-8}\,\mathrm{cm} ^3\,\mathrm{s} ^{-1}\,\mathrm{K} ^{-\frac{1}{2}},\end{aligned} $$(44)

T 0 , K = 5.6 × 10 4 K , $$ \begin{aligned}&T_{0,\mathrm{K} } =5.6 \times 10^4\,\mathrm{K} ,\end{aligned} $$(45)

β Na = 2.0 × 10 9 cm 3 s 1 K 1 2 , $$ \begin{aligned}&\beta _\mathrm{Na} =2.0 \times 10^{-9}\, \mathrm{cm} ^3\,\mathrm{s} ^{-1}\,\mathrm{K} ^{-\frac{1}{2}},\end{aligned} $$(46)

T 0 , Na = 5.6 × 10 4 K . $$ \begin{aligned}&T_{0,\mathrm{Na} } =5.6 \times 10^4\,\mathrm{K} . \end{aligned} $$(47)

5. Conclusions

We detail an extension of the ionization model of Marchand et al. (2021) to include the thermal ionization of one gas species. It is possible to increase the number of ionized species by adding their contribution in the same manner, at the price of a higher numerical cost. We have presented the examples of K, Na, and H. Both the chemical abundances and the resistivities show a behavior consistent with previous work. This method is then a powerful tool to self-consistently calculate the ionization of the dust-grain mixture and the nonideal MHD resistivities in hydrodynamical simulations, faster than a complete chemical network. The method is flexible and valid in a wide variety of environments, including star formation and protoplanetary disks.


1

For example, for a machine precision of 10−16, if |ϵ − 1/Θ|< 10−10, then 1 − (ϵΘ)2 > 10−6.

Acknowledgments

We thank the referee for their insightful comments that helped improve the manuscript. We thank Jérémy E. Cohen for his insight on the solving of linear systems. P. M. acknowledges financial support by the Kathryn W. Davis Postdoctoral Fellowship of the American Museum of Natural History. U. L. acknowledges financial support from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130). M.-M. M. L. acknowledges partial support from NSF grant AST18-15461.

References

  1. Aller, L. H. 1961, The Abundances of the Elements (New York: Interscience Publ) [Google Scholar]
  2. Ashton, A., & Hayhurst, A. 1973, Combust. Flame, 21, 69 [CrossRef] [Google Scholar]
  3. Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132 [Google Scholar]
  4. Desch, S. J., & Turner, N. J. 2015, ApJ, 811, 156 [NASA ADS] [CrossRef] [Google Scholar]
  5. Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803 [NASA ADS] [CrossRef] [Google Scholar]
  6. Ebel, D. S., Hubbard, A., McNally, C., et al. 2012, Meteorit. Planetary Sci. Suppl., 75, 5387 [NASA ADS] [Google Scholar]
  7. Getman, K. V., & Feigelson, E. D. 2021, ApJ, 916, 32 [NASA ADS] [CrossRef] [Google Scholar]
  8. Hollander, T., Kalff, P. J., & Alkemade, C. T. J. 1963, J. Chem. Phys., 39, 2558 [NASA ADS] [CrossRef] [Google Scholar]
  9. Hollenbach, D. J., & Salpeter, E. E. 1969, J. Chem. Phys., 50, 4157 [NASA ADS] [CrossRef] [Google Scholar]
  10. Hubbard, A., McNally, C. P., & Mac Low, M.-M. 2012, ApJ, 761, 58 [NASA ADS] [CrossRef] [Google Scholar]
  11. Joung, M. K. R., Mac Low, M.-M., & Ebel, D. S. 2004, ApJ, 606, 532 [NASA ADS] [CrossRef] [Google Scholar]
  12. Koga, S., Tsukamoto, Y., Okuzumi, S., & Machida, M. N. 2019, MNRAS, 484, 2119 [Google Scholar]
  13. Kunz, M. W., & Mouschovias, T. C. 2010, MNRAS, 408, 322 [NASA ADS] [CrossRef] [Google Scholar]
  14. Larson, R. B. 1969, MNRAS, 145, 271 [Google Scholar]
  15. Lenzuni, P., Gail, H.-P., & Henning, T. 1995, ApJ, 447, 848 [NASA ADS] [CrossRef] [Google Scholar]
  16. Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180 [Google Scholar]
  17. Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Marchand, P., Commerçon, B., & Chabrier, G. 2018, A&A, 619, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Marchand, P., Guillet, V., Lebreuilly, U., & Mac Low, M. M. 2021, A&A, 649, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [Google Scholar]
  21. McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. McNally, C. P., Hubbard, A., Mac Low, M.-M., Ebel, D. S., & D’Alessio, P. 2013, ApJ, 767, L2 [NASA ADS] [CrossRef] [Google Scholar]
  23. McNally, C. P., Hubbard, A., Yang, C.-C., & Mac Low, M.-M. 2014, ApJ, 791, 62 [NASA ADS] [CrossRef] [Google Scholar]
  24. Padovani, M., Hennebelle, P., & Galli, D. 2013, A&A, 560, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Pneuman, G. W., & Mitchell, T. P. 1965, Icarus, 4, 494 [NASA ADS] [CrossRef] [Google Scholar]
  26. Schofield, K., & Sugden, S. T. 1965, Some Observations on the Ionization of Alkali and Alkaline-earth Elements in Hydrogen Flames (Cambridge: University of Cambridge), 10 [Google Scholar]
  27. Shui, V. H. 1977, Phys. Fluids, 20, 32 [NASA ADS] [CrossRef] [Google Scholar]
  28. Spitzer, L., Jr 1949, Leaflet Astron. Soc. Pac., 5, 336 [Google Scholar]
  29. Teyssier, R. 2002, A&A, 385, 337 [CrossRef] [EDP Sciences] [Google Scholar]
  30. Umebayashi, T., & Nakano, T. 1990, MNRAS, 243, 103 [NASA ADS] [CrossRef] [Google Scholar]
  31. Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037 [Google Scholar]

Appendix A: Jacobian components

The components of the normalized Jacobian matrix required for the solution of equation (36) are

J 1 , 1 = f 1 ψ = ψ 2 ( ϵ + ϵ 0 ) Θ e ψ , $$ \begin{aligned} \mathbb{J} _{1,1} = \frac{\partial f_1}{\partial \psi } =&\frac{\psi -2}{(\epsilon^\prime +\epsilon _0) \Theta e^{\psi }},\end{aligned} $$(A.1)

J 1 , 2 = f 1 ϵ = 1 ψ ( ϵ + ϵ 0 ) 2 Θ e ψ , $$ \begin{aligned} \mathbb{J} _{1,2} = \frac{\partial f_1}{\partial \epsilon^\prime } =&-\frac{1-\psi }{(\epsilon^\prime +\epsilon _0)^2 \Theta e^{\psi }},\end{aligned} $$(A.2)

J 1 , 3 = f 1 n i = 0 , $$ \begin{aligned} \mathbb{J} _{1,3} = \frac{\partial f_1}{\partial n_\mathrm{i} } =&0,\end{aligned} $$(A.3)

J 1 , 4 = f 1 n s = 0 , $$ \begin{aligned} \mathbb{J} _{1,4} = \frac{\partial f_1}{\partial n_\mathrm{s} } =&0, \end{aligned} $$(A.4)

J 2 , 1 = f 2 ψ = 1 n i + n s n k τ k , $$ \begin{aligned} \mathbb{J} _{2,1} = \frac{\partial f_2}{\partial \psi } =&\frac{1}{n_\mathrm{i} +n_\mathrm{s} }\sum n_k \tau _k,\end{aligned} $$(A.5)

J 2 , 2 = f 2 ϵ = n i + q is n s n i + n s + 1 n i + n s n k Z k ϵ , $$ \begin{aligned} \mathbb{J} _{2,2} = \frac{\partial f_2}{\partial \epsilon^\prime } =&-\frac{n_\mathrm{i} +q_\mathrm{is} n_\mathrm{s} }{n_\mathrm{i} +n_\mathrm{s} }+\frac{1}{n_\mathrm{i} +n_\mathrm{s} }\sum n_k \frac{\partial Z_k}{\partial \epsilon^\prime },\end{aligned} $$(A.6)

J 2 , 3 = f 2 n i = ( ϵ + ϵ 0 ) n s 1 q is ( n i + n s ) 2 1 ( n i + n s ) 2 n k Z k , $$ \begin{aligned} \mathbb{J} _{2,3} = \frac{\partial f_2}{\partial n_\mathrm{i} } =&- (\epsilon^\prime +\epsilon _0) n_\mathrm{s} \frac{1-q_\mathrm{is} }{(n_\mathrm{i} +n_\mathrm{s} )^2} - \frac{1}{(n_\mathrm{i} +n_\mathrm{s} )^2} \sum n_k Z_k,\end{aligned} $$(A.7)

J 2 , 4 = f 2 n s = ( ϵ + ϵ 0 ) n i q is 1 ( n i + n s ) 2 1 ( n i + n s ) 2 n k Z k , $$ \begin{aligned} \mathbb{J} _{2,4} = \frac{\partial f_2}{\partial n_\mathrm{s} } =&- (\epsilon^\prime +\epsilon _0) n_\mathrm{i} \frac{q_\mathrm{is} -1}{(n_\mathrm{i} +n_\mathrm{s} )^2} - \frac{1}{(n_\mathrm{i} +n_\mathrm{s} )^2} \sum n_k Z_k, \end{aligned} $$(A.8)

J 3 , 1 = f 3 ψ = n i v i ( ζ + k s , i n s ) n H n k π a k 2 , $$ \begin{aligned} \mathbb{J} _{3,1} = \frac{\partial f_3}{\partial \psi } =&\frac{n_\mathrm{i} v_\mathrm{i} }{(\zeta +k_\mathrm{s,i} n_\mathrm{s} )n_\mathrm{H} } \sum n_k \pi a_k^2,\end{aligned} $$(A.9)

J 3 , 2 = f 3 ϵ = σ v ie n i ( n i + q is n s ) ( ζ + k s , i n s ) n H n i v i ( ζ + k s , i n s ) n H n k π a k 2 J k ϵ , $$ \begin{aligned} \mathbb{J} _{3,2} = \frac{\partial f_3}{\partial \epsilon^\prime } =&-\frac{\langle \sigma v \rangle _\mathrm{ie} n_\mathrm{i} (n_\mathrm{i} +q_\mathrm{is} n_\mathrm{s} )}{(\zeta +k_\mathrm{s,i} n_\mathrm{s} )n_\mathrm{H} } \nonumber \\&- \frac{n_\mathrm{i} v_\mathrm{i} }{(\zeta +k_\mathrm{s,i} n_\mathrm{s} )n_\mathrm{H} } \sum n_k \pi a_k^2 \frac{\partial J_k}{\partial \epsilon^\prime },\end{aligned} $$(A.10)

J 3 , 3 = f 3 n i = 1 n H σ v ie ( ϵ + ϵ 0 ) ( 2 n i + q is n s ) ( ζ + k s , i n s ) n H v i n k π a k 2 J ( τ k ) ( ζ + k s , i n s ) n H k i , s ( n s , 0 n s ) ( ζ + k s , i n s ) n H , $$ \begin{aligned} \mathbb{J} _{3,3} = \frac{\partial f_3}{\partial n_\mathrm{i} } =&-\frac{1}{n_\mathrm{H} } - \frac{\langle \sigma v \rangle _\mathrm{ie} (\epsilon^\prime +\epsilon _0) (2n_\mathrm{i} +q_\mathrm{is} n_\mathrm{s} )}{(\zeta +k_\mathrm{s,i} n_\mathrm{s} )n_\mathrm{H} } \nonumber \\&-\frac{v_\mathrm{i} \sum n_k \pi a_k^2 \langle \tilde{J}(\tau _k) \rangle }{(\zeta +k_\mathrm{s,i} n_\mathrm{s} )n_\mathrm{H} } - \frac{k_\mathrm{i,s} (n_{\mathrm{s} ,0}-n_\mathrm{s} )}{(\zeta +k_\mathrm{s,i} n_\mathrm{s} )n_\mathrm{H} }, \end{aligned} $$(A.11)

J 3 , 4 = f 3 n s = σ v ie ( ϵ + ϵ 0 ) n i ( q is ζ k s , i n i ) ( ζ + k s , i n s ) 2 n H + k s , i n i v i n k π a k 2 J ( τ k ) ( ζ + k s , i n s ) 2 n H + k i , s n i ( ζ + k s , i n s , 0 ) ( ζ + k s , i n s ) 2 n H $$ \begin{aligned} \mathbb{J} _{3,4} = \frac{\partial f_3}{\partial n_\mathrm{s} } =&- \frac{\langle \sigma v \rangle _\mathrm{ie} (\epsilon^\prime +\epsilon _0)n_\mathrm{i} (q_\mathrm{is} \zeta -k_\mathrm{s,i} n_\mathrm{i} )}{(\zeta +k_\mathrm{s,i} n_\mathrm{s} )^2 n_\mathrm{H} }\nonumber \\&+\frac{k_\mathrm{s,i} n_\mathrm{i} v_\mathrm{i} \sum n_k \pi a_k^2 \langle \tilde{J}(\tau _k) \rangle }{(\zeta +k_\mathrm{s,i} n_\mathrm{s} )^2n_\mathrm{H} }\nonumber \\&+\frac{k_\mathrm{i,s} n_\mathrm{i} (\zeta +k_\mathrm{s,i} n_{\mathrm{s} ,0})}{(\zeta +k_\mathrm{s,i} n_\mathrm{s} )^2n_\mathrm{H} } \end{aligned} $$(A.12)

J 4 , 1 = f 4 ψ = n s v s ( ζ + β T 1 2 e T 0 T n H ) n s , 0 n k π a k 2 , $$ \begin{aligned} \mathbb{J} _{4,1} = \frac{\partial f_4}{\partial \psi } =&\frac{n_\mathrm{s} v_\mathrm{s} }{(\zeta +\beta T^\frac{1}{2}e^{-\frac{T_0}{T}}n_\mathrm{H} )n_{\mathrm{s} ,0}} \sum n_k \pi a_k^2,\end{aligned} $$(A.13)

J 4 , 2 = f 4 ϵ = σ v se n s ( n i + q is n s ) ( ζ + β T 1 2 e T 0 T n H ) n s , 0 n s v s ( ζ + β T 1 2 e T 0 T n H ) n s , 0 n k π a k 2 J k ϵ , $$ \begin{aligned} \mathbb{J} _{4,2} = \frac{\partial f_4}{\partial \epsilon^\prime } =&-\frac{\langle \sigma v \rangle _\mathrm{se} n_\mathrm{s} (n_\mathrm{i} +q_\mathrm{is} n_\mathrm{s} )}{(\zeta +\beta T^\frac{1}{2}e^{-\frac{T_0}{T}}n_\mathrm{H} )n_{\mathrm{s} ,0}} \nonumber \\&- \frac{n_\mathrm{s} v_\mathrm{s} }{(\zeta +\beta T^\frac{1}{2}e^{-\frac{T_0}{T}}n_\mathrm{H} )n_{\mathrm{s} ,0}} \sum n_k \pi a_k^2 \frac{\partial J_{k}}{\partial \epsilon^\prime },\end{aligned} $$(A.14)

J 4 , 3 = f 4 n i = ( 1 n s n s , 0 ) k i , s ζ + β T 1 2 e T 0 T n H + σ v se ( ϵ + ϵ 0 ) n s + k s , i n s ( ζ + β T 1 2 e T 0 T n H ) n s , 0 , $$ \begin{aligned} \mathbb{J} _{4,3} = \frac{\partial f_4}{\partial n_\mathrm{i} } =&\left(1-\frac{n_\mathrm{s} }{n_{\mathrm{s} ,0}}\right)\frac{k_\mathrm{i,s} }{\zeta +\beta T^\frac{1}{2}e^{-\frac{T_0}{T}}n_\mathrm{H} } \nonumber \\&+ \frac{-\langle \sigma v \rangle _\mathrm{se} (\epsilon^\prime +\epsilon _0) n_\mathrm{s} +k_\mathrm{s,i} n_\mathrm{s} }{(\zeta +\beta T^\frac{1}{2}e^{-\frac{T_0}{T}}n_\mathrm{H} ) n_{\mathrm{s} ,0}},\end{aligned} $$(A.15)

J 4 , 4 = f 4 n s = 1 n s , 0 k i , s n i + σ v se ( ϵ + ϵ 0 ) ( n i + 2 q is n s ) ( ζ + β T 1 2 e T 0 T n H ) n s , 0 v s n k π a k 2 J ( τ k ) ( ζ + β T 1 2 e T 0 T n H ) n s , 0 k s , i ( n H n s , 0 n i ) ( ζ + β T 1 2 e T 0 T n H ) n s , 0 , $$ \begin{aligned} \mathbb{J} _{4,4} = \frac{\partial f_4}{\partial n_\mathrm{s} } =&-\frac{1}{n_{\mathrm{s} ,0}}-\frac{k_\mathrm{i,s} n_\mathrm{i} +\langle \sigma v \rangle _\mathrm{se} (\epsilon^\prime +\epsilon _0)(n_\mathrm{i} +2q_\mathrm{is} n_\mathrm{s} )}{(\zeta +\beta T^\frac{1}{2}e^{-\frac{T_0}{T}}n_\mathrm{H} )n_{\mathrm{s} ,0}} \nonumber \\&- \frac{v_\mathrm{s} \sum n_k \pi a_k^2 \langle \tilde{J}(\tau _k) \rangle }{(\zeta +\beta T^\frac{1}{2}e^{-\frac{T_0}{T}}n_\mathrm{H} )n_{\mathrm{s} ,0}} - \frac{k_\mathrm{s,i} (n_\mathrm{H} -n_{\mathrm{s} ,0}-n_\mathrm{i} )}{(\zeta +\beta T^\frac{1}{2}e^{-\frac{T_0}{T}}n_\mathrm{H} ) n_{\mathrm{s} ,0}}, \end{aligned} $$(A.16)

with

Z k ϵ = ( ϵ + ϵ 0 ) 2 Θ 3 α k 4 ( ϵ + ϵ 0 ) Θ 2 α k Θ ( 1 + α k ( ϵ + ϵ 0 ) Θ + ( ϵ + ϵ 0 ) 2 Θ 2 ) 2 , $$ \begin{aligned} \frac{\partial Z_k}{\partial \epsilon^\prime } =&\frac{-(\epsilon^\prime +\epsilon _0)^2 \Theta ^3 \alpha _k - 4(\epsilon^\prime +\epsilon _0) \Theta ^2 - \alpha _k \Theta }{(1+\alpha _k (\epsilon^\prime +\epsilon _0) \Theta +(\epsilon^\prime +\epsilon _0)^2 \Theta ^2)^2},\end{aligned} $$(A.17)

J k ϵ = 2 τ k ( α k 1 ) ( ϵ + ϵ 0 ) 2 Θ 3 + Θ + 2 ( ϵ + ϵ 0 ) Θ 2 ( 1 + α k ( ϵ + ϵ 0 ) Θ + ( ϵ + ϵ 0 ) 2 Θ 2 ) 2 . $$ \begin{aligned} \frac{\partial J_k}{\partial \epsilon^\prime } =&\frac{2}{\tau _k} \frac{(\alpha _k-1)(\epsilon^\prime +\epsilon _0)^2 \Theta ^3 + \Theta + 2(\epsilon^\prime +\epsilon _0) \Theta ^2}{(1+\alpha _k (\epsilon^\prime +\epsilon _0) \Theta +(\epsilon^\prime +\epsilon _0)^2 \Theta ^2)^2}. \end{aligned} $$(A.18)

All Tables

Table 1.

Constants depending on the choice of species s for sodium, potassium, and hydrogen.

All Figures

thumbnail Fig. 1.

Test of our method for protostellar collapse conditions. Evolution of the average charge of grains (left axis) for several grain sizes (color lines), and the fractional abundance (right axis) of cosmic-ray ionized ions ni/nH (solid line), thermally ionized K+ ions ns/nH (dotted line), and electrons ne/nH (dashed line).

In the text
thumbnail Fig. 2.

Same as Fig. 1, but for Na. The only ions present are the cosmic-ray ionized ones i and Na (ion s), which is both cosmic-ray and thermally ionized.

In the text
thumbnail Fig. 3.

Same as Fig. 1, but for H. The only ions presents are the cosmic-ray ionized ones i and H (ions s), which is both cosmic-ray and thermally ionized.

In the text
thumbnail Fig. 4.

Evolution of the Ohmic (yellow), Hall (blue), and ambipolar (dark red) resistivities for the thermal ionization of K. The dashed line represents the Hall resistivity in negative values. The thin dotted lines are the resistivities in the absence of thermal ionization.

In the text
thumbnail Fig. 5.

Same as Fig. 4, but for Na.

In the text
thumbnail Fig. 6.

Same as Fig. 4, but for H.

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.