2.1. Cosmic ray ionisation

When a cosmic ray (or an X-ray) strikes an atom or molecule, it has the ability to remove an electron which creates an ion. Two ion species are modelled: species i with mass $m_{\text{i}_{\text{R}}}$ , which represents the light elements (e.g. hydrogen and helium compounds), and species I with mass $m_{\text{I}_{\text{R}}}$ , which represents the metals. The light ion mass is calculated from the given abundances or mass fractions of hydrogen and helium assuming all the hydrogen is molecular hydrogen, and the metals are modelled as a single species, where the mass is a free parameter.

Dust grains have the ability to absorb electrons, which gives the grains a negative charge or lose an electron to give them a positive charge. Nicil includes two dust models: fixed radius and MRN size distribution. For a fixed grain size, $a_{\text{g}}$ , the number density is proportional to the total number density, n (Keith & Wardle Reference Keith and Wardle2014),

(1) $$\begin{equation} n_{\text{g}} = \frac{m_{\text{n}}}{m_{\text{g}}}f_{\text{dg}} n, \end{equation}$$

where $f_{\text{dg}}$ is the dust-to-gas mass ratio and $m_{\text{n}}$ is the mass of a neutral particle.

The number density for the MRN grain distribution follows a power law, viz.,

(2) $$\begin{equation} \frac{\text{d}n_{\text{g}}}{\text{d}a} = An_{\text{H}}a^{-3.5}, \end{equation}$$

where n _H is the number density of the hydrogen nucleus, n _g(a) is the number density of grains with a radius smaller than a, and A = 1.5 × 10⁻²⁵ cm^2.5 (Draine & Lee Reference Draine and Lee1984). For simplicity, we assume n _H ≈ n, where n is the total number density. Following Nakano et al. (Reference Nakano, Nishi and Umebayashi2002), the given range of a _n < a < a _x is divided into N bins of fixed width Δlog a. Each bin has a characteristic grain size, a_j , given by its log-average grain size, and the number density of each bin is calculated using (2).

The mass of a grain particle is given by

(3) $$\begin{equation} m_{\rm g} = \frac{4}{3}\pi a_{\rm g}^3 \rho _{\rm b}, \end{equation}$$

where ρ_b is the grain bulk density.

The ion and electron number densities calculated by this process vary as (e.g. Umebayashi & Nakano Reference Umebayashi and Nakano1980; Fujii, Okuzumi, & Inutsuka Reference Fujii, Okuzumi and Inutsuka2011)

(4) $$\begin{eqnarray} \frac{{\rm d}n_s}{{\rm d}t} &=& \zeta n_{\rm n} - k_{{\rm e}s}n_s n_{{\rm e}_{\rm R}} \nonumber \\ &&- \sum _{Z=-1}^1 \sum _j k_{s{\rm g}}(Z,a_j)n_s n_{\rm g}(Z,a_j), \end{eqnarray}$$

(5) $$\begin{eqnarray} \frac{{\rm d}n_{{\rm e}_{\rm R}}}{{\rm d}t}&=& \zeta n_{\rm n} - \sum _s k_{{\rm e}s}n_s n_{{\rm e}_{\rm R}} \nonumber \\ &&- \sum _{Z=-1}^1 \sum _j k_{\rm eg}(Z,a_j)n_{{\rm e}_{\rm R}}n_{\rm g}(Z,a_j), \end{eqnarray}$$

(6) $$\begin{eqnarray} \frac{{\rm d}n_{\rm g}(Z,a_j)}{{\rm d}t}&=& \ -\sum _s \left[k_{s{\rm g}}(Z,a_j)n_s n_{\rm g}(Z,a_j) \right. \nonumber \\ &&- \left.k_{s{\rm g}}(Z-1,a_j)n_s n_{\rm g}(Z-1,a_j) \right] \nonumber \\ &&- k_{{\rm eg}}(Z,a_j)n_{\rm e} n_{\rm g}(Z,a_j) \nonumber \\ &&+ k_{{\rm eg}}(Z+1,a_j)n_{\rm e} n_{\rm g}(Z+1,a_j), \end{eqnarray}$$

where s ∈ {i_R, I_R}, ζ is the ionisation rate, k_ij are the charge capture rates, and we set n _g(Z = ±2, a_j ) ≡ 0 since their populations will be small compared the singly charged grains. Although we only include two ion species for simplicity and performance, it is possible to expand this network to include an arbitrary number of molecular ions.

The ion–electron charge capture rates are given by (Umebayashi & Nakano Reference Umebayashi and Nakano1990)

(7) $$\begin{eqnarray} k_{\rm ei} &=& \left[ 3.5 X \left(\frac{T}{300}\right)^{-0.7} + 4.5 Y \left(\frac{T}{300}\right)^{-0.67}\right] \nonumber \\ &&\times 10^{-12}\, {\rm cm}^3\,{\rm s}^{-1} \end{eqnarray}$$

(8) $$\begin{eqnarray} k_{\rm eI} = 2.8\times 10^{-12}\left(\frac{T}{300}\right)^{-0.86} {\rm cm}^3\,{\rm s}^{-1} , \end{eqnarray}$$

where X and Y are the mass fractions of hydrogen and helium, respectively, used to determine the mass of the light ion. For grains, the charge capture rates are (Fujii et al. Reference Fujii, Okuzumi and Inutsuka2011)

(9) $$\begin{eqnarray} k_{s{\rm g}}(Z,a_j) &=& a_j^2 \sqrt{\frac{8\pi k_{\rm B}T}{m_s}}\left\lbrace \begin{array}{l l}\left(1-\frac{Z e^2}{a_j k_{\rm B}T}\right) & \quad{\rm if }\quad Z \le 0\\ \exp \left(\frac{-Z e^2}{a_jk_{\rm B}T}\right) & \quad{\rm if }\quad Z > 0 \end{array}\right. , \end{eqnarray}$$

(10) $$\begin{eqnarray} k_{\rm eg}(Z,a_j) &=&a_j^2 \sqrt{\frac{8\pi k_{\rm B}T}{m_{\rm e}}} \left\lbrace \begin{array}{l l}\exp \left(\frac{Z e^2}{a_j k_{\rm B}T}\right) & \quad{\rm if }\quad Z < 0 \\ \left(1+\frac{Z e^2}{a_j k_{\rm B}T}\right) & \quad{\rm if }\quad Z \ge 0 \end{array}\right., \end{eqnarray}$$

where k _B is the Boltzmann constant, e is the electron charge, m _e is the mass of an electron, Z is the charge on the grains, and T is the temperature (where the gas and dust are assumed to be in thermal equilibrium).

From charge neutrality,

(11) $$\begin{equation} \sum _s n_s - n_{{\rm e}_{\rm R}} + \sum _{Z=-1}^1 \sum _j Z n_{\rm g}(Z,a_j) = 0, \end{equation}$$

and from conservation of grains,

(12) $$\begin{equation} n_{\rm g}(a_j) = \sum _{Z=-1}^1 n_{\rm g}(Z,a_j), \end{equation}$$

where n _g(a_j ) is the total grain number density as calculated in (1) or (2).

Following Keith & Wardle (Reference Keith and Wardle2014), we assume that the system is approximately in a steady state system (i.e. $\frac{{\rm d}}{{\rm d}t} \approx 0$ ). The set of equations, (4)–(6), (11), and (12) represents an over-subscription, thus we remove $\frac{{\rm d}n_{{\rm e}_{\rm R}}}{{\rm d}t}$ and $\frac{{\rm d}n_{\rm g}(Z=0,a_j)}{{\rm d}t}$ in favour of charge neutrality and conservation of grains. The number densities are calculated by iteratively solving

(13) $$\begin{equation} \bm {n}^{t+1} = \bm {n}^t - \mathbb {J}^{-1}(\bm {n}^t) \bm {f}(\bm {n}^t) = \bm {n}^t - \bm {x}^t, \end{equation}$$

where t is the iteration number, $\mathbb {J}(\bm {n}^t)$ is the Jacobian, and in the default case, n ^t = {n^t _iR , n^t _IR , n^t _e, n _g ^t (Z = −1), n^t _g(Z = 0), n _g ^t (Z = −1)} and f ( n ^t ) = 0 are $\lbrace \frac{{\rm d}n_{{\rm i}_{\rm R}}}{{\rm d}t}$ , $\frac{{\rm d}n_{{\rm I}_{\rm R}}}{{\rm d}t}$ , charge neutrality, $\frac{{\rm d}n_{\rm g}(Z=-1)}{{\rm d}t}$ , grain conservation, $\frac{{\rm d}n_{\rm g}(Z=1)}{{\rm d}t}\rbrace$ . Rather than solving for inverse of the Jacobian, we use LU decomposition of the Jacobian to solve for x ^t in $\mathbb {J}(\bm {n}^t) \bm {x}^t = \bm {f}(\bm {n}^t)$ .

The number densities can span several orders of magnitude, and this method becomes unstable for n _g ≫ n _e. Although this will likely not be important in most physical applications, an alternative method using average grain charges is described in Appendix D; this method is only used after the above method fails to iterate to a solution within a set number of iterations.

2.2. Thermal ionisation

At high gas temperatures, thermal ionisation can remove electrons. For each species, j, the ion number densities can be calculated using the Saha equation,

(14) $$\begin{equation} \frac{n_{{\rm e}_{\rm T}} n_{{\rm i}_{\rm T},j,k+1}}{n_{{\rm i}_{\rm T},j,k}} = \frac{2g_{j,k+1}}{g_{j,k}} \left( \frac{2\pi m_{\rm e} k_{\rm B} T}{h^2}\right)^{3/2} \exp {\left( -\frac{\chi _{j,k+1}}{k_{\rm B}T} \right)}, \end{equation}$$

where n _{iT, j, k} is the number density of species j which has been ionised k times, χ _{j, k} is the ionisation potential, g _{j, k} is the statistical weight of level k, and h is Planck’s constant. We assume that each atom can be ionised at most twice, thus the total number density of species j will be n _{iT, j} = n _{iT, j, 0} + n _{iT, j, 1} + n _{iT, j, 2}, and the number of electrons contributed from species j is n _{eT, j} = n _{iT, j, 1} + 2n _{iT, j, 2}. Since the number density of each ionisation level of each species can be written as a function of electron number density only, the total electron number density from thermal ionisation is

(15) $$\begin{equation} n_{{\rm e}_{\rm T}} = \sum _j \left[n_{{\rm i}_{\rm T},j,1}\left(n_{{\rm e}_{\rm T}} \right) + 2n_{{\rm i}_{\rm T},j,2}\left(n_{{\rm e}_{\rm T}} \right)\right], \end{equation}$$

which can then be solved iteratively using the Newton-Raphson method. The singly and doubly ionised ions, with number densities n _{iT, 1} and n _{iT, 2}, respectively, are tracked separately since they have different charges.

The average ion mass is given by Draine & Sutin (Reference Draine and Sutin1987):

(16) $$\begin{equation} m_{{\rm i}_{\rm T}} = \left(\frac{n_{{\rm i}_{\rm T},1} + n_{{\rm i}_{\rm T},2}}{\sum _{j,k} \frac{ n_{{\rm i}_{\rm T},j,k} }{ \sqrt{m_{{\rm i}_{\rm T},j}} }}\right)^2, \end{equation}$$

where m _{iT, j} is the mass of a particle of species j; note that all ionisation levels are assumed to have the same mass.

In their analytical discussion, Keith & Wardle (Reference Keith and Wardle2014) account for the possibility that dust grains absorb electrons. Similar to cosmic ray ionisation, this would lead to a reduction in the electron number density and a population of negatively charged grains. However, electron absorption by grains is not currently included in our thermal ionisation model. At high temperatures, the grains absorb less than a few percent of the electrons, thus do not contribute meaningfully to n _eT . At low temperatures, Z _iT = 0 since n _eT = 0, but numerical round-off error leads to Z _iT ≈ 0, which then gives non-sensical results of the electron number density. Given the physical irrelevance and numerical difficulties at high and low temperatures, respectively, electron absorption is not included in this version of Nicil.

2.3. Charged species populations

Both cosmic ray and thermal ionisation are calculated independently of one another for numerical stability and efficiency. This results in two disconnected free electron populations, and the total electron number density is given by n _e = n _eR + n _eT . As will be shown in Section 3, one ionisation process typically dominates the other, thus our method is reasonable.

The ionisation processes yield seven charged populations:

1. 1. free electrons with mass m _e, number density n _e, and charge Z _e ≡ −1,

2. 2. light ions from cosmic ray ionisation with mass m _iR , number density n _iR , and charge Z _iR = 1,

3. 3. metallic ions from cosmic ray ionisation with mass m _IR , number density n _IR , and charge Z _IR = 1,

4. 4. singly ionised ions from thermal ionisation with mass m _iT , number density n _{iT, 1}, and charge Z _{iT, 1} = 1,

5. 5. doubly ionised ions from thermal ionisation with mass m _iT , number density n _{iT, 2}, and charge Z _{iT, 2} = 2,

6. 6. charged grains from cosmic ray ionisation with mass m _g, number density n _g(Z = −1), and charge Z _g = −1.

7. 7. charged grains from cosmic ray ionisation with mass m _g, number density n _g(Z = 1), and charge Z _g = +1.

These species are defined with subscripts e, i_R, I_R, i_{T, 1}, i_{T, 2}, and g, respectively, where the latter will be used generally for all grain populations. The density of neutral particles is

(17) $$\begin{eqnarray} \rho _{\rm n} &=& \rho - \left(\rho _{{\rm i}_{\rm R}} + \rho _{{\rm I}_{\rm R}} + \rho _{{\rm i}_{\rm T}} + \rho _{\rm e}\right), \nonumber\\ &=& \rho - \left[n_{{\rm i}_{\rm R}}m_{{\rm i}_{\rm R}} + n_{{\rm I}_{\rm R}}m_{{\rm I}_{\rm R}} + \left(n_{{\rm i}_{\rm T},1} + n_{{\rm i}_{\rm T},2}\right)m_{{\rm i}_{\rm T}} + n_{\rm e}m_{\rm e}\right]. \end{eqnarray}$$

2.4. Conductivities

The conductivities are calculated under the assumption that the fluid evolves on a timescale longer than the collision timescale of any charged particle. The collisional frequencies, ν, are empirically calculated rates for charged species. The electron-ion rate is (Pandey & Wardle Reference Pandey and Wardle2008)

(18) $$\begin{equation} \nu _{\rm ei} = 51 \ {\rm s}^{-1} \left(\frac{n_{\rm e}}{{\rm cm}^{-3}}\right)\left(\frac{T}{{\rm K}}\right)^{-3/2}, \end{equation}$$

and the ion–electron rate is $\nu _{\rm ie} = \frac{\rho _{\rm e}}{\rho _{\rm i}}\nu _{\rm ei}$ . Note that ν_eiR = ν_eIR = ν_{eiT, 1} = ν_{eiT, 2} ≡ ν_ei because there is no dependence on ion properties, and ν_iRe ≠ ν_iTe.

The plasma-neutral collisional frequency is

(19) $$\begin{equation} \nu _{j{\rm n}} = \frac{\left< \sigma v\right>_{j{\rm n}}}{m_{\rm n}+m_j}\rho _{\rm n}, \end{equation}$$

where < σv > _jn is the rate coefficient for the momentum transfer by the collision of particle of species j with the neutrals. For electron-neutral collisions, it is assumed that the neutrals are comprised of hydrogen and helium, since these two elements should dominate any physical system. The rate coefficient is then given by

(20) $$\begin{equation} \left< \sigma v\right>_{\rm en} = X_{{\rm H}_2} \left< \sigma v\right>_{{\rm e-H}_2} + X_{\rm H} \left< \sigma v\right>_{{\rm e-H}} + Y \left< \sigma v\right>_{\rm e-He}, \end{equation}$$

where X _H2 , X _H, and Y are the mass fractions of molecular and atomic hydrogen and helium, respectively, and X _H2 + X _H ≡ X; X and Y are free parameters if use_massfrac=.true., otherwise they are calculated from the given abundances. Following Pinto & Galli (Reference Pinto and Galli2008),

(21) $$\begin{eqnarray} \left< \sigma v\right>_{{\rm e-H}_2} &=& \left(0.535+0.203\theta -0.163\theta ^2+0.050\theta ^3\right) \times 10^{-9} \nonumber \\ && \times\, T^{1/2} \ {\rm cm}^3 \ {\rm s}^{-1} , \end{eqnarray}$$

(22) $$\begin{eqnarray} \left< \sigma v\right>_{{\rm e-H}} &=& \left(2.841+0.093\theta -0.245\theta ^2+0.089\theta ^3\right) \times 10^{-9} \nonumber \\ && \times\, T^{1/2} \ {\rm cm}^3 \ {\rm s}^{-1} , \end{eqnarray}$$

(23) $$\begin{eqnarray} \left< \sigma v\right>_{\rm e-He} &=& 0.428 \times 10^{-9} T^{1/2} \ {\rm cm}^3 \ {\rm s}^{-1}, \end{eqnarray}$$

where θ ≡ log T for T in K and the drift velocity is assumed to be zero.

The ion-neutral rate is (Pinto & Galli Reference Pinto and Galli2008)

(24) $$\begin{eqnarray} \left< \sigma v\right>_{\rm in} &=& 2.81\times 10^{-9} \ {\rm cm}^3 \ {\rm s}^{-1} \nonumber \\ &&\times \left|Z_{\rm i}\right|^{1/2}\left[ X_{{\rm H}_2}\left(\frac{p_{{\rm H}_2}}{{\rm \AA}^3}\right)^{1/2}\left(\frac{\mu _{{\rm i-H}_2}}{m_{\rm p}}\right)^{-1/2} \right. \nonumber \\ &&+ \left. X_{{\rm H}}\left(\frac{p_{{\rm H}}}{{\rm \AA}^3}\right)^{1/2}\left(\frac{\mu _{{\rm i-H}}}{m_{\rm p}}\right)^{-1/2} \right. \nonumber \\ &&+ \left. Y \left(\frac{p_{\rm He}}{{\rm \AA}^3}\right)^{1/2}\left(\frac{\mu _{\rm i-He}}{m_{\rm p}}\right)^{-1/2} \right], \end{eqnarray}$$

where p _H2 , p _H, and p _He are the polarisabilities (Osterbrock Reference Osterbrock1961). This term is calculated independently for each ion due to the ion mass dependence in the reduced masses, μ_i-H2 , μ_i-H, and μ_i-He.

For grain-neutral collisions, the rate coefficient is (Wardle & Ng Reference Wardle and Ng1999; Pinto & Galli Reference Pinto and Galli2008)

(25) $$\begin{equation} \left< \sigma v\right>_{\rm gn} = \pi a_{\rm g}^2\delta _{\rm gn}\sqrt{\frac{128k_{\rm B}T}{9\pi m_{\rm n}}}, \end{equation}$$

where δ_gn is the Epstein coefficient, and

(26) $$\begin{equation} m_{\rm n} = \left(\frac{X}{m_{{\rm H}_2}}+ \frac{Y}{m_{\rm He}} + \sum _k \frac{Z_k}{m_k}\right)^{-1} \end{equation}$$

is the mass of the neutral particle, where we sum over all included metals, k (if use_massfrac=.false.), which have mass fractions Z_k and masses m_k . This is the only instance in this paper where Z does not refer to charge. Since m _n is a characteristic mass, we assume that all the hydrogen is molecular hydrogen, which is reasonable for low temperatures.

The Hall parameter for species j is given by

(27) $$\begin{equation} \beta _j = \frac{|Z_j|eB}{m_j c}\frac{1}{\nu _{j{\rm n}}}. \end{equation}$$

A modified form is presented in Appendix E.

Finally, the Ohmic, Hall, and Pedersen conductivities are given by (e.g. Wardle & Ng Reference Wardle and Ng1999; Wardle Reference Wardle2007)

(28) $$\begin{eqnarray} \sigma _{\rm O} &=& \frac{ec}{B}\sum _j n_j|Z_j|\beta _j, \end{eqnarray}$$

(29) $$\begin{eqnarray} \sigma _{\rm H} &=& \frac{ec}{B}\sum _j \frac{n_j Z_j}{1+ \beta _j^2}, \end{eqnarray}$$

(30) $$\begin{eqnarray} \sigma _{\rm P} &=& \frac{ec}{B}\sum _j \frac{n_j |Z_j| \beta _j}{1+ \beta _j^2}. \end{eqnarray}$$

We explicitly note that σ_O and σ_P are positive, whereas σ_H can be positive or negative. The total conductivity perpendicular to the magnetic field is

(31) $$\begin{equation} \sigma _\bot = \sqrt{\sigma _{\rm H}^2 + \sigma _{\rm P}^2}. \end{equation}$$

2.5. Non-ideal MHD coefficients

The induction equation in MHD is

(32) $$\begin{equation} \frac{{\rm d} \bm {B}}{{\rm d} t} = \left(\bm {B}\cdot \bm {\nabla }\right)\bm {v}-\bm {B}\left(\bm {\nabla }\cdot \bm {v}\right) + \left.\frac{{\rm d} \bm {B}}{{\rm d} t}\right|_{\rm non-ideal}, \end{equation}$$

where $ \left.\frac{{\rm d} \bm {B}}{{\rm d} t}\right|_{\rm non-ideal}$ is the non-ideal MHD term, which is the sum of Ohmic resistivity (OR), the Hall effect (HE), and ambipolar diffusion (AD) terms; these terms are given by

(33) $$\begin{eqnarray} \left.\frac{{\rm d} \bm {B}}{{\rm d} t}\right|_{\rm OR} &=& -\bm {\nabla } \times \left[ \eta _{\rm OR} \left(\bm {\nabla }\times \bm {B}\right)\right], \end{eqnarray}$$

(34) $$\begin{eqnarray} \left.\frac{{\rm d} \bm {B}}{{\rm d} t}\right|_{\rm HE} &=& -\bm {\nabla } \times \left[ \eta _{\rm HE} \left(\bm {\nabla }\times \bm {B}\right)\times \bm {\hat{B}}\right], \end{eqnarray}$$

(35) $$\begin{eqnarray} \left.\frac{{\rm d} \bm {B}}{{\rm d} t}\right|_{\rm AD} &=& \bm {\nabla } \times \left\lbrace \eta _{\rm AD}\left[\left(\bm {\nabla }\times \bm {B}\right)\times \bm {\hat{B}}\right]\times \bm {\hat{B}}\right\rbrace , \end{eqnarray}$$

where the use of the magnetic unit vector, $\bm {\hat{B}}$ , is to ensure that all three coefficients, η, have units of area per time.

From conservation of energy, the contribution to internal energy, u, from the non-ideal MHD terms (Wurster, Price, & Ayliffe Reference Wurster, Price and Ayliffe2014) is

(36) $$\begin{eqnarray} \left.\frac{{\rm d} u}{{\rm d} t}\right|_{\rm OR} &=& - \frac{\eta _{\rm OR}}{\rho }\bm {J}\cdot \bm {J} = - \frac{\eta _{\rm OR}}{\rho } J^2, \end{eqnarray}$$

(37) $$\begin{eqnarray} \left.\frac{{\rm d} u}{{\rm d} t}\right|_{\rm HE} &=& - \frac{\eta _{\rm HE}}{\rho }\left(\bm {J}\times \hat{\bm {B}}\right)\cdot \bm {J} = 0, \end{eqnarray}$$

(38) $$\begin{eqnarray} \left.\frac{{\rm d} u}{{\rm d} t}\right|_{\rm AD} &=& - \frac{\eta _{\rm AD}}{\rho }\left[\left(\bm {J}\times \hat{\bm {B}}\right)\times \hat{\bm {B}}\right]\cdot \bm {J} \nonumber\\ &=& - \frac{\eta _{\rm AD}}{\rho }\left[\left(\bm {J}\cdot \hat{\bm {B}}\right)^2 - J^2\right]. \end{eqnarray}$$

Next, we invoke the strong coupling approximation, which allows the medium to be treated using the single fluid approximation. It is assumed that the ion pressure and momentum are negligible compared to that of the neutrals, i.e. ρ ~ ρ_n and ρ_i ≪ ρ, where ρ, ρ_n and ρ_i are the total, neutral, and ion mass densities, respectively.

The general form of the coefficients (Wardle Reference Wardle2007) is

(39) $$\begin{eqnarray} \eta _{\rm OR} &=& \frac{c^2}{4\pi \sigma _{\rm O}}, \end{eqnarray}$$

(40) $$\begin{eqnarray} \eta _{\rm HE} &=& \frac{c^2}{4\pi \sigma _\bot }\frac{\sigma _{\rm H}}{\sigma _\bot }, \end{eqnarray}$$

(41) $$\begin{eqnarray} \eta _{\rm AD} &=& \frac{c^2}{4\pi \sigma _\bot }\frac{\sigma _{\rm P}}{\sigma _\bot } - \eta _{\rm OR} \nonumber \\ &=& \frac{c^2}{4\pi \sigma _{\rm O}}\frac{\sigma _{\rm O}\sigma _{\rm P} - \sigma _\bot ^2}{\sigma _\bot ^2} \equiv \frac{c^2}{4\pi \sigma _{\rm O}}\frac{\sigma _{\rm A}}{\sigma _\bot ^2}. \end{eqnarray}$$

The value of η depends on the microphysics of the model, and η_HE can be either positive or negative, whereas η_OR and η_AD are positive (e.g. Wardle & Ng Reference Wardle and Ng1999). To prevent η_AD ≲ 0 due to numerical round-off when σ_Oσ_P ≈ σ² _⊥, we use

(42) $$\begin{eqnarray} \sigma _{\rm O}\sigma _{\rm P} - \sigma _\bot ^2 &\equiv & \sigma _{\rm A} \nonumber \\ &=& \left(\frac{ec}{B}\right)^2\sum _{j>k} \left[ \frac{ n_j|Z_j|\beta _j}{1+\beta _j^2}\frac{n_k|Z_k|\beta _k}{1+\beta _k^2}\right. \nonumber \\ &&\times \, \left.\left(\frac{Z_k\beta _k}{|Z_k|} - \frac{Z_j\beta _j}{|Z_j|}\right)^2\right], \end{eqnarray}$$

where all pairs of charged species, j and k, are summed over (M. Wardle, private communication 2014).

The non-ideal MHD terms constrain the timestep by

(43) $$\begin{equation} {\rm d}t < C_{\rm non-ideal}\frac{h^2}{|\eta |}, \end{equation}$$

where h is the particle smoothing length or the cell size, η = max (η_OR, |η_HE|, η_AD) and C _non-ideal < 1 is a positive coefficient analogous to the Courant number. Stability tests by Bai (Reference Bai2014) suggest $C_{\rm non\hbox{-}ideal} \lesssim \frac{1}{6}$ is required, whilst wave calculations and tests by Wurster et al. (Reference Wurster, Price and Bate2016) suggest $C_{\rm non\hbox{-}ideal} \lesssim \frac{1}{2\pi }$ is required for stability. Thus, stability can be achieved by choosing a value of C _non-ideal appropriate for a conditionally stable algorithm.

5.1. Parameter modifications to nicil.F90

The parameters which can be modified by the user are listed at the top of nicil.F90 between the lines labelled ‘Input Parameters,’ and ‘End of Input Parameters.’ These parameters are defined as public, thus can also be modified by other programmes in the user’s code.

If the user wishes to add or remove elements if use_massfrac=.false., then this is done in the subroutine nicil_initial_species. Modify nelements to be the new number of elements, and add the characteristics of the new elements to the relevant arrays, using the current elements as a template.

5.2. Modifications to user’s initialisation subroutine

Nicil contains several variables that require initialising. When defining the variables and subroutines in the existing initialisation subroutine, include

• use nicil, only : nicil_initialise

• integer :: ierr .

Once the units are initialised by the existing code, include

• call nicil_initialise(utime,udist, &umass, unit_Bfield,ierr)

• if (ierr > 0) & call abort_subroutine_of_user

where utime, udist, umass, and unit_Bfield are unit conversions from CGS to code units. Note that the code unit of temperature is assumed to be Kelvin. If ierr > 0, then an error occurred during initialisation, and the user’s programme should be immediately aborted, using the user’s preferred method (e.g. by calling abort_subroutine_of_user).

The number densities are calculated iteratively, thus the user must define and save these arrays within the body of their code; the number density array for cosmic rays (thermal ionisation) requires four values (one value) for every cell/particle. The user is required to initialise these arrays to zero. For simplicity, they can be treated (i.e. saved and updated) analogously to the user’s existing density array.

5.3. Modifications to user’s density loop

In many codes, the density is calculated prior to calculating the forces and magnetic fields. At the same time, Nicil should calculate the number densities, which do not depend on any neighbours. To calculate the ith element of these arrays, nR and neT, include

• call nicil_get_ion_n(rho(i),T(i),&

• nR(1:4,i),neT(i),ierr)

• if (ierr/=0) then

• call nicil_translate_error(ierr)

• if (ierr > 0) call abort_subroutine_of_user

• end if

in the density loop, where rho(i) and T(i) are the current values of density and temperature, respectively, and nR(1:4,i) and neT(i) are the characteristic number densities from cosmic rays and the electron number density from thermal ionisation, respectively. If an error has occurred, nicil_get_ion_n will return ierr ≠ 0, which will then be passed into nicil_translate_error. If ierr is returned from nicil_translate_error as a negative number, then the error is not fatal; the error message will be written to file if warn_verbose = .true., otherwise all non-fatal errors will be suppressed. If a fatal error message is returned, then it will be printed to file prior to the user executing the existing code’s subroutine to terminate the programme.

The following subroutines and variables need to be imported and defined:

• use nicil, only : nicil_get_ion_n, &

• nicil_translate_error

• integer :: ierr .

5.4. Modifications to user’s magnetic field loop

In the subroutine where the magnetic field is calculated, define

• use nicil, only : nicil_get_eta, &

• nicil_translate_error

• integer:: ierr

• real:: eta_ohmi, eta_halli, eta_ambii

Then, for each particle/cell i, include

• call nicil_get_eta(eta_ohmi, &

• eta_halli,eta_ambii,Bi,rho(i), &T(i),nR(1:4,i),neT(i),ierr)

• if (ierr/=0) then

• call nicil_translate_error(ierr)

• if (ierr > 0) & call abort_subroutine_of_user

• end if

where eta_ohmi, eta_halli, and eta_ambii are the calculated (i.e. output) non-ideal MHD coefficients, and Bi is the magnitude of the particle/cell’s magnetic field.

With the non-ideal MHD coefficients now calculated, the user can calculate $\left.\frac{{\rm d} \bm {B}}{{\rm d} t}\right|_{\rm non\hbox{-}ideal}$ using the existing routines within the user’s code.

5.5. Optional modifications to the user’s code

5.5.3. Ion velocity

The ion velocity is given by

(47) $$\begin{equation} \bm {v}_{\rm ion} = \bm {v} + \eta _{\rm AD}\left(\bm {\nabla }\times \hat{B}\right)\times \hat{B}, \end{equation}$$

which can optionally be calculated by calling

• call nicil_get_vion(eta_ambii,&

• vxi,vyi,vzi,Bxi,Byi,Bzi,jcurrenti,&

• vioni,ierr)

where vxi, vyi, vzi, Bxi, Byi, Bzi are the components of the velocity and magnetic field of particle/cell i, and vioni is an output array containing the three components of the ion velocity. Similar to nicil_get_dt, it is optimal to call this array in the i-loop after the non-ideal MHD coefficients are calculated for particle i. Note that ierr is not initialised in this subroutine, thus must either be passed in from nicil_get_eta, or reinitialised before calling this routine. The drift velocity is calculated by

(48) $$\begin{equation} \bm {v}_{\rm drift} = \bm {v} - \bm {v}_{\rm ion}. \end{equation}$$

If warn_verbose = .true., then a warning will be printed if the ion and neutral velocity differ by more than 10%; recall that v _drift ≈ 0 is assumed in rate coefficients given in (21)–(23).