1 INTRODUCTION
Astrophysical magnetic fields are often modelled using ideal magnetohydrodynamics (MHD), which assumes that the gas is fully ionised and that the electrons are tied to the magnetic field and have zero resistivity. However, in molecular cloud cores, for example, this has been known to be a poor approximation to the true conditions since at least Mestel & Spitzer (Reference Mestel and Spitzer1956), where the ionisation levels can be as low as $n_{\text{e}}/ n_{\text{H}_{2}} = 10^{14}$ (Nakano & Umebayashi Reference Nakano and Umebayashi1986; Umebayashi & Nakano Reference Umebayashi and Nakano1990).
When a medium is only partially ionised, such as in a molecular cloud, three nonideal MHD effects occur:

1. Ohmic resistivity: the drift between electrons and ions/neutrals; neither ions nor electrons are tied to the magnetic field,

2. Hall effect: ionelectron drift; only electrons are tied to the magnetic field, and

3. Ambipolar diffusion: ionneutral drift; both ions and electrons are tied to the magnetic field.
These environmentally dependent terms add resistivity to the electrons, and ultimately affect the evolution of the magnetic fiels (e.g. Wardle & Ng Reference Wardle and Ng1999; Nakano, Nishi, & Umebayashi Reference Nakano, Nishi and Umebayashi2002; Tassis & Mouschovias Reference Tassis and Mouschovias2007; Wardle Reference Wardle2007; Pandey & Wardle Reference Pandey and Wardle2008; Keith & Wardle Reference Keith and Wardle2014). Recent studies, both idealised (e.g. Bai Reference Bai2014, Reference Bai2015) and realistic (e.g. Shu et al. Reference Shu, Galli, Lizano and Cai2006; Mellon & Li Reference Mellon and Li2009; Duffin & Pudritz Reference Duffin and Pudritz2009; Dapp & Basu Reference Dapp and Basu2010; Krasnopolsky, Li, & Shang Reference Krasnopolsky, Li and Shang2010; Machida, Inutsuka, & Matsumoto Reference Machida, Inutsuka and Matsumoto2011; Li, Krasnopolsky, & Shang Reference Li, Krasnopolsky and Shang2011; Dapp, Basu, & Kunz Reference Dapp, Basu and Kunz2012; Tomida et al. Reference Tomida, Tomisaka, Matsumoto, Hori, Okuzumi, Machida and Saigo2013; Tomida, Okuzumi, & Machida Reference Tomida, Okuzumi and Machida2015; Tsukamoto et al. Reference Tsukamoto, Iwasaki, Okuzumi, Machida and Inutsuka2015a, Reference Tsukamoto, Iwasaki, Okuzumi, Machida and Inutsuka2015b; Wurster, Price, & Bate Reference Wurster, Price and Bate2016) have included some or all of the nonideal MHD effects. In each study, however, the authors have been required to write and implement their own algorithms to model the nonideal MHD effects.
Multiple algorithms existing within a community provides a selfconsistent check on the processes, and this helps to validate the results. However, writing and implementing the nonideal MHD algorithms is a nontrivial process, both physically and numerically. Thus, we introduce Nicil: NonIdeal MHD Coefficients and Ionisation Library, which is a prewritten, Fortran90 module that selfconsistently calculates the nonideal MHD coefficients and is ready for download and implementation into existing codes. Early version of Nicil are used in Wurster et al. (Reference Wurster, Price and Bate2016) and Wurster, Price, & Bate (in preparation). We are aware that a similar standalone code has recently been published by Marchand et al. (Reference Marchand, Masson, Chabrier, Hennebelle, Commerçon and Vaytet2016).
This paper is a usermanual for Nicil, and is organised as follows. Section 2 describes the algorithms that Nicil uses. We perform a simple test of Nicil in Section 3, and provide an astrophysical example in Section 4. Section 5 explains how to implement Nicil into an existing numerical code, and we conclude in Section 6. Throughout the text, we will refer to the ‘author’ as the author of Nicil, the ‘user’ as the person who is implementing Nicil into the existing code, and the ‘user’s code’ or ‘existing code’ as the code into which Nicil is being implemented.
2 ALGORITHMS
The complete Nicil library can be downloaded at www.bitbucket.org/jameswurster/nicil. This is a free library under the GNU license agreement: free to use, modify, and share, does not come with a warranty, and this paper must be cited if Nicil or any modified version thereof is used in a study. The important files are listed and summarised in Appendix A. The important parameters which the user can modify, along with their default values, are listed in Appendix B.
Nicil includes various processes that can be included or excluded at the user’s discretion. Major processes are cosmic ray and thermal ionisation (one or both can be selected) and the assumption of the grain size distribution (fixed or MRN distribution; Mathis, Rumpl, & Nordsieck Reference Mathis, Rumpl and Nordsieck1977). By default, all three nonideal MHD coefficients are selfconsistently calculated, however, Nicil includes the option of calculating only selected terms or using fixed coefficients (see Appendix C for a discussion of fixed coefficients).
The remainder of this section describes the processes that are used to calculate the nonideal MHD coefficients. We define all variables, but for clarity, we do not state the default values in the main text.
2.1. Cosmic ray ionisation
When a cosmic ray (or an Xray) 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),
where $f_{\text{dg}}$ is the dusttogas 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.,
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^{−25} 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 logaverage grain size, and the number density of each bin is calculated using (2).
The mass of a grain particle is given by
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)
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)
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)
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,
and from conservation of grains,
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 oversubscription, 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
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,
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
which can then be solved iteratively using the NewtonRaphson 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):
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 roundoff error leads to Z _{iT } ≈ 0, which then gives nonsensical 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. free electrons with mass m _{e}, number density n _{e}, and charge Z _{e} ≡ −1,

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

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

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

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

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

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
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 electronion rate is (Pandey & Wardle Reference Pandey and Wardle2008)
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 plasmaneutral collisional frequency is
where < σv > _{ jn} is the rate coefficient for the momentum transfer by the collision of particle of species j with the neutrals. For electronneutral 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
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),
where θ ≡ log T for T in K and the drift velocity is assumed to be zero.
The ionneutral rate is (Pinto & Galli Reference Pinto and Galli2008)
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, μ_{iH2 }, μ_{iH}, and μ_{iHe}.
For grainneutral collisions, the rate coefficient is (Wardle & Ng Reference Wardle and Ng1999; Pinto & Galli Reference Pinto and Galli2008)
where δ_{gn} is the Epstein coefficient, and
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
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)
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
2.5. Nonideal MHD coefficients
The induction equation in MHD is
where $ \left.\frac{{\rm d} \bm {B}}{{\rm d} t}\right_{\rm nonideal}$ is the nonideal MHD term, which is the sum of Ohmic resistivity (OR), the Hall effect (HE), and ambipolar diffusion (AD) terms; these terms are given by
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 nonideal MHD terms (Wurster, Price, & Ayliffe Reference Wurster, Price and Ayliffe2014) is
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
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 roundoff when σ_{O}σ_{P} ≈ σ^{2} _{⊥}, we use
where all pairs of charged species, j and k, are summed over (M. Wardle, private communication 2014).
The nonideal MHD terms constrain the timestep by
where h is the particle smoothing length or the cell size, η = max (η_{OR}, η_{HE}, η_{AD}) and C _{nonideal} < 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 _{nonideal} appropriate for a conditionally stable algorithm.
3 TESTS PROGRAMMES AND RESULTS
The Nicil library includes two standalone test codes, which can be used to test the effects of varying the parameters. Both can be compiled by typing make in the NICIL directory. This yields the executables nicil_ex_eta and nicil_ex_sph. The latter is a simple SPH code that is primarily used to demonstrate how Nicil is to be implemented into an existing code.
The test programme nicil_ex_eta will calculate the nonideal MHD coefficients, and output them and the constituent components that are required for their calculation, including number densities of all the charged species and elements, the conductivities, and the coefficients. By default, it outputs three sets of data: one over a range of densities assuming a constant temperature, one over a range of temperatures assuming a constant density, and one where temperature and density are related through the barotropic equation of state (Machida, Inutsuka, & Matsumoto Reference Machida, Inutsuka and Matsumoto2006):
where T _{0} = 10 K, n is the total density, n _{1} = 10^{11}, n _{2} = 10^{16}, and n _{3} = 10^{21} cm^{−3}, Γ_{1} = 0.4, Γ_{2} = −0.3, and Γ_{3} = 0.56667; the fixed temperature and density for the first two data sets can be modified in src/ nicil_ex_eta.F90. When a magnetic field is required for the constant density and temperature plots, for the purpose of illustration, we use the relation used in Wardle & Ng (Reference Wardle and Ng1999):
When a magnetic field is required for the barotropic equation of state, we use (Li et al. Reference Li, Krasnopolsky and Shang2011)
Figure 1 shows the output of nicil_ex_eta using the default options listed in Table A2; the lefthand plot shows the results using T = 30 K and T ≡ T(n) as a function of density, and the righthand plot shows the results using ρ = 10^{−13} g cm^{−3} and n ≡ n(T) as a function of temperature.
For constant T = 30 K, the trends for the ion and electron number densities and the conductivities are in approximate agreement with Wardle & Ng (Reference Wardle and Ng1999); differences arise as a result of different initial conditions and assumptions. The elemental number densities are zero, as is expected at such a low temperature. At this temperature, the conductivities and nonideal MHD coefficients are strongly dependent on density, whilst the grain charge and ion and electron number densities essentially have two states: one at high density and one at low density.
For constant ρ = 10^{−13} g cm^{−3}, which is characteristic in discs around protostars, cosmic ray ionisation is the dominant ionisation source for T ≲ 600 K; for T ≳ 600 K, thermal ionisation is the dominant source of electrons. The changeover from the dominance of cosmic ray ionisation to thermal ionisation is abrupt, confirming that the two processes can be calculated independently without loss of information.
These tests were performed using Version 1.2.1 of Nicil, which can be found in commit 8220ec8 from 2016 August 2. The plots can be reproduced using the included Python graphing script, plot_results.py, which calls GNUplot. It is recommended that the user becomes familiar with the effect of the various parameters and processes prior to implementing Nicil into the existing code. This can be done by modifying the parameters (see Section 5.1), plotting them using plot_results.py, and then comparing the new graphs to those made with the default parameters.
4 TEST IN A COLLAPSING MOLECULAR CLOUD
For an astrophysical example, we model the collapse of a rotating, 1 M_{⊙} cloud of gas using the 3D smoothed particle MHD code Phantom with the inclusion of selfgravity. The cloud has an initial radius of 4 × 10^{16} cm and is initially threaded with a magnetic field of strength 163 μG (five critical masstoflux units) which is antialigned with the angular momentum vector. When the maximum density surpasses ρ_{crit} = 10^{−10} g cm^{−3} and a given set of criteria are met, the densest gas particle is replaced with a sink particle and its neighbours within r _{acc} = 6.7 AU are immediately accreted onto it. This setup is the same as in Wurster et al. (Reference Wurster, Price and Bate2016), with 10^{6} particles in the cloud.
We ran four tests: Ideal MHD, Nicil with the default settings, using MRN grain size distribution with five bins, and using only cosmic ray ionisation. The faceon gas column density at 1.07 t _{ff} (where the freefall time is t _{ff} = 2.4 × 10^{4} yr) is shown in Figure 2 for the ideal, default, and MRN models. No discernible disc forms in the ideal MHD model, which demonstrates the magnetic braking catastrophe (e.g. Allen, Li, & Shu (Reference Allen, Li and Shu2003); Price & Bate (Reference Price and Bate2007); Mellon & Li (Reference Mellon and Li2008); Hennebelle & Fromang (Reference Hennebelle and Fromang2008); Wurster et al. Reference Wurster, Price and Bate2016). When the default parameters are used, a weak disc forms around the sink particle. A lowmass disc may also be forming in the MRN model, however, it is difficult to draw conclusions since such a disc will be heavily influenced by the sink particle. At this time, the temperature is less than a few hundred Kelvin, thus thermal ionisation plays a minimal role. As expected, the cosmic rayonly model yields similar results to the default model.
The radial profiles in the disc of the magnetic field, temperature, and nonideal MHD coefficients are shown in Figure 3 for the ideal MHD, default, and MRN models. Gas is assumed to be ‘in the disc’, if ρ > ρ_{disc, min} = 10^{−13} g cm^{−3}.
Given that the maximum temperature is typically ≲ 200 K, thermal ionisation plays a negligible role. The coefficients are typically larger in magnitude for the default model, with η_{HE} > η_{AD} > η_{OR}. The negative Hall coefficient is responsible for increasing the magnetic field from the ideal MHD to the MRN to the default model.
This MRN grain model includes five bins, and takes 55% longer to reach the t = 1.07 t _{ff} than the default model. Thus, the user must be cautious about performance when using the MRN model.
5 IMPLEMENTATION INTO EXISTING CODE
Nicil is a Fortran90 module that is contained in one file, src/ nicil.F90. To include Nicil in the user’s code, copy this file into the existing code’s source directory and include it in the Makefile, ensuring that nicil.F90 is compiled in double precision; see src/Makefile for details.
The modifications required to the user’s code are summarised below. The Nicil library includes a simple SPH programme, src/ nicil_ex_sph.F90, which can be used as an example.
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 nonfatal 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) nonideal MHD coefficients, and Bi is the magnitude of the particle/cell’s magnetic field.
With the nonideal 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.1. Timesteps
For each particle/cell, the nonideal MHD timestep can be calculated by calling

call nimhd_get_dt(dtohmi,dthalli, &

dtambii,hi, &eta_ohmi,eta_halli,eta_ambii)
where hi is the smoothing length of the particle or the size of the cell, and dtohmi, dthalli, and dtambii are the three output timesteps. These are returned independently in case the user employs a timestepping scheme such as super timestepping (Alexiades, Amiez, & Gremaud Reference Alexiades, Amiez and Gremaud1996), which requires dtohmi and dtambii to be treated differently from dthalli. The timesteps can be calculated with minimal extra cost if the modifications are included in the iloop after the nonideal MHD coefficients are already calculated.
5.5.2. Internal energy
Nicil also includes a subroutine that evolves the internal energy of the particle or cell. This is called by

call nimhd_get_dudt(dudtnonideal, &

eta_ohmi,eta_ambii,rhoi,jcurrenti, &Bxyzi)
where Bxyzi is an array containing the three components of the magnetic field, i.e. (Bxi,Byi,Bzi), and dudtnonideal is the calculated (i.e. output) change in energy for the particle/cell caused by the nonideal MHD effects.
5.5.3. Ion velocity
The ion velocity is given by
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 iloop after the nonideal 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
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).
5.5.4. Removing ifstatements
Given the versatile design of Nicil, the source code contains several ifstatements whose argument is a logical operator set prior to runtime. Since these ifstatements will be called repeatedly throughout a calculation, their calls may decrease performance of the code if not properly optimised by the compiler. Thus, src/hardcode_ifs.py can be executed in the src/ directory to remove them.
Specifically, this script will read through nicil_source.F90, determine the value of the logical operators, then rewrite nicil.F90 such that every time an ifstatement is encountered with a logical operator, the ifstatement will be removed, as will its contents if the logical is set to .false.. Note that any changes made to nicil.F90 will be overwritten, thus all parameters must be set in nicil_source.F90 if this script is to be used. The files nicil.F90 and nicil_source.F90 are identical in the repository.
6 CONCLUSION
We have introduced NICIL: NonIdeal MHD Coefficients and Ionisation Library. This a Fortran90 module that is fully parameterisable to allow the user to determine which processes to include and the values of the parameters. We have described its algorithms, including the cosmic ray and thermal ionisation processes, and how the conductivities and nonideal MHD coefficients are calculated. We have summarised how to implement Nicil into an existing code.
The library contains two test codes, one of which outputs the nonideal MHD coefficients and its constituent components. We have described the results using the default values, and showed how the components are affected by temperature and density.
Nicil is a selfcontained code that is freely available, and is ready to be implemented in existing codes.
ACKNOWLEDGEMENTS
We would like to thank Daniel Price and Mathew Bate for useful discussions and for their testing and debugging efforts. JW acknowledges support from the Australian Research Council (ARC) Discovery Projects Grant DP130102078 and from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007–2013 grant agreement no. 339248).
A LIST OF FILES
Table A1 lists and summarises the important files that are included in the Nicil library.
B FREE PARAMETERS
Table A2 lists the free parameters, the default settings, and the reference for the default value (where applicable).
For use_massfrac = .false., abundances are used in the thermal ionisation calculation. Five elements are included by default, and they, their abundances, and their first and second ionisation potentials are listed in Table A3. The abundance, x_{j} is related to the logarithmic abundance X_{j} by x_{j} = 10^{ Xj }/(∑_{ i }10^{ Xi }).
We initially assume that all of the hydrogen is molecular hydrogen, H_{2}. Using the local temperature and the Saha equation, the number densities of both molecular and atomic hydrogen are determined, viz.,
where χ_{ H − H 2 } = 4.476 eV is the dissociation potential and n _{H} + n _{H2 } = A _{H} n _{0}/2, where A _{H} is the abundance of hydrogen as given in Table A3 and n _{0} is the total number density. The mass fractions, X _{H2 } and X _{H}, are then calculated for use in the rate coefficients. For thermal ionisation, we only allow molecular hydrogen to be singly ionised, with an ionisation potential of 15.60 eV.
C CONSTANT NONIDEAL MHD COEFFICIENTS
In some scenarios, such as test problems, it may be useful to use semiconstant coefficients. To implement this, set eta_constant=.true.; the coefficients then become
where B is the magnitude of the magnetic field and $v_{\rm A} = \frac{B}{\sqrt{4\pi \rho }}$ is the Alfvén velocity which are selfconsistently calculated, as done during tests of ambipolar diffusion in e.g. Mac Low et al. (Reference Low, Norman, Konigl and Wardle1995), Choi, Kim, & Wiita (Reference Choi, Kim and Wiita2009), and Wurster et al. (Reference Wurster, Price and Ayliffe2014). C is a free parameter to be input by the user. The default values are given in Table A4.
To calculate the coefficients from fixed parameters that are characteristic of the environment, set ${\tt eta\_const\_calc = .true.}$ , and the coefficients become
where the free parameters are the fixed electron number density n _{e, 0}, the fixed ion density ρ_{i, 0}, the fixed neutral density ρ_{n, 0}, the powerlaw exponent α (α = 0 for molecular cloud densities, and α = 0.5 for lowdensity regime; c.f. Mac Low et al. Reference Low, Norman, Konigl and Wardle1995), the collisional coupling coefficient for ambipolar diffusion γ_{AD} and s _{H} = ±1 is the sign of the Hall coefficient. Invoking the strong coupling approximation, we set ρ_{n} = ρ in (A7). The default values are given in Table A4.
Since these forms of the nonideal MHD coefficients are not selfconsistently calculated, it is recommended that eta_constant=.true. be used for testing purposes only.
D APPROXIMATING NUMBER DENSITIES USING AVERAGE GRAIN CHARGES
If NICIL fails to calculate the number densities from cosmic rays using the Jacobian method describe in Section 2.1, then the following method is invoked to approximate the number densities. This method is used with caution since it makes several assumptions. First, it is assumed that n _{n} ≈ n. Second, it is assumed that recombination rate is small compared to the other processes, thus k _{es } = 0. Third, the grains of each size are treated as a single population with an average grain charge of $\bar{Z} < 0$ . From (4) and (5), these simplifications allow the ion and electron number densities to be written as
where s ∈ {i_{R}, I_{R}}. Then, from charge neutrality,
the NewtonRaphson method can be used to solve for $\bar{Z}$ . Recall that n _{g}(a_{j} ) is the total grain number density of grains of size a_{j} .
Using $\bar{Z}$ , the ion and electron number densities are calculated using (A8) and (A9), respectively. Using (6), charge neutrality and grain conservation, the grain number densities are calculated viz.,
where k ^{+} ≡ k(Z = 1, a_{j} ), k ^{−} ≡ k(Z = −1, a_{j} ) and k ^{0} ≡ k(Z = 0, a_{j} ).
E MODIFIED HALL PARAMETER
A modified version of the Hall parameter was used in Wurster et al. (Reference Wurster, Price and Bate2016):
With these modifications, the coefficient for Ohmic resistivity from Pandey & Wardle (Reference Pandey and Wardle2008) and Keith & Wardle (Reference Keith and Wardle2014) can be recovered under the assumption β_{i} ≪ β_{e}. By default, these forms are not used, but are included for legacy reasons.