# Fast methods to track grain coagulation and ionization.
## II. Extension to thermal ionization

Pierre Marchand, Vincent Guillet, Ugo Lebreuilly, Mordecai-Mark Mac Low

► To cite this version:

Pierre Marchand, Vincent Guillet, Ugo Lebreuilly, Mordecai-Mark Mac Low. Fast methods to track grain coagulation and ionization. II. Extension to thermal ionization. Astronomy and Astrophysics - A&A, 2022, 666, pp.A27. 10.1051/0004-6361/202142551. hal-03587594

HAL Id: hal-03587594

https://hal.science/hal-03587594

Submitted on 9 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fast methods for tracking grain coagulation and ionization

## II. Extension to thermal ionization

P. Marchand<fnref n="1" />,<fnref n="2" />, V. Guillet<fnref n="3" />,<fnref n="4" />, U. Lebreuilly<fnref n="5" />, and M.-M. Mac Low<fnref n="2" />,⁶

Received 29 October 2021 / Accepted 18 February 2022

### ABSTRACT

Thermal ionization is a critical process at temperatures $T > 10^3$ K, particularly during star formation. An increase in ionization leads to a decrease in nonideal magnetohydrodynamics (MHD) resistivities, which has a significant impact on protoplanetary disks and protostar formation. We developed an extension of the fast computational ionization method presented in our recent paper to include thermal ionization. The model can be used to inexpensively calculate the density of ions and electrons and the electric charge of each size of grains for an arbitrary size distribution. This tool should be particularly useful for the self-consistent calculation of nonideal MHD resistivities in multidimensional simulations, especially of protostellar collapse and protoplanetary disks.

**Key words.** dust, extinction – stars: formation – ISM: abundances

## 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 $H_2$ molecules absorbs energy, allowing a rapid second collapse that leads to the creation of the protostar when all $H_2$ 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 temper-

atures >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 chondrule 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 $n_i$ and $n_s$. 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 $\mu_i = m_i/m_H$. Species *s*, however, undergoes both cosmic-ray and thermal ionization. We define $n_{s,0}$ as the total abundance of species *s* (both neutral and charged), so that $n_{s,0}$ is an upper bound for $n_s$. We also consider an arbitrary size distribution of dust grains.

<footnote n="1">¹ Institut de Recherche en Astrophysique et Planétologie, Université Paul, Sabatier Toulouse 3, 118 Rte de Narbonne, 31062 Toulouse, France</footnote>

<footnote n="2">² Department of Astrophysics, American Museum of Natural History, Central Park West at 79th Street, New York, NY 10024, USA
e-mail: pierre.marchand.astr@gmail.com</footnote>

<footnote n="3">³ Université Paris-Saclay, CNRS, Institut d'astrophysique spatiale, 91405 Orsay, France</footnote>

<footnote n="4">⁴ Laboratoire Univers et Particules de Montpellier, Université de Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France</footnote>

<footnote n="5">⁵ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette, France</footnote>
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, \tau_k)(J_i(Z, \tau_k) + J_s(Z, \tau_k)) = f(Z + 1, \tau_k)J_e(Z + 1, \tau_k), \quad (1)$$

where

$$\tau_k = \frac{a_k k_B T}{e^2} \qquad (2)$$

is the reduced temperature of a grain of radius $a_k$ and temperature $T$, while $k_B$ is the Boltzmann constant. Furthermore, $f(Z, \tau_k)$ is the distribution function of grain charges for a grain of reduced temperature $\tau_k$, and $J_i$, $J_s$, and $J_e$ are the fluxes of species i, s, and electrons onto the grains, respectively. The fluxes are (Draine & Sutin 1987)

$$J_j(Z, \tau_k) = n_j s_j v_j \pi a_k^2 \tilde{J}(Ze/q_j, \tau_k), \quad (3)$$

with $n_j$ being the abundance of species j (for $j = i, s$, or $e$), $s_j$ being the sticking probability of species j on grains, $m_j$ and $q_j$ being the mass and charge of species j, and $v_j = (8k_B T/\pi m_j)^{1/2}$ being the thermal speed of j. The polarization factor of grains for species j is $\tilde{J}$, which depends on the relative signs of Z and $q_j$.

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 $\tau_k \ll 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, \tau_k) + f(0, \tau_k) + f(1, \tau_k) = 1. \qquad (4)$$

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

$$f(-1, \tau_k) = f(0, \tau_k) \frac{J_e(0, \tau_k)}{J_i(-1, \tau_k) + J_s(-1, \tau_k)}, \quad (5)$$

$$f(1, \tau_k) = f(0, \tau_k) \frac{J_i(0, \tau_k) + J_s(0, \tau_k)}{J_e(1, \tau_k)}. \quad (6)$$

As in Paper I, we have

$$\tilde{J}(0, \tau_k) \approx (\pi/2\tau_k)^{\frac{1}{2}}, \quad (7)$$

$$\tilde{J}(-1, \tau_k) \approx (2/\tau_k), \quad (8)$$

where $\tilde{J}(0, \tau_k)$ appears in $J_e(0, \tau_k)$, $J_i(0, \tau_k)$, and $J_s(0, \tau_k)$, while $\tilde{J}(-1, \tau_k)$ appears in $J_e(1, \tau_k)$, $J_i(-1, \tau_k)$, and $J_s(-1, \tau_k)$. We can then solve the equation system (4)–(6). With $s_i = s_s = 1$, we obtain

$$f(0, \tau_k) = \frac{1}{1 + \frac{1}{\alpha_k} \left[ \epsilon \Theta + \frac{1}{\epsilon \Theta} \right]}, \quad (9)$$

$$f(-1, \tau_k) = \frac{\epsilon^2 \Theta^2}{1 + \alpha_k \epsilon \Theta + \epsilon^2 \Theta^2}, \quad (10)$$

$$f(1, \tau_k) = \frac{1}{1 + \alpha_k \epsilon \Theta + \epsilon^2 \Theta^2}, \quad (11)$$

where $\Theta = s_e(\mu_i m_H / m_e)^{1/2}$, $q_{is} = (\mu_i / \mu_s)^{1/2}$, and

$$\epsilon = \frac{n_e}{n_i + q_{is} n_s}. \quad (12)$$

The main difference with Paper I is the appearance of the term $q_{is}$. While $n_i + n_s$ is the total abundance of ions, $n_i + q_{is} n_s$ is

an effective abundance that reflects the relative flux of ions onto grains. The average charge of grains

$$Z_k = \sum_{Z=-1}^{Z=1} Z f(Z, \tau_k) \qquad (13)$$

and the grain-ion recombination enhancement factor

$$\langle \tilde{J}(\tau_k) \rangle = \sum_{Z=-1}^{Z=1} \tilde{J}(Ze/q_i, \tau_k) f(Z, \tau_k) \qquad (14)$$

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

$$Z_k = \frac{1 - \epsilon^2 \Theta^2}{1 + \alpha_k \epsilon \Theta + \epsilon^2 \Theta^2}, \quad (15)$$

$$\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}, \quad (16)$$

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

For the larger grains ($\tau_k \gg 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 $\psi$. For $\psi < 0$, $e^\psi$ represents the repulsion of the flux of electrons by the negatively charged grains, while $1-\psi$ characterizes the attraction of the flux of ions. The flux equilibrium can thus be written as

$$(n_e v_e s_e) e^\psi = (n_i v_i s_i + n_s v_s s_s)(1 - \psi). \quad (17)$$

Introducing the same notations as above, we can write

$$\epsilon = \frac{1 - \psi}{\Theta e^{\psi}}, \quad (18)$$

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

$$Z_k = \psi \tau_k, \quad (19)$$

$$\langle \tilde{J}(\tau_k) \rangle = (1 - \psi). \quad (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 + \sum n_k Z_k = 0. \quad (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 $\mu_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

$$\zeta(n_H - n_{s,0} - n_i) + k_{s,i} (n_H - n_{s,0} - n_i) n_s = \langle \sigma v \rangle_{ie} n_e n_i + n_i v_i \sum n_k \pi a_k^2 J_k + k_{i,s} (n_{s,0} - n_s) n_i, \quad (22)$$
where $\zeta$ is the cosmic-ray ionization rate, $k_{s,i}$ and $k_{i,s}$ are the chemical reaction rates of species $s \to i$ and $i \to s$, respectively, and we consider the recombination rate of ions $i$ with electrons to be $\langle \sigma v \rangle_{ie} = 2 \times 10^{-7} (T/300)^{1/3} \text{ cm}^3 \text{ s}^{-1}$, based on the recombination rate of HCO$^+$ taken from the UMIST database (McElroy et al. 2013). The terms of the form $k_{a,b} n_a n_c$ are the transformation of species $c$ to species $b$, through the chemical reaction with species $a$. Hence the term $(n_H - n_{s,0} - n_i)$ represents the total abundance of neutral species that can be ionized into ion $i$, as $(n_{s,0} - n_s)$ 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)

$$ \left( \frac{dn_s}{dt} \right)_{\text{thermal ionization}} = \beta_s n_H T^{1/2} e^{-T_{0,s}/T} (n_{s,0} - n_s), \quad (23) $$

with the values of $\beta$ and $T_0$ 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

$$ \begin{aligned} \zeta(n_{s,0} - n_s) + \left(\frac{dn_s}{dt}\right)_{\text{thermal ionization}} &+ k_{i,s}(n_{s,0} - n_s)n_i \\ = \langle \sigma v \rangle_{\text{se}} n_e n_s + n_s v_s \sum n_k \pi a_k^2 J_k + k_{s,i}(n_H - n_{s,0} - n_i)n_s. \end{aligned} \quad (24) $$

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

$$ \frac{4}{3} \pi n_H \lambda_D^3 \gg 1, \quad (25) $$

with

$$ \lambda_D = \sqrt{\frac{k_B T}{4 \pi n_i e^2}}. \quad (26) $$

### 3. Numerical implementation and tests

Equations (18), (21), (22), and (24) need to be solved for $\psi$, $\epsilon$, $n_i$, and $n_s$. 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 $\psi$ as a function of $\epsilon$. 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, $\psi$ converges toward zero by a negative value and becomes very small in an absolute value. At high temperature, $n_s$ overwhelmingly dominates $n_i$ due to the thermal ionization. Therefore, the numerical implementation has to be robust for the cases $|\psi| \ll 1$ and $n_s/n_i \gg 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 $n_s$ grows much larger than $n_i$. We therefore normalized Eq. (21) by the total number of ions $n_i + n_s$, we used both the cosmic-ray ionization rate and the chemical reaction rate $s \to 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 $n_i \gg n_s$, $\epsilon$ converges toward $1/\Theta$ as density increases (see Fig. 2 of Paper I). When the difference between $\epsilon$ and $1/\Theta$ becomes close to machine precision, the term $1 - (\epsilon\Theta)^2$ reaches a lower bound<fnref n="1" /> 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 $\epsilon$ by $\epsilon' + \epsilon_0$, with $\epsilon_0 = 1/\Theta$. This substitution has to be made in all the equations. The new variable to find is $\epsilon'$ which converges toward zero instead of $1/\Theta$. The term $1 - (\epsilon\Theta)^2$ in Eq. (15) is then mathematically equal to and can be replaced by $-2\epsilon'\Theta - (\epsilon'\Theta)^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(\psi, \epsilon', n_i, n_s) = (f_1, f_2, f_3, f_4)$, with

$$ f_1 = \frac{1 - \psi}{(\epsilon' + \epsilon_0)\Theta e^{\psi}} - 1 = 0, \quad (27) $$

$$ f_2 = 1 - (\epsilon' + \epsilon_0) \frac{n_i + q_{is} n_s}{n_i + n_s} + \frac{1}{n_i + n_s} \sum n_k Z_k = 0, \quad (28) $$

$$ \begin{aligned} f_3 = 1 - & \frac{n_{s,0} + n_i}{n_H} - \frac{\langle \sigma v \rangle_{\text{ic}} (\epsilon' + \epsilon_0) n_i (n_i + q_{is} n_s)}{(\zeta + k_{s,i} n_s) n_H} \\ & - \frac{n_i v_i \sum n_k \pi a_k^2 \langle \tilde{J}(\tau_k) \rangle}{(\zeta + k_{s,i} n_s) n_H} - \frac{k_{i,s} n_i (n_{s,0} - n_s)}{(\zeta + k_{s,i} n_s) n_H} = 0, \end{aligned} \quad (29) $$

$$ \begin{aligned} f_4 = & \left(1 - \frac{n_s}{n_{s,0}}\right) \left(1 + \frac{k_{i,s} n_i}{\zeta + \beta T^{\frac{1}{2}} e^{-\frac{T_0}{T}} n_H}\right) \\ & - \frac{\langle \sigma v \rangle_{\text{se}} (\epsilon' + \epsilon_0) n_s (n_i + q_{is} n_s)}{(\zeta + \beta T^{\frac{1}{2}} e^{-\frac{T_0}{T}} n_H) n_{s,0}} \\ & - \frac{n_s v_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_H) n_{s,0}} \\ & - \frac{k_{s,i} n_s (n_H - n_{s,0} - n_i)}{(\zeta + \beta T^{\frac{1}{2}} e^{-\frac{T_0}{T}} n_H) n_{s,0}} = 0, \end{aligned} \quad (30) $$

where

$$ Z_k = \psi \tau_k + \frac{-2\epsilon'\Theta - (\epsilon'\Theta)^2}{1 + \alpha_k (\epsilon' + \epsilon_0)\Theta + (\epsilon' + \epsilon_0)^2\Theta^2}, \quad (31) $$

$$ \langle \tilde{J}(\tau_k) \rangle = (1 - \psi) + \frac{\frac{2}{\tau_k}((\epsilon' + \epsilon_0)^2\Theta^2 + (\epsilon' + \epsilon_0)\Theta) + \alpha_k(\epsilon' + \epsilon_0)\Theta + 1}{(\epsilon' + \epsilon_0)^2\Theta^2 + \alpha_k(\epsilon' + \epsilon_0)\Theta + 1}, \quad (32) $$

$$ \epsilon = \frac{n_e}{n_s} = \epsilon' + \epsilon_0, \quad (33) $$

$$ \epsilon_0 = \frac{1}{\Theta}, \quad (34) $$

<footnote n="1">¹ For example, for a machine precision of $10^{-16}$, if $|\epsilon - 1/\Theta| < 10^{-10}$, then $1 - (\epsilon\Theta)^2 > 10^{-6}$.</footnote>
<table><caption>**Table 1.** Constants depending on the choice of species *s* for sodium, potassium, and hydrogen.</caption><thead><tr><th>Species</th><th>μ<sub>s</sub></th><th>n<sub>s,0</sub>/n<sub>H</sub></th><th>β<sub>s</sub> (cm<sup>3</sup> s<sup>-1</sup> K<sup>-1/2</sup>)</th><th>T<sub>0,s</sub> (K)</th><th>k<sub>i,s</sub> (cm<sup>3</sup> s<sup>-1</sup>)</th><th>k<sub>s,i</sub> (cm<sup>3</sup> s<sup>-1</sup>)</th><th>⟨σv⟩<sub>se</sub> (cm<sup>3</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Na</td><td>22.99</td><td>3.1 × 10<sup>-9</sup></td><td>1.4 × 10<sup>-15</sup></td><td>6.0 × 10<sup>4</sup></td><td>6.2 × 10<sup>-9</sup></td><td>0</td><td>2.78 × 10<sup>-12</sup>(T/300)<sup>-0.68</sup></td></tr><tr><td>K</td><td>39.09</td><td>2.2 × 10<sup>-10</sup></td><td>6.5 × 10<sup>-15</sup></td><td>5.1 × 10<sup>4</sup></td><td>6.2 × 10<sup>-9</sup></td><td>0</td><td>2.78 × 10<sup>-12</sup>(T/300)<sup>-0.68</sup></td></tr><tr><td>H</td><td>1.0</td><td>1.0</td><td>2.0 × 10<sup>-10</sup></td><td>15.8 × 10<sup>4</sup></td><td>3.7 × 10<sup>-14</sup></td><td>3.8 × 10<sup>-9</sup></td><td>3.5 × 10<sup>-12</sup>(T/300)<sup>-0.75</sup></td></tr></tbody></table>

**Notes.** The values of $\beta_s$ and $T_{0,s}$ are taken from Pneuman & Mitchell (1965) (see Sect. 4 for an important discussion about those rates), while $n_{s,0}/n_H$ are the same as in Umebayashi & Nakano (1990). The values of $k_{i,s}$, $k_{s,i}$, and $\langle\sigma v\rangle_{se}$ are taken from the UMIST database (McElroy et al. 2013), by summing the rates over all reactions involving those species in the reduced network of Marchand et al. (2016).

and

$$v_{i,s} = \left( \frac{8k_B T}{\pi \mu_{i,s} m_H} \right)^{1/2} . \qquad (35)$$

Here, $v_{i,s}$ and $\mu_{i,s}$ stand for either $v_i$ and $\mu_i$ or $v_s$ and $\mu_s$.

### 3.3. Solving the system

Similarly to Paper I, we used a Newton-Raphson method to solve the equation system. Let $\mathbf{X} = (\psi, \epsilon, n_i, n_s)$. Starting from an educated guess $\mathbf{X}_0$, we iterated

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

until $\|\mathbf{F}(\mathbf{X}_n)\| < \delta \ll 1$. The matrix $J$ is the Jacobian of the system defined by $J_{i,j} = \partial f_i / \partial X_j$. 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 ($n_H < 10^7$ cm$^{-3}$) and low temperature (typically 10 K), a good starting point is

$$\psi = \psi_0, \qquad (37)$$

$$\epsilon = 0.9999, \qquad (38)$$

$$n_i = 10^{-7} n_H, \qquad (39)$$

$$n_s = 10^{-4} n_i, \qquad (40)$$

where $\psi_0$ is the solution of

$$\frac{1 - \psi}{\Theta e^{\psi}} = 1. \qquad (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, $n_s$ 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 $n_H = 10^4-10^{25}$ cm$^{-3}$, starting from the lowest density and increasing $n_H$ 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 \left( 1 + \left[ \frac{n_H}{n_1} \right]^{0.8} \right)^{\frac{1}{2}} \left( 1 + \left[ \frac{n_H}{n_2} \right] \right)^{-0.3} \left( 1 + \left[ \frac{n_H}{n_3} \right] \right)^{-\frac{1.7}{3}}, \qquad (42)$$

with $T_0 = 10$ K, $n_1 = 10^{11}$ cm$^{-3}$, $n_2 = 10^{16}$ cm$^{-3}$, and $n_3 = 10^{21}$ 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 $a_{min} = 5$ nm and maximum grain size $a_{max} = 250$ nm, sampled by 26 bins. The grain bulk density is $\rho = 2.9$ g cm$^{-3}$ and the dust-to-gas mass ratio is 1%, which is typical of the interstellar medium (Bohlin et al. 1978). We set $\zeta = 5 \times 10^{-17}$ s$^{-1}$ (Padovani et al. 2013), $s_e = 0.5$ (Umebayashi & Nakano 1990), and $\mu_i = 25$ (close to the molecular mass of Mg, Fe, or HCO$^+$, Marchand et al. 2016).

In this test, the temperature exceeds several $10^3$ 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 $n_H \approx 10^{16}-10^{19}$ 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 $n_{s,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 $n_{s,0}$ in the table. Here, we present the evolution of $n_i/n_H$, $n_s/n_H$, $n_e/n_H$ and the average grain charge for several bins. The results are displayed in Figs. 1–3, respectively.

These figures can be compared to Fig. 7 of Marchand et al. (2016). Although the abundances of Na$^+$ and K$^+$ seem overestimated compared to $n_e$, 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 $n_H = 10^{18}$ cm$^{-3}$ and $T = 1600$ K, and they saturate around $n_H = 10^{22}$ cm$^{-3}$ and $T = 1.5 \times 10^4$ K at the maximum fractional abundance of their respective species. Furthermore, H has a similar behavior, starting its thermal ionization at $n_H = 10^{20}$ cm$^{-3}$ and $T = 2650$ K, turning virtually all neutrals into ions by $n_H = 10^{23}$ cm$^{-3}$ and $T = 7 \times 10^4$ K. The main difference with Fig. 7 of Marchand et al. (2016) is the earlier rise of electron density at $n_H \approx 10^{16}$ 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 $n_{H^+}/n_H > 0.1$ reached at
<figure><img src="image_2.png" /><figcaption>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 n<sub>i</sub>/n<sub>H</sub> (solid line), thermally ionized K<sup>+</sup> ions n<sub>s</sub>/n<sub>H</sub> (dotted line), and electrons n<sub>e</sub>/n<sub>H</sub> (dashed line).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Same as Fig. 1, but for H. The only ions present are the cosmic-ray ionized ones i and H (ions s), which is both cosmic-ray and thermally ionized.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Same as Fig. 4, but for Na.</figcaption></figure>

$n_H > 10^{23} \text{ 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 4–6 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 \times 10^{-7} \text{ G} n_H^{1/2}, \qquad (43)$$

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 $n_H \ge 10^{12} \text{ cm}^{-3}$.

The nonideal MHD terms significantly affect the MHD evolution of protostellar collapse and protoplanetary disks for resistivities larger than $\approx 10^{18} \text{ cm}^2 \text{ 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).
<figure><img src="image_7.png" /><figcaption>Fig. 6. Same as Fig. 4, but for H.</figcaption></figure>

## 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 $\zeta$ 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 & Sugden 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 H$_2$. 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:

$$ \beta_K = 1.0 \times 10^{-8} \text{ cm}^3 \text{s}^{-1} \text{ K}^{-\frac{1}{2}}, \qquad (44) $$

$$ T_{0,K} = 5.6 \times 10^4 \text{ K}, \qquad (45) $$

$$ \beta_{\text{Na}} = 2.0 \times 10^{-9} \text{ cm}^3 \text{s}^{-1} \text{ K}^{-\frac{1}{2}}, \qquad (46) $$

$$ T_{0,\text{Na}} = 5.6 \times 10^4 \text{ K}. \qquad (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.

*Acknowledgements.* 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

Aller, L. H. 1961, The Abundances of the Elements (New York: Interscience Publ)

Ashton, A., & Hayhurst, A. 1973, Combust. Flame, 21, 69

Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132

Desch, S. J., & Turner, N. J. 2015, ApJ, 811, 156

Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803

Ebel, D. S., Hubbard, A., McNally, C., et al. 2012, Meteorit. Planetary Sci. Suppl., 75, 5387

Getman, K. V., & Feigelson, E. D. 2021, ApJ, 916, 32

Hollander, T., Kalff, P. J., & Alkemade, C. T. J. 1963, J. Chem. Phys., 39, 2558

Hollenbach, D. J., & Salpeter, E. E. 1969, J. Chem. Phys., 50, 4157

Hubbard, A., McNally, C. P., & Mac Low, M.-M. 2012, ApJ, 761, 58

Joung, M. K. R., Mac Low, M.-M., & Ebel, D. S. 2004, ApJ, 606, 532

Koga, S., Tsukamoto, Y., Okuzumi, S., & Machida, M. N. 2019, MNRAS, 484, 2119

Kunz, M. W., & Mouschovias, T. C. 2010, MNRAS, 408, 322

Larson, R. B. 1969, MNRAS, 145, 271

Lenzuni, P., Gail, H.-P., & Henning, T. 1995, ApJ, 447, 848

Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180

Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18

Marchand, P., Commerçon, B., & Chabrier, G. 2018, A&A, 619, A37

Marchand, P., Guillet, V., Lebreuilly, U., & Mac Low, M. M. 2021, A&A, 649, A50

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36

McNally, C. P., Hubbard, A., Mac Low, M.-M., Ebel, D. S., & D’Alessio, P. 2013, ApJ, 767, L2

McNally, C. P., Hubbard, A., Yang, C.-C., & Mac Low, M.-M. 2014, ApJ, 791, 62

Padovani, M., Hennebelle, P., & Galli, D. 2013, A&A, 560, A114

Pneuman, G. W., & Mitchell, T. P. 1965, Icarus, 4, 494

Schofield, K., & Sugden, T. M. 1965, Some Observations on the Ionization of Alkali and Alkaline-earth Elements in Hydrogen Flames (Cambridge: University of Cambridge), 10

Shui, V. H. 1977, Phys. Fluids, 20, 32

Spitzer, L., Jr 1949, Leaflet Astron. Soc. Pac., 5, 336

Teyssier, R. 2002, A&A, 385, 337

Umebayashi, T., & Nakano, T. 1990, MNRAS, 243, 103

Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037
Appendix A: Jacobian components

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

$$
\mathcal{J}_{1,1} = \frac{\partial f_1}{\partial \psi} = \frac{\psi - 2}{(\epsilon' + \epsilon_0)\Theta e^{\psi}}, \quad (A.1)
$$

$$
\mathcal{J}_{1,2} = \frac{\partial f_1}{\partial \epsilon'} = - \frac{1 - \psi}{(\epsilon' + \epsilon_0)^2 \Theta e^{\psi}}, \quad (A.2)
$$

$$
\mathcal{J}_{1,3} = \frac{\partial f_1}{\partial n_i} = 0, \qquad (A.3)
$$

$$
\mathcal{J}_{1,4} = \frac{\partial f_1}{\partial n_s} = 0, \qquad (A.4)
$$

$$
\mathcal{J}_{2,1} = \frac{\partial f_2}{\partial \psi} = \frac{1}{n_i + n_s} \sum n_k \tau_k, \qquad (A.5)
$$

$$
\mathcal{J}_{2,2} = \frac{\partial f_2}{\partial \epsilon'} = - \frac{n_i + q_{is} n_s}{n_i + n_s} + \frac{1}{n_i + n_s} \sum n_k \frac{\partial Z_k}{\partial \epsilon'}, \quad (A.6)
$$

$$
\mathcal{J}_{2,3} = \frac{\partial f_2}{\partial n_i} = - (\epsilon' + \epsilon_0) n_s \frac{1 - q_{is}}{(n_i + n_s)^2} - \frac{1}{(n_i + n_s)^2} \sum n_k Z_k, \quad (A.7)
$$

$$
\mathcal{J}_{2,4} = \frac{\partial f_2}{\partial n_s} = - (\epsilon' + \epsilon_0) n_i \frac{q_{is} - 1}{(n_i + n_s)^2} - \frac{1}{(n_i + n_s)^2} \sum n_k Z_k, \quad (A.8)
$$

$$
\mathcal{J}_{3,1} = \frac{\partial f_3}{\partial \psi} = \frac{n_i v_i}{(\zeta + k_{s,i} n_s) n_H} \sum n_k \pi a_k^2, \quad (A.9)
$$

$$
\mathcal{J}_{3,2} = \frac{\partial f_3}{\partial \epsilon'} = - \frac{\langle \sigma v \rangle_{\text{ie}} n_i (n_i + q_{\text{is}} n_s)}{(\zeta + k_{s,i} n_s) n_H} \\
\qquad - \frac{n_i v_i}{(\zeta + k_{s,i} n_s) n_H} \sum n_k \pi a_k^2 \frac{\partial J_k}{\partial \epsilon'}, \quad (A.10)
$$

$$
\mathcal{J}_{3,3} = \frac{\partial f_3}{\partial n_i} = - \frac{1}{n_H} - \frac{\langle \sigma v \rangle_{\text{ie}} (\epsilon' + \epsilon_0) (2n_i + q_{\text{is}} n_s)}{(\zeta + k_{s,i} n_s) n_H} \\
- \frac{v_i \sum n_k \pi a_k^2 \langle \mathcal{J}(\tau_k) \rangle}{(\zeta + k_{s,i} n_s) n_H} - \frac{k_{i,s} (n_{s,0} - n_s)}{(\zeta + k_{s,i} n_s) n_H}, \quad (A.11)
$$

$$
\begin{equation}
\begin{split}
\mathcal{J}_{3,4} = \frac{\partial f_3}{\partial n_s} = & - \frac{\langle \sigma v \rangle_{\text{ie}} (\epsilon' + \epsilon_0) n_i (q_{\text{is}} \zeta - k_{s,i} n_i)}{(\zeta + k_{s,i} n_s)^2 n_H} \\
& + \frac{k_{s,i} n_i v_i \sum n_k \pi a_k^2 \langle \mathcal{J}(\tau_k) \rangle}{(\zeta + k_{s,i} n_s)^2 n_H} \\
& + \frac{k_{i,s} n_i (\zeta + k_{s,i} n_{s,0})}{(\zeta + k_{s,i} n_s)^2 n_H}
\end{split}
\tag{A.12}
\end{equation}
$$

$$
\mathcal{J}_{4,1} = \frac{\partial f_4}{\partial \psi} = \frac{n_s v_s}{(\zeta + \beta T^{\frac{1}{2}} e^{-\frac{T_0}{T}} n_H) n_{s,0}} \sum n_k \pi a_k^2, \quad (A.13)
$$

$$
\begin{equation}
\begin{split}
\mathcal{J}_{4,2} = \frac{\partial f_4}{\partial \epsilon'} = & - \frac{\langle \sigma v \rangle_{\text{se}} n_s (n_i + q_{\text{is}} n_s)}{(\zeta + \beta T^{\frac{1}{2}} e^{-\frac{T_0}{T}} n_H) n_{s,0}} \\
& - \frac{n_s v_s}{(\zeta + \beta T^{\frac{1}{2}} e^{-\frac{T_0}{T}} n_H) n_{s,0}} \sum n_k \pi a_k^2 \frac{\partial J_k}{\partial \epsilon'}
\end{split}
\tag{A.14}
\end{equation}
$$

$$
\mathcal{J}_{4,3} = \frac{\partial f_4}{\partial n_i} = \left(1 - \frac{n_s}{n_{s,0}}\right) \frac{k_{i,s}}{\zeta + \beta T^{\frac{1}{2}} e^{-\frac{T_0}{T}} n_H} \\
+ \frac{-\langle \sigma v \rangle_{\text{se}} (\epsilon' + \epsilon_0) n_s + k_{s,i} n_s}{(\zeta + \beta T^{\frac{1}{2}} e^{-\frac{T_0}{T}} n_H) n_{s,0}}, \quad (A.15)
$$

$$
\begin{equation}
\begin{split}
\mathcal{J}_{4,4} = \frac{\partial f_4}{\partial n_s} = & - \frac{1}{n_{s,0}} - \frac{k_{i,s} n_i + \langle \sigma v \rangle_{\text{se}} (\epsilon' + \epsilon_0) (n_i + 2q_{is} n_s)}{(\zeta + \beta T^{\frac{1}{2}} e^{-\frac{T_0}{T}} n_H) n_{s,0}} \\
& - \frac{v_s \sum n_k \pi a_k^2 \langle \mathcal{J}(\tau_k) \rangle}{(\zeta + \beta T^{\frac{1}{2}} e^{-\frac{T_0}{T}} n_H) n_{s,0}} - \frac{k_{s,i} (n_H - n_{s,0} - n_i)}{(\zeta + \beta T^{\frac{1}{2}} e^{-\frac{T_0}{T}} n_H) n_{s,0}}
\end{split}
\tag{A.16}
\end{equation}
$$

with

$$
\frac{\partial Z_k}{\partial \epsilon'} = \frac{-(\epsilon' + \epsilon_0)^2 \Theta^3 \alpha_k - 4(\epsilon' + \epsilon_0)\Theta^2 - \alpha_k \Theta}{(1 + \alpha_k(\epsilon' + \epsilon_0)\Theta + (\epsilon' + \epsilon_0)^2 \Theta^2)^2}, \quad (A.17)
$$

$$
\frac{\partial J_k}{\partial \epsilon'} = \frac{2}{\tau_k} \frac{(\alpha_k - 1)(\epsilon' + \epsilon_0)^2 \Theta^3 + \Theta + 2(\epsilon' + \epsilon_0)\Theta^2}{(1 + \alpha_k(\epsilon' + \epsilon_0)\Theta + (\epsilon' + \epsilon_0)^2 \Theta^2)^2}. \quad (A.18)
$$