The challenges of modelling microphysics: ambipolar diffusion, chemistry, and cosmic rays in MHD shocks

TL;DR

This paper develops a 1D non-ideal MHD code to study the impact of microphysics, including cosmic rays and dust, on magnetic shock waves in astrophysical environments, highlighting their critical role.

Contribution

It introduces a new non-ideal MHD simulation framework that incorporates complex microphysical effects to better understand astrophysical shock phenomena.

Findings

Cosmic rays significantly influence shock dynamics.

Dust microphysics alters shock structure.

Microphysics inclusion is essential for realistic models.

Abstract

From molecular clouds to protoplanetary disks, non-ideal magnetic effects are important in many astrophysical environments. Indeed, in star and disk formation processes, it has become clear that these effects are critical to the evolution of the system. The efficacy of non-ideal effects are, however, determined by the complex interplay between magnetic fields, ionising radiation, cosmic rays, microphysics, and chemistry. In order to understand these key microphysical parameters, we present a one-dimensional non-ideal magnetohydrodynamics code and apply it to a model of a time-dependent, oblique, magnetic shock wave. By varying the microphysical ingredients of the model, we find that cosmic rays and dust play a major role, and that, despite the uncertainties, the inclusion of microphysics is essential to obtain a realistic outcome in magnetic astrophysical simulations.

Tables

Table 1. Table 1: Coefficients used in Eq. ( 25 ) to fit the amount of heat deposited per cosmic ray ionisation following. Note that units are eV.

$i$	$c_i$	$i$	$c_i$	$i$	$c_i$
0	6.882876	2	-0.532834	4	-0.016907
1	2.231421	3	0.146966	5	0.000642

Table 2. Table 2: List of reactions in our reduced network, the rate coefficient symbol, and reference to the text. Symbols g ( Z < 0 ) 𝑍 0 (Z<0) and g ( Z > 0 ) 𝑍 0 (Z>0) indicate grains with negative and positive charge, respectively, while g ( Z + 1 ) 𝑍 1 (Z+1) and g ( Z − 1 ) 𝑍 1 (Z-1) indicate the reactant grain plus or minus one charge.

H2		$\to$	Mg+ + e-	$k_{H_2}$	$k_{H_2}=\zeta$
Mg+	+ e-	$\to$	H2	$k_{\mathrm{rec}}$	Eq. (26)
Mg+	+ g$\left(Z>0\right)$	$\to$	H2 + g$\left(Z+1\right)$	$k_{i,j}^{+}$	Eq. (33)
Mg+	+ g0	$\to$	H2 + g+	$k_{i,j}^0$	Eq. (34)
Mg+	+ g	$\to$	H2 + g$\left(Z+1\right)$	$k_{i,j}^{-}$	Eq. (31)
e-	+ g$\left(Z>0\right)$	$\to$	g$\left(Z-1\right)$	$k_{i,j}^{-}$	Eq. (31)
e-	+ g0	$\to$	g-	$k_{i,j}^0$	Eq. (34)
e-	+ g	$\to$	g$\left(Z-1\right)$	$k_{i,j}^{+}$	Eq. (33)
g-	+ g+	$\to$	g0 + g0	$k_{i,j}^{-}$	Eq. (31)
g–	+ g++	$\to$	g0 + g0	$k_{i,j}^{-}$	Eq. (31)
g++	+ g0	$\to$	g+ + g+	$k_{i,j}^0$	Eq. (34)
g–	+ g0	$\to$	g- + g-	$k_{i,j}^0$	Eq. (34)
g++	+ g-	$\to$	g0 + g+	$k_{i,j}^{-}$	Eq. (31)
g–	+ g+	$\to$	g0 + g-	$k_{i,j}^{-}$	Eq. (31)

Table 3. Table 3: Coefficient for the electron sticking S ​ ( T , Z ) 𝑆 𝑇 𝑍 S(T,Z) . See text for further details and its implementation in the rates.f90 file.

$Z$	$c_{Z,1}$	$c_{Z,2}$	$c_{Z,3}$
$-4$	-0.41953296	0.37771378	0.06950703
$-3$	-0.41819111	0.34261771	0.18453260
$-2$	-0.40908288	0.30922328	0.20339583
$-1$	-0.39206456	0.26730250	0.16134104
$0$	-0.36848365	0.22498770	0.05100280
$1$	-0.34350388	0.21428663	-0.17200160
$2$	-0.33385135	0.31951247	-0.60550265
$3$	-0.35192170	0.55369825	-1.17743478
$4$	-0.28404485	0.18312184	-0.70825975

Table 4. Table 4: Description and parameters of the different models along with the corresponding default/comparison. See the text for additional details.

Model name	Description	Reference parameter or model
reference	see Sect. 9.1	—
reference_noth	without radiative cooling or CR heating	reference
ideal	ideal MHD	reference
ideal_noth	ideal without radiative cooling or CR heating	ideal
adtab	${\eta }_{\mathrm{AD}}$ from equilibrium tables	time-dependent chemistry
arec	alternative recombination rate (Eq. 50)	recombination rate from Eq. (26)
stick01	electron sticking $S\left(T\right)=0.1$	fit to Bai (2011); see Sect. 4.3
stick1	electron sticking $S\left(T\right)=1$	fit to Bai (2011); see Sect. 4.3
N5e20	initial effective column density $N_{\mathrm{eff},1}=5\times {10}^{20}$ cm-2	$N_{\mathrm{eff},1}=5\times {10}^{21}$ cm-2
N5e22	initial effective column density $N_{\mathrm{eff},1}=5\times {10}^{22}$ cm-2	$N_{\mathrm{eff},1}=5\times {10}^{21}$ cm-2
cr$x$	constant cosmic ray ionisation rate $\zeta ={10}^{-x}$	cosmic ray attenuation as in Sect. 5
crlow	“low” CRs; $b_1=1.327\times {10}^{-12}$ s-1 and $b_2=-0.211$	“high” CRs; $b_1=5.34\times {10}^{-5}$ s-1 and $b_2=-0.384$
bulk10	dust bulk density ${\rho }_0=10$ g cm-3	${\rho }_0=3$ g cm-3
d2g1e4	dust to gas mass ratio $𝒟={10}^{-4}$	$𝒟={10}^{-2}$
d2g_step	$𝒟={10}^{-6}$ initially in the upstream region, $𝒟={10}^{-2}$ otherwise	$𝒟={10}^{-2}$
nogg	no grain-grain chemistry	grain-grain chemistry included
pexp25	$p=-2.5$ in $\phi \left(a\right)\propto a^p$	$p=-3.5$
pexp50	$p=-5$ in $\phi \left(a\right)\propto a^p$	$p=-3.5$
amin1e6	$a_{\min }={10}^{-6}$ cm in $\phi \left(a\right)\propto a^p$	$a_{\min }={10}^{-7}$ cm
cool01	cooling rate multiplied by $0.1$	standard cooling
cool10	cooling rate multiplied by $10$	standard cooling
noheat	no cosmic ray heating	cosmic ray heating included
nochem1e$x$	chemistry not solved; constant ionisation fraction $f_i={10}^{-x}$	time-dependent chemistry, ionisation fraction

Table 5. Table 5: Initial conditions for the MHD shock tube model.

Variable	Left state	Right state	Units
$n$	${10}^4$	$4\times {10}^4$	cm-3
$T$	${10}^3$	$10$	K
$B_x$	${10}^{-4}$	${10}^{-4}$	G
$B_y$	${10}^{-4}$	$2\times {10}^{-4}$	G
$B_z$	$0$	$0$	G
$v_x$	${10}^6$	$10$	cm s-1
$v_y$	$0$	$10$	cm s-1
$v_z$	$0$	$0$	cm s-1
$𝒟$	${10}^{-2}$	${10}^{-2}$	—
$f_i$	${10}^{-7}$	${10}^{-7}$	—

Table 6. Table 6: Coefficients for fitting function Eq. ( 56 ) applied to reaction rates that involve dust grains (Sect.  4.2 ). In the first column, for simplicity, we indicate only the reactants. Coefficients in the table use the notation a ​ ( b ) = a × 10 b 𝑎 𝑏 𝑎 superscript 10 𝑏 a(b)=a\times 10^{b} and have units of cm 3  s -1 .

Reactants	$c_0$	$c_1$	$c_2$	$c_3$	$c_4$	$c_5$
Mg+ + g	$-7.418\left(0\right)$	$2.331\left(-4\right)$	$1.619\left(-2\right)$	$-9.891\left(-3\right)$	$4.072\left(-3\right)$	$-3.404\left(-4\right)$
Mg+ + g-	$-5.268\left(0\right)$	$-5.931\left(-1\right)$	$1.688\left(-1\right)$	$-1.089\left(-1\right)$	$2.762\left(-2\right)$	$-2.041\left(-3\right)$
g+ + g--	$-6.110\left(0\right)$	$-5.995\left(-1\right)$	$1.778\left(-1\right)$	$-1.133\left(-1\right)$	$2.844\left(-2\right)$	$-2.091\left(-3\right)$
e- + g-	$-2.997\left(1\right)$	$3.028\left(1\right)$	$-1.622\left(1\right)$	$4.445\left(0\right)$	$-5.975\left(-1\right)$	$3.122\left(-2\right)$
e- + g+	$-2.946\left(0\right)$	$-5.931\left(-1\right)$	$1.688\left(-1\right)$	$-1.089\left(-1\right)$	$2.762\left(-2\right)$	$-2.041\left(-3\right)$
g + g++	$-8.310\left(0\right)$	$4.465\left(-3\right)$	$1.149\left(-2\right)$	$-7.213\left(-3\right)$	$3.670\left(-3\right)$	$-3.233\left(-4\right)$
e- + g++	$-2.705\left(0\right)$	$-6.448\left(-1\right)$	$2.285\left(-1\right)$	$-1.302\left(-1\right)$	$2.942\left(-2\right)$	$-2.013\left(-3\right)$
g+ + g-	$-6.266\left(0\right)$	$-5.478\left(-1\right)$	$1.109\left(-1\right)$	$-8.398\left(-2\right)$	$2.403\left(-2\right)$	$-1.883\left(-3\right)$
g++ + g--	$-6.266\left(0\right)$	$-5.478\left(-1\right)$	$1.109\left(-1\right)$	$-8.398\left(-2\right)$	$2.403\left(-2\right)$	$-1.883\left(-3\right)$
g++ + g-	$-6.146\left(0\right)$	$-5.890\left(-1\right)$	$1.646\left(-1\right)$	$-1.079\left(-1\right)$	$2.770\left(-2\right)$	$-2.062\left(-3\right)$
g + g--	$-8.310\left(0\right)$	$4.465\left(-3\right)$	$1.149\left(-2\right)$	$-7.213\left(-3\right)$	$3.670\left(-3\right)$	$-3.233\left(-4\right)$
Mg+ + g+	$-2.985\left(1\right)$	$2.590\left(1\right)$	$-1.336\left(1\right)$	$3.574\left(0\right)$	$-4.728\left(-1\right)$	$2.443\left(-2\right)$
Mg+ + g--	$-5.027\left(0\right)$	$-6.449\left(-1\right)$	$2.285\left(-1\right)$	$-1.302\left(-1\right)$	$2.941\left(-2\right)$	$-2.013\left(-3\right)$
e- + g	$-5.096\left(0\right)$	$2.384\left(-4\right)$	$1.618\left(-2\right)$	$-9.888\left(-3\right)$	$4.071\left(-3\right)$	$-3.403\left(-4\right)$

Table 7. Table 7: List of reactions in UN90 considering the assumptions made (see text). Since the evolution of neutral species is not tracked, they are omitted as products, leading to a lack of mass conservation. Note, however, that the charge is instead correctly conserved. Reactions 1, 6, and 13 assume that O + is instantaneously converted into OH + , i.e., m + , while reaction 23 assumes that H + 3 superscript subscript absent 3 {}_{3}^{+} is immediately formed from H + 2 superscript subscript absent 2 {}_{2}^{+} , and any H produced in reaction 24 is omitted.

1	H+	+	O	$\to$	m+
2	H+	+	O2	$\to$	m+
3	H+	+	M	$\to$	M+
4	He+	+	H2	$\to$	H+
5	He+	+	CO	$\to$	C+
6	He+	+	O2	$\to$	m+
7	H3+	+	CO	$\to$	m+
8	H3+	+	O	$\to$	m+
9	H3+	+	O2	$\to$	m+
10	H3+	+	M	$\to$	M+
11	C+	+	H2	$\to$	m+
12	C+	+	O2	$\to$	m+
13	C+	+	O2	$\to$	m+
14	C+	+	M	$\to$	M+
15	m+	+	M	$\to$	M+
16	H+	+	e-	$\to$	no products
17	He+	+	e-	$\to$	no products
18	H3+	+	e-	$\to$	no products
19	H3+	+	e-	$\to$	no products
20	C+	+	e-	$\to$	no products
21	m+	+	e-	$\to$	no products
22	M+	+	e-	$\to$	no products
23	H2			$\to$	H3+	+	e-
24	He			$\to$	He+	+	e-
25	H2			$\to$	H+	+	e-

Equations

$$\frac{\partial \mathbf{U}}{\partial t}=-\frac{\partial \mathbf{F}(\mathbf{U})}{\partial x}+\mathbf{S}(\mathbf{U}),$$

$${\partial }_t\rho$$

$${\partial }_t\left[\rho v_x\right]$$

$${\partial }_t\left[\rho v_y\right]$$

$${\partial }_t\left[\rho v_z\right]$$

$${\partial }_t{B}_x$$

$${\partial }_t{B}_y$$

$${\partial }_t{B}_z$$

$$+\frac{{\eta }_{\mathrm{AD}}}{{\mathbf{B}}^2}\left(F_{B,z}{B}_x-F_{B,x}{B}_z\right)],$$

$${\partial }_{t}E$$

$${\partial }_t\zeta$$

$${\partial }_t\left[\rho X_i\right]$$

$$P^{\ast }=P+\frac{{\mathbf{B}}^2}{8\pi },$$

$$P=(\gamma -1)\left(E-\frac{\rho {\mathbf{v}}^2}{2}-\frac{{\mathbf{B}}^2}{8\pi }\right),$$

$$P=\frac{\rho k_{\mathrm{B}}}{\mu m_{\mathrm{p}}}T,$$

$$F_{B,x}$$

$$F_{B,y}$$

$$F_{B,z}$$

$${\eta }_{\mathrm{AD}}=c^2\left(\frac{{\sigma }_{\mathrm{P}}}{{\sigma }_{\mathrm{P}}^2+{\sigma }_{\mathrm{H}}^2}-\frac{1}{{\sigma }_{\parallel }}\right),$$

$${\mathbf{F}}_i=\frac{{\mathbf{F}}_{i+1/2}-{\mathbf{F}}_{i-1/2}}{\Delta x},$$

$${\mathbf{F}}_{i+1/2}=\frac{{\alpha }^{+}{\mathbf{F}}_i+{\alpha }^{-}{\mathbf{F}}_{i+1}-{\alpha }^{+}{\alpha }^{-}\left({\mathbf{U}}_{i+1}-{\mathbf{U}}_i\right)}{{\alpha }^{+}+{\alpha }^{-}},$$

$${\alpha }^{\pm }=\max \left|\pm {\lambda }_i^{\pm },\pm {\lambda }_{i+1}^{\pm },0\right|,$$

$$c_{\mathrm{f}}^2=\frac{1}{2}\left(\theta +\sqrt{{\theta }^2-4c_{\mathrm{s}}^2\frac{B_x^2}{4\pi \rho }}\right),$$

$${\Lambda }_{\mathrm{chem}}=4.2\times 10^{-31}n({\mathrm{H}}_2)T^{3.3}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3}{\mathrm{s}}^{-1},$$

$${\Gamma }_{\mathrm{CR}}=\zeta Q(n_{\mathrm{H}})n({\mathrm{H}}_2)\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3}{\mathrm{s}}^{-1},$$

$$Q(n_{\mathrm{H}})=\sum _{i=0}^5{c}_i\log {\left(n_{\mathrm{H}}\right)}^i\mathrm{eV},$$

$$k_{\mathrm{rec}}(T)=k_0{\left[\sqrt{\frac{T}{T_0}}{\left(1+\sqrt{\frac{T}{T_0}}\right)}^{1-b}{\left(1+\sqrt{\frac{T}{T_1}}\right)}^{1+b}\right]}^{-1},$$

$${\mathcal{P}}_i$$

$$X_i\rho {\mathcal{L}}_i$$

$$⟨k(T)⟩=\frac{{\int }_{a_{\min }}^{a_{\max }}\phi (a)k(a,T)\mathrm{d}a}{{\int }_{a_{\min }}^{a_{\max }}\phi (a)\mathrm{d}a}$$

$$⟨k(T)⟩=\frac{{\int }_{a_{\min }}^{a_{\max }}{\int }_{a_{\min }}^{a_{\max }}\phi (a)k(a,a^{\prime },T)\mathrm{d}a\phi (a^{\prime })\mathrm{d}{a}^{\prime }}{{\int }_{a_{\min }}^{a_{\max }}{\int }_{a_{\min }}^{a_{\max }}\phi (a)\mathrm{d}a\phi (a^{\prime })\mathrm{d}{a}^{\prime }}$$

$$k_{i,j}^{-}(a_s,T)$$

$$v_g=\sqrt{\frac{8k_{\mathrm{B}}T}{{\mu }_{i,j}{m}_p}},$$

$$k_{i,j}^{+}(a_s,T)$$

$$k_{i,j}^0(a_s,T)=\pi v_{\mathrm{g}}{a}_s^2\left(1+\sqrt{\frac{\pi q^2{Z}_i^2}{2a_s{k}_{\mathrm{B}}T}}\right)S(T).$$

$$\Delta x_{\mathrm{eff}}=\frac{\Delta x}{\cos \vartheta \cos \varphi },$$

$$N_{\mathrm{eff},i}=N_{\mathrm{eff},1}+\sum _{j=2}^i\Delta x_{\mathrm{eff},j}{n}_j,$$

Full text

††pagerange: The challenges of modelling microphysics: ambipolar diffusion, chemistry, and cosmic rays in MHD shocks–H††pubyear: 2099

1 Introduction

Magnetic fields play a key role in determining the structure and evolution of many astrophysical environments. For example, in star forming regions, magnetic fields stabilise against gravity, influence the shape of molecular filaments, and play an important role in large-scale turbulence (see e.g. Padoan & Nordlund, 2002; McKee & Ostriker, 2007; Hennebelle & Chabrier, 2008; Federrath et al., 2010). In protoplanetary disks, magnetic fields can strongly influence the evolution of gas and dust (e.g. Balbus & Hawley, 1998; Armitage, 2011; Flock et al., 2016; Xu & Bai, 2016), remove angular momentum via outflows (e.g. Pudritz & Norman, 1983), affect the dynamical behaviour of planetesimals (e.g. Gressel et al., 2011), and even contribute to heating of (exo)planet atmospheres or to reduce atmospheric loss (e.g. Batygin et al., 2013; Cohen et al., 2014; Rogers & Showman, 2014; Dong et al., 2018).

How magnetic fields couple to gas and dust in these different contexts depends on the level of ionisation at each location and time. In weakly ionised conditions, such as can be found in molecular clouds and protoplanetary disks, one cannot generally assume that ideal magnetohydrodynamics (MHD) applies. If the magnetic diffusion timescale is comparable to the dynamical timescale, the coupling between magnetic fields and dynamics is regulated by microphysical processes that determine how the ionised matter is “felt” by the magnetic fields (e.g. Mestel & Spitzer, 1956; Wardle & Ng, 1999; Smith & Rosen, 2003; Duffin & Pudritz, 2008; Tomida et al., 2015). This can significantly affect the structure (e.g. outflow launching, disk formation) and dynamics (e.g. magnetic braking) of proto-stellar systems (see, e.g., Vaytet et al. 2018 and references therein).

Unfortunately, there remain significant uncertainties in certain aspects of astrophysically-relevant microphysics and chemistry (e.g. reaction rates, electron-grain sticking coefficients; Nishi et al. 1991; Bai 2011). These uncertainties could naturally affect models that include microphysics and lead to different outcomes when different ingredients are employed in, for example, the chemistry (e.g. Egan & Charnley, 1996; Ilgner & Nelson, 2006; Marchand et al., 2016; Wurster, 2016; Dzyurkevich et al., 2017), the dust content (e.g. Nishi et al., 1991; Okuzumi, 2009; Ivlev et al., 2016; Zhao et al., 2016), or cosmic rays (Padovani et al., 2013).

Our goal in this study is to determine which physical parameters are relevant for the evolution of the gas when self-consistently evolving time-dependent MHD alongside microphysics and chemistry (e.g. Kunz & Mouschovias, 2009; Xu & Bai, 2016), rather than by post-processing simulation snapshots (e.g. Padovani et al., 2013; Dzyurkevich et al., 2017), or employing pre-computed equilibrium tables (e.g. Gressel et al., 2011; Marchand et al., 2016). For this reason, we adopt a well-established, relatively simple framework for our experiments, i.e., a time-dependent, oblique, magnetic shock wave set in an environment that resembles the conditions of a pre-stellar core/dense molecular cloud (Lesaffre et al., 2004; Chen & Ostriker, 2012; Hollenbach et al., 2013; Flower & Pineau des Forêts, 2015; Holdship et al., 2017; Nesterenok, 2018). In this particular setting, the dominant non-ideal effect is ambipolar diffusion (Draine, 1980; Smith & Rosen, 2003; Duffin & Pudritz, 2008).

There have been several other studies of non-ideal MHD shocks in dusty plasmas using both steady-state (e.g. Pilipp & Hartquist, 1994; Wardle, 1998; Chapman & Wardle, 2006) and time-dependent approaches (e.g. van Loo et al., 2009; Ashmore et al., 2010; Van Loo et al., 2013) that examine different microphysical effects. Expanding upon these works, this paper aims to compare the importance of several microphysical ingredients and, in particular, what role cosmic rays play in determining the structure and evolution of MHD shocks.

We have developed, applied, and made publicly available a 1D, non-ideal MHD code and pre-processor111https://bitbucket.org/tgrassi/lemongrab/ to explore how common, simplifying assumptions about the microphysics affect the results relative to a full treatment of the problem. By varying the ingredients included in the experiments, we also identify simplifications to the chemistry and microphysics that do not strongly affect the results and are therefore worth exploring for possible use in large-scale, multi-dimensional, non-ideal MHD simulations.

Our self-consistent, albeit simplified, formulation of the problem also makes it possible to identify feedback processes in the chemistry/microphysics which are responsible for non-linear responses to variations of the external or internal parameters. For instance, as discussed in Sect. 9, the indirect effect of cosmic rays on the gas temperature via ambipolar diffusion heating.

The paper is structured as follows: First, in Sect. 2 we introduce the equations of non-ideal MHD and describe their implementation. In Sects. 3 to 6, we discuss the details and assumptions made for cooling and heating, chemistry, cosmic rays, and non-ideal microphysics. We verify the results produced by the code (described in Sect. 7) with a set of well-established tests in Sect. 8 before investigating how varying the microphysical ingredients affect the results in Sect. 9. We conclude in Sect. 10.

2 Methods: Non-ideal MHD 1D code

To test the effects of microphysics in a physically motivated, non-linear and time-evolving environment, we developed a 1D, time-implicit MHD code. The code evolves the physical quantities, $\mathbf{U}$, forward in time via

[TABLE]

where $\mathbf{F}$ are the fluxes, $\mathbf{S}$ the sources and sinks, and both are functions of $\mathbf{U}$. We define $\mathbf{U}=(\rho,\rho v_{x},\rho v_{y},\rho v_{z},B_{x},B_{y},B_{z},E,\rho X_{i})$, where $\rho$ is the mass density, $v_{i}$ and $B_{i}$ are the $i$th component of velocity and magnetic field, respectively, $E$ is the total energy density, and $X_{i}$ are the mass fractions of chemical species. $\mathbf{U}$ is defined at the centre of each cell.

Assuming that the ions and neutrals can be represented by a single fluid222cf. the more complex and numerically challenging multiple fluid approach (e.g. Ciolek & Roberge, 2002; Falle, 2003; O’Sullivan & Downes, 2006). (Shu et al., 1987; Choi et al., 2009; Masson et al., 2012), Eq. (1) can be explicitly written, including ambipolar diffusion terms, as333We employ Gaussian cgs units, i.e., the permeability $\mu_{0}=1$.:

[TABLE]

where $\mathbf{B}$ is the modulus of the magnetic field444$\mathbf{B}^{2}={B_{x}^{2}+B_{y}^{2}+B_{z}^{2}}$, and where $\partial_{t}\equiv\partial/\partial t$ and $\partial_{x}\equiv\partial/\partial x$ . Each chemical species is advected and their chemistry evolved according to Eq. (11); production ($\mathcal{P}_{i}$) and loss ($\mathcal{L}_{i}$) terms are discussed in Sect. 4.

The total pressure is

[TABLE]

while we assume an ideal equation of state for the thermal pressure

[TABLE]

and related to the temperature $T$, needed by the chemistry, via the ideal gas law

[TABLE]

where $k_{\rm B}$ is Boltzmann’s constant, $\gamma$ the adiabatic index, $\mu$ the mean molecular weight, and $m_{\rm p}$ the mass of the proton. The Lorentz force components are

[TABLE]

Spatial derivatives are evaluated using 2nd-order finite differences.

The ambipolar diffusion resistivity is given by

[TABLE]

where $c$ is the speed of light, $\sigma_{\parallel}$, $\sigma_{\rm P}$, and $\sigma_{\rm H}$ are the parallel, Pedersen and Hall conductivities, respectively, and will be discussed in Sect. 6.

The temporal evolution of the cosmic ray ionisation rate in each cell, $\zeta$, is calculated using Eq. (10), and discussed in detail in Sect. 5. Cosmic ray heating ($\Gamma_{\rm CR}$), as well as chemical cooling ($\Lambda_{\rm chem}$), are described in Sect. 3.

In Eqs. (2) to (11), since $\partial_{y}=\partial_{z}=0$, the solenoidal condition ($\nabla\cdot\mathbf{B}=0$) therefore requires that $\partial_{x}B_{x}=0$, which is guaranteed in the code by construction.

The complex, non-linear interplay between the myriad of different processes described by Eqs. (2) to (11) is summarised in Fig. 1: chemistry affects MHD via the resistivity coefficients ($\eta_{\rm AD}$), while MHD affects the energy ($E$), density ($\rho$), and magnetic field ($B$) evolution. Magnetic field and density determine the effective column density seen by cosmic rays ($N_{\rm eff}$), while chemistry depends on temperature via the reaction rate coefficients $k(T)$, on density, and on the cosmic ray ionisation rate ($\zeta$). Cosmic rays also affect temperature via direct heating ($\Gamma$).

2.1 HLL solver

We linearise the spatial derivatives on the right-hand side (RHS) of Eq. (1) using a standard HLL method (Harten et al., 1983) in order to numerically calculate the fluxes. The flux in the $i$th cell is given by

[TABLE]

where $i\pm 1/2$ denotes the quantity evaluated at the cell’s interface. These are defined as

[TABLE]

where

[TABLE]

$\lambda^{\pm}_{i}=v_{x}\pm c_{\rm f}$ are the eigenvalues of the Riemann problem at the cell interfaces, and the fast magnetosonic velocity evaluated at $i$ is

[TABLE]

with $\theta=v_{\rm A}^{2}+c_{\rm s}^{2}$, the speed of sound $c_{\rm s}^{2}=\gamma P/\rho$, and the Alfvén speed $v_{\rm A}^{2}=\mathbf{B}^{2}/(4\pi\rho)$.

2.2 Implicit time integration

We employ the DLSODES solver (Hindmarsh, 1983; Hindmarsh et al., 2005) to integrate the system (Eqs. (2) to (11)) forward in time. This approach avoids the need to explicitly define a time-step using a standard Courant condition; only absolute and relative tolerances of the individual quantities are needed (see below). The DLSODES solver has access to the RHS of Eq. (1) for all grid points and, for the sake of simplicity, we use the internally-generated Jacobian555For additional details, refer to the solver documentation contained in the opkdmain.f file.. Our approach is fully implicit in time and does not require any operator splitting to solve, e.g., the chemistry or cooling alongside the MHD. The code has been successfully validated against a set of standard numerical experiments which are discussed in Appendix A.

The accuracy of the method is determined by absolute ($\varepsilon_{\rm atol}$) and relative ($\varepsilon_{\rm rtol}$) tolerances666Tolerances are employed by the solver to compute the local error associated with the quantity $y$ as $\varepsilon_{\rm loc}=\varepsilon_{\rm rtol}|y|+\varepsilon_{\rm atol}$. Smaller tolerances increase the accuracy of the calculation, but could considerably increase the computational time. defined for each variable and each cell. We set777Units of $\varepsilon_{\rm atol}$ and $\varepsilon_{\rm rtol}$ are in code units, that assumes cgs. $\varepsilon_{\rm rtol}=10^{-8}$ for all variables, while $\varepsilon_{\rm atol}=10^{-30}$ for density ($\rho$), $\varepsilon_{\rm atol}=10^{-25}$ for the momentum ($\rho v_{i}$) and energy ($E$), $\varepsilon_{\rm atol}=10^{-10}$ for the magnetic field ($B_{i}$), and $\varepsilon_{\rm atol}=10^{-20}$ for the cosmic ray ionisation rate ($\zeta$). The chemistry (i.e. $\rho X_{i}$), meanwhile, uses $\varepsilon_{\rm atol}=10^{-30}$ for all species. In principle, the solver allows for different tolerances in different cells, but we find this unnecessary in the current study.

3 Methods: Cooling and Heating

In addition to the usual ambipolar MHD heating and cooling processes stated explicitly in Eq. (9), we also include a simplified radiative chemical cooling and direct heating from cosmic rays. The former is taken from Eq. (11) of Smith & Mac Low (1997):

[TABLE]

where $n({\rm H_{2}})$ is the molecular hydrogen number density in cm*-3*, and $T$ is the gas temperature in K. The cooling function employed here is accurate enough given the chemistry model we adopt for the current investigation (see Sect. 4). In Sect. 9, however, we test the effects of varying the strength of the cooling term.

We model the cosmic ray heating as

[TABLE]

where $\zeta$ is the cosmic ray ionisation rate (see Sect. 5), $n_{\rm H}$ is the number density of H nuclei ($n_{\rm H}=2n_{\rm H_{2}}$ in our chemical network), and $Q(n_{\rm H})$ is the heat deposited in the gas per ionisation event, taken from Fig. 2 of Galli & Padovani (2015) (see also Glassgold et al., 2012) and fit here using

[TABLE]

using the coefficients found in Tab. 1. The fitting function is valid in the range $n_{\rm H}=10^{2}$ to $10^{10}$ cm*-3*.

4 Methods: chemistry

To maintain a reasonable level of control over the many parameters in our model, we employ a simplified chemical network that follows the approach of Fujii et al. (2011) (see their Fig. 2 and our Tab. 2). Despite this reduced model, in Sect. 8, we demonstrate that this network is capable of reproducing the results of a few more complicated chemical networks. Our model assumes that the ionisation of H2 produces a cascade of fast reactions that lead immediately to e- and Mg+ as products, where the latter is a proxy for all cations. Analogously, Mg+ quickly recombines with electrons (and negatively charged grains) to reform H2. In our model, the molecular hydrogen ionisation rate coefficient $k_{\rm H_{2}}$ is equal to the cosmic-ray ionisation rate $\zeta$ (see Sect. 5). The rate of Mg+ recombination is obtained from Verner & Ferland (1996) in the form888http://www.pa.uky.edu/~verner/rec.html:

[TABLE]

with $k_{0}=1.92\times 10^{-11}$ cm3 s*-1*, $b=0.3028$, $T_{0}=4.849\times 10^{2}$ K, and $T_{1}=5.89\times 10^{6}$ K. Since Mg+ is a proxy for all positive ions, $k_{\rm rec}$ is therefore an effective recombination rate. In Sect. 9, we will examine the effects of varying $k_{\rm rec}$.

4.1 Differential equations for chemistry

Differential equations for the production rate $\mathcal{P}_{i}$ and loss rate $\rho X_{i}\mathcal{L}_{i}$ of the $i$th species are solved simultaneously with the equations of MHD in a single system, and are defined999Eq. (27) and Eq. (28) represent a standard set of differential equations for chemistry, but where the species abundances are given by their mass density instead of the more typical number density. by

[TABLE]

where $k_{r_{1},r_{2}}$ is the reaction rate coefficient between species $r_{1}$ and $r_{2}$, while $m_{i}$ and $X_{i}$ are the mass and the mass fraction of the $i$th species, such that $n_{i}m_{i}=\rho X_{i}$.

Together with the chemical network (defined in the previous Section), Eqs. (27) and (28) conserve the total number density, but not the total mass, because an H2 molecule is instantaneously converted into an Mg+ atom that is 24 times more massive. In principle, for the standalone chemical network, this does not represent a problem because the number density is conserved by construction. However, the hydrodynamics advects the mass density of the species, and it is therefore crucial to ensure conservation of mass. To avoid this issue, we define the mass of Mg+ as $m_{\rm Mg^{+}}^{\prime}=m_{\rm H_{2}}-m_{\rm e^{-}}$. When using the actual mass of Mg+, we find that the relative error on total mass conservation can be as large as $10^{-4}$ (instead of $\lesssim 10^{-7}$), while the error on global charge can reach $10^{-2}$ in the worst cases (instead of $\lesssim 10^{-8}$). We therefore use $m_{\rm Mg^{+}}^{\prime}$ for our models. We also note that a non-reduced chemical network will not, in general, be affected by this problem, since the mass will be correctly conserved by each reaction.

4.2 Grain chemistry

In order to determine the fraction of charged species to compute the resistivity coefficients (see Sect. 6), we include dust grains that can recombine with electrons and exchange charge with cations and amongst themselves (see Tab. 2). We integrate the grain size distribution $\varphi(a)$ over size range $a_{\rm min}$ to $a_{\rm max}$ to provide averaged reaction rate coefficients, $k(a,T)$, that are functions of the grain size $a$:

[TABLE]

for particle-grain interactions, and

[TABLE]

for grain-grain interactions, where $k(a,a^{\prime},T)$ are the rate coefficients for collisions of grains with radius $a$ and $a^{\prime}$, respectively.

Following Draine & Sutin (1987), for reactants with opposite charge ($Z_{i}Z_{j}<0$, e.g. electrons and positively charged grains), we include the Coulomb factor and the charge focusing due to polarisation as

[TABLE]

where $q$ is the elemental charge, $Z_{i}$ is the charge of the particle (e.g. electrons have $Z_{i}=-1$), $Z_{j}$ is the charge of the grain101010When two grains interact, $Z_{i}$ and $Z_{j}$ are the charge counts of the grains. Since Draine & Sutin (1987) consider only an interaction between a conducting sphere and a test charge (see their sect. II.a), we assume that $j$ is always the larger collision partner, i.e., a grain in the grain-particle collision, and the smaller grain in grain-grain interactions., $a_{s}=a_{i}+a_{j}$ is the sum of the grain sizes that reduces to $a_{s}=a_{j}$ when $i$ is a particle, $S(T)$ is the sticking coefficient (to be discussed below), and

[TABLE]

is the gas thermal velocity, $\mu_{i,j}=m_{i}m_{j}/(m_{i}+m_{j})$ is the reduced mass of the two species, and $m_{\rm p}$ the proton mass.

Analogously, for reactants with repulsive charges ($Z_{i}Z_{j}>0$, e.g., a cation and a positively charged grain), the rate is

[TABLE]

with $\theta_{v}=Z_{j}^{3/2}/[Z_{i}(\sqrt{Z_{i}}+\sqrt{Z_{j}})]$.

Finally, for $Z_{i}Z_{j}=0$ (e.g. a cation and a neutral grain), we have

[TABLE]

For grain-particle interactions, we compute $\langle k_{i,j}^{\pm,0}(T)\rangle$ by using Eq. (29) together with Eq. (31), Eq. (33), or Eq. (34), and assuming $a_{s}=a_{j}$. Grain-grain interactions are modelled analogously but using $a_{s}=a_{i}+a_{j}$ and Eq. (30) instead of Eq. (29). All of the rate coefficients discussed in this Section are shown for comparison in Fig. 2.

To speed-up code execution, we pre-compute the grain rate coefficients as a function of the temperature and apply a linear fitting function in logarithmic space at run-time. The grain size distribution properties are discussed in Appendix B, while fitting functions for reactions involving grains are given in Appendix C.

4.3 Sticking coefficient

To model the electron-grain sticking coefficient, $S(T)$, in Eq. (31) we refer to the appendix of Nishi et al. (1991) as well as Bai (2011). Electrons, because of their excess energy, only stick with some probability when they encounter a grain. A key parameter is $D$, the depth of the potential well between electrons and grains due to the polarisation interaction. We do not implement the equations reported there, but we fit Fig. 6 in Bai (2011) assuming $D=1$ eV and $a=0.1$ $\mu$m. The fit is $\log[S(T,Z)]=c_{Z,1}\log(T)^{2}+c_{Z,2}\log(T)+c_{Z,3}$, where $Z$ is the grain charge and the coefficients are listed in Tab. 3. Using $a=0.1$ $\mu$m (instead of integrating over the grain size distribution) may lead to errors, but since the choice of $D$ is arbitrary, we consider this fit accurate enough for the aims of this work. Moreover, comparing the two panels from Fig. 6 of Bai (2011), we note that the sticking coefficient with neutral grains (i.e. the main electron sticking route in our model) is almost size-independent. Nevertheless, in Sect. 9, we will vary the sticking coefficient and demonstrate its impact on the results.

As for cation-grain and grain-grain sticking coefficients, following Draine & Sutin (1987), we set $S(T)=1$ for the interaction between positive ions and any type of grain; in Eq. (31), Eq. (33), and Eq. (34), $S(T)<1$ only when the grain partner is an electron.

5 Methods: cosmic rays

Being charged particles, cosmic rays follow helical trajectories around magnetic field lines as they propagate. As a consequence, in the presence of a magnetic field, the effective column density $N_{\rm eff}$ “seen” by a cosmic ray can be much larger than the line of sight column density, especially if the magnetic field is not laminar (Padovani et al., 2013).

We consistently compute the propagation of cosmic rays following the approach of Padovani et al. (2018), where the cosmic ray ionisation rate of H2, $k_{\rm H_{2}}=\zeta$, is a function of the effective column density travelled by the particle, $\zeta\equiv f(N_{\rm eff})$, and is described in Appendix F of Padovani et al. (2018). Since we consider that our shock wave occurs in the vicinity of a pre-stellar core, prior to the cosmic rays entering our simulation domain, we assume they are partially attenuated by the surrounding medium. Thus, following Ivlev et al. (2015), we assume that the initial effective column density experienced by the cosmic rays is $N_{{\rm eff},1}=5\times 10^{21}$ cm*-2*. This value is evaluated at the centre of the first cell (thus, the “$1$” subscript). Since cosmic rays gyrate around magnetic field lines, and given the periodicity of the simulation domain along $y$- and $z$-directions, we compute the effective distance travelled as

[TABLE]

where $\vartheta=\arctan(B_{z}/B_{x})$ and $\phi=\arctan(B_{y}/\sqrt{B_{z}^{2}+B_{x}^{2}})$, while $B_{x}$, $B_{y}$, and $B_{z}$ are evaluated at the cell interface. We then compute the effective column density at the centre of the $i$th cell as

[TABLE]

where $\Delta x_{{\rm eff},j}$ and $n_{j}$ are calculated at the interface between cells $j-1$ and $j$ using a linear interpolation. Once the effective column density is available, we can retrieve the H2 ionisation rate $\zeta_{i}$ at the centre of the $i$th cell.

Within the column density range we are interested in here ($10^{20}<N<10^{23}$ cm*-2*), $\zeta_{i}$ is fairly represented by a power-law

[TABLE]

the coefficients $b_{1}$ and $b_{2}$ are discussed below.

Since $\cos{(\arctan{x})}=\left({x^{2}+1}\right)^{-1/2}$, using the definitions of $\vartheta$ and $\varphi$, we can rewrite Eq. (35) as

[TABLE]

and thus

[TABLE]

Substituting this into Eq. (37), and assuming that $\Delta x$ and $B_{x}$ are constant in time, we then differentiate with respect to time and obtain

[TABLE]

We note that the $\zeta_{i}$ in a given cell depends (non-trivially) on the densities and magnetic field values of all the cells from the first to the $i$th , i.e., on $3\times i$ variables. In practice, this considerably reduces the internal time step of DLSODES, since the number of variables that $\zeta_{i}$ depends on is large and the Jacobian becomes considerably less sparse. In fact, by using a constant $\zeta$, the integration time can be reduced by a factor of approximately one hundred.

In principle, when computing the propagation of cosmic rays, one should account for the effects of magnetic focusing and mirroring (Cesarsky & Volk, 1978; Desch et al., 2004; Padovani & Galli, 2011). Focusing and mirroring mechanisms act to amplify and reducing the cosmic ray flux, respectively, and could be important in regions of star formation. However, Silsbee et al. (2018) has demonstrated that these two effects nearly cancel each other out when the magnetic field strength has a single peak along the field lines, which is indeed the case in this work. Therefore, in the following, we choose to neglect mirroring and focusing effects.

The propagation of cosmic rays can also be affected by scattering due to self-generated Alfvén waves (Skilling & Strong, 1976; Hartquist et al., 1978), but this mechanism is only important at the edges and the more diffuse parts of a molecular cloud, and thus we can safely neglect it in this work.

5.1 Lower and upper bounds of $\zeta_{i}$

The cosmic ray ionisation rate at a given column density $N$ is given by

[TABLE]

where $j(E)$ is the cosmic ray differential flux (also called the spectrum), $\Phi$ is a multiplicity factor accounting for ionisation by secondary electrons, and $\sigma^{\rm ion}$ is the ionisation cross section. Since $\sigma^{\rm ion}$ is known to peak at low energies, the maximum contribution to $\zeta$ comes from cosmic rays in the energy range $10~{}{\rm MeV}\lesssim E\lesssim 1~{}{\rm GeV}$ (Padovani et al. 2009).

The most recent Voyager 1 data release (Cummings et al., 2016) leads to the conclusion that no upturn is expected in the interstellar proton spectrum down to energies of at least 3 MeV. The corresponding ionisation rate, however, is more than a factor of 10 smaller than estimates from observations in diffuse clouds (Indriolo et al., 2015; Neufeld & Wolfire, 2017). For this reason, as in Padovani et al. (2018), we consider two different models for the cosmic ray proton spectrum: a “low” spectrum, obtained by extrapolating the Voyager 1 data to low energies, and a “high” spectrum. The latter can be considered as an upper bound to the actual average galactic cosmic ray spectrum and provides an upper limit to the values of $\zeta$ estimated for diffuse clouds. The resulting ionisation rates and their comparison to observations is discussed in Ivlev et al. (2015).

The values of $b_{1}$ and $b_{2}$ in Eq. (37) that are needed to reproduce the two trends in $\zeta$ are $b_{1}=1.327\times 10^{-12}$ s*-1* and $b_{2}=-0.211$ for the “low” case, and $b_{1}=5.34\times 10^{-5}$ s*-1* and $b_{2}=-0.384$ for the “high” case. Note that the validity of the fit is limited to $10^{20}<N<10^{23}$ cm*-2*.

6 Methods: Non-ideal MHD coefficients

Ambipolar diffusion is controlled by the $\eta_{\rm AD}$ resistivity coefficient, appearing in Eqs. (7) to (9) and defined in Eq. (18). The resistivity coefficient depends on conductivities (i.e. $\sigma_{\parallel}$, $\sigma_{\rm P}$, and $\sigma_{\rm H}$) that are functions of temperature, magnetic field, and species abundances. Given the reduced chemistry that we include in our model (Sect. 4), we assume that the only interactions that affect the conductivity are collisions between charged particles (electrons, cations, and grains) and molecular hydrogen. In principle, if we were to follow Pinto & Galli (2008), each charged species could exchange momentum with any other species in the gas, but, since the momentum transfer is dominated by the interaction between charged species and H2 (which in our model is the main neutral component of the gas), we do not explicitly include all interactions. For the sake of completeness, however, we report all the possible interactions from Pinto & Galli (2008) in Appendix E.

6.1 Conductivities

The three conductivities (parallel, Pedersen, and Hall) are given by (e.g. Pinto et al., 2008):

[TABLE]

where $q$ is the elementary charge, $m_{i}$ is the mass of a charged particle, $qZ_{i}$ is its charge, $\rho_{i}$ its mass density, and $\beta_{i,n}$ is the Hall parameter which takes into account the interaction between charged particles and neutral species (in our case H2). The sum is over electrons, cations, and charged dust grains.

The Hall parameter for collisions between the $i$th charged particle (gas or dust) and a neutral species is given by

[TABLE]

where $\rho_{n}$ is the neutral gas mass density, $m_{n}$ its mass, and $R_{i,n}(T)$ is the momentum exchange rate coefficient, which is described in the next Section.

6.2 Momentum transfer rate coefficients

6.2.1 Charged grains–H2

To model the interaction between charged dust grains and molecular hydrogen, we follow Sect. 6 of Pinto & Galli (2008). We employ their Eq. (25) when the condition in their Eq. (23) is satisfied (the hard sphere approximation rate; $R_{\rm hs}$), otherwise we use their Eq. (A.3) (the Langevin rate; $R_{\rm L}$). Assuming that, at a critical grain size $a_{\rm c}$, Eq. (25) and Eq. (A.3) are equal, we can write

[TABLE]

where $\delta=1.3$ is taken from Liu et al. (2003) and $\alpha_{\rm pol}=8.06\times 10^{-25}$ cm3 is the polarisability of molecular hydrogen (Pinto & Galli, 2008).

Adopting a grain size distribution $\varphi(a)\propto a^{p}$ over size range $a_{\rm min}$ to $a_{\rm max}$ (Sect. 4.2), and the Langevin rate $R_{L}$ is size-independent, the rate coefficient then becomes

[TABLE]

where $m_{\rm H_{2}}$ is the mass of molecular hydrogen and, since the mass of the grain is larger than the mass of H2, the reduced mass is $\mu\sim m_{\rm H2}$. A simplified expression for Eq. (47) evaluated for the parameters stated above is reported in Appendix D.

We note that Eq. (47) is valid only when $a_{\rm min}\leq a_{\rm c}(T)\leq a_{\rm max}$. Substituting $Z=1$, $\delta=1.3$, and $\alpha_{\rm pol}=8.06\times 10^{-25}$ cm3 into Eq. (46), we find that $a_{\rm c}(T=1\,{\rm K})=1.71\times 10^{-7}$ cm and that it decreases as $a_{\rm c}\propto T^{-2.5}$, hence the validity of Eq. (47) is only critical for $a_{\rm min}$ because $a_{\rm max}\gg a_{\rm c}$ when $T\geq 1$ K. Conversely, using the same relation, $T_{\rm c}=8.6$ K is the critical temperature corresponding to $a_{\rm min}=10^{-7}$ cm. Therefore, when $T\geq T_{\rm c}$, Eq. (47) reduces to the second term on the right-hand side only (i.e. the normalised integral over $R_{\rm hs}$). We report the total rate and individual terms in Fig. 3. Since $Z$ does not strongly affect $\langle R_{g,n}(T,Z)\rangle$, we do not discuss the behaviour of larger $Z$ here, although we do include $Z=\pm 2$ in the Figure for comparison.

6.2.2 Electrons–H2

At low energies ($\lesssim 1$ eV), the collisional rate between electrons and molecular hydrogen deviates significantly from the Langevin approximation. Therefore, we employ the fit from Pinto & Galli (2008) based on the cross-section obtained by comparing theoretical and experimental data:

[TABLE]

where all variables are in cgs units. More details can be found in Pinto & Galli (2008).

6.2.3 Cations–H2

Analogously, the rate for collisions between positive ions and molecular hydrogen is also taken from Pinto & Galli (2008, A.3)

[TABLE]

and is the same as discussed in Sect. 6.2.1 when $a<a_{\rm c}$.

7 Code structure

In this Section, we provide a brief overview of the publicly available111111https://bitbucket.org/tgrassi/lemongrab/ code, lemongrab, developed for this study.

Following the approach of Krome (Grassi et al., 2014), the code consists of a Python pre-processor that computes the chemical reaction rates including dust grains, plus the momentum exchange cross-sections for the resistivity coefficient, and then writes optimised Fortran code that contains the MHD solver and other physics modules.

In contrast to Krome, in lemongrab, the Fortran files are directly modified by Python via specific directives that are recognised by the pre-processor as writable code blocks. The first stage is controlled by main.py, which creates an instance of the chemical network class (network.py) from an external file containing the reaction rate coefficients, and parses them into a set of objects according to the reaction class (reaction.py) and the species class (species.py). The reaction class also integrates any reaction rate coefficients that depend on the grain size distribution. Common variables, such as the grain size range $a_{min}$ to $a_{max}$, the exponent of the power law $\varphi(a)\propto a^{p}$, and the bulk density $\rho_{0}$, are defined in common.py. Dynamically-generated rate coefficient functions are written to rates.f90. The momentum exchange coefficients are also pre-computed by the pre-processor and supplied to the Fortran code (in nonideal.f90) via linear fitting to a logarithmically-spaced table in temperature. Finally, the Python pre-processor automatically generates and places the right-hand side of the chemical differential equations in odechem.f90.

The core of the second stage is the MHD solver (ode.f90), which is called by test.f90. The ode.f90 file contains the call to DLSODES which solves the complete set of differential equations, i.e., Eqs. (2) to (11).

The ode.f90 file also supplies the chemical differential equations (odechem.f90), the chemical reaction fluxes (fluxes.f90), and the rate coefficient constants (rates.f90) to DLSODES. Moreover, ode.f90 also accesses nonideal.f90, where the non-ideal coefficients calculation routine is contained, cooling.f90 and heating.f90 for cooling and heating processes, and the cosmic-ray propagation functions in crays.f90. The initial conditions for the MHD shock are stored in input.dat, while variables that remain constant during the simulation (e.g. $B_{x}$) are stored in commons.f90 to help the compiler in optimise the calculation.

A single variable n(physical_variables, cells) is used to represent the main data structure. It includes the values of the physical variables, $\mathbf{U}$ (Eq. 1), for all cells. This approach allows DLSODES to solve Eqs. (2) to (11) for all cells simultaneously and without any operator-splitting.

8 Code benchmark

In this Section, we compare our resistivity calculations to Marchand et al. (2016, hereinafter M16), which explores the behaviour of the non-ideal MHD coefficients using a zero-dimensional barotropic collapse problem and equilibrium chemistry. This benchmark was chosen because the source code is publicly available121212https://bitbucket.org/pmarchan/chemistry/ and the setup applies to the physical regimes we are studying here, making it an ideal target for comparison.

We have, meanwhile, also successfully tested the MHD implementation in lemongrab against two additional standard benchmarks. The results can be found in Appendix A.

Similar to Umebayashi & Nakano (1990, hereinafter UN90), the background model of M16 is a zero-dimensional collapsing cloud with a temperature determined using the equation of state from Machida et al. (2006), but also defined in Eq. (9) of M16 and using the parameters listed in their Eq. (10). The chemical evolution of the collapse is modelled from $n_{\rm H}=10^{2}-10^{25}$ cm*-3*, i.e., up until the formation of the second core. They set the ionisation rate to a constant $\zeta=10^{-17}$ s*-1*. In contrast to UN90, the dust in their model is given a power-law size distribution $\varphi(a)\propto a^{-3.5}$. The distribution is normalised in order to obtain the same total surface area as the fiducial distribution of Kunz & Mouschovias (2009), as shown in their Eqs. (16) and (17). We use the results reported in Figs. 3 and 5 of M16 as our reference for comparison.

The test consists of running the chemistry forward in time until equilibrium is reached at each density and then calculating the ambipolar diffusion ($\eta_{\rm AD}$), Ohmic resistivity ($\eta_{\rm O}$), and Hall resistivity ($\eta_{\rm H}$) coefficients. We initialise the temperature and magnetic field as functions of number density following M16. In our case, we limit the density range to $1-10^{12}$ cm*-3* because we miss some physical processes that are important at higher densities/temperatures, e.g., grain sublimation (see Fig. 2 of M16). Grains are permitted have a charge from $Z=-2$ to $Z=2$.

In order to reproduce their results, we use a grain size distribution $\varphi\propto a^{-3.5}$ with size range $a_{\rm min}=1.81\times 10^{-6}$ to $a_{\rm max}=9.049\times 10^{-5}$ cm, and a dust-to-gas ratio of $\mathcal{D}=0.1$. We also adopt the recombination rate coefficient from M16

[TABLE]

instead of our Eq. (26). In this test we assume sticking coefficient $S=1$ for all the rates involving grains, except for the electron-neutral grain attachment where $S=0.6$. Moreover, we do not use the Mg+ reduced mass, but the actual one.

We report out results in Fig. 4, where we find good agreement with M16 in both chemical abundances and resistivity coefficients. See also Appendix F for further details.

9 The effects of microphysics on the structure of magnetic shocks

The aim of this study is to understand the effects of chemistry, dust microphysics, cosmic rays, and cooling/heating on the evolution of astrophysical magnetic shocks. In this Section, we analyse these effects in detail by varying the physical ingredients and parameters of a reference model. The different models and their characteristics are reported in Tab. 4.

9.1 Reference model

All of the tests presented in this Section are based on a 1-D MHD reference shock tube model with a box size of $L_{\rm box}=3\times 10^{17}$ cm ($\simeq 0.1$ pc $\simeq 2\times 10^{4}$ AU) and $1024$ linearly-spaced grid points. The shock moves from left to right (see Fig. 5 and Tab. 5), with initial left state density $n_{\rm L}=10^{4}$ cm*-3*, velocity components $v_{x,{\rm L}}=10^{6}$ cm s*-1* and $v_{y,{\rm L}}=0$ cm s*-1*, magnetic field $B_{y,{\rm L}}=10^{-4}$ G, and temperature $T_{\rm L}=10^{3}$ K. The unperturbed right side has $n_{\rm R}=4\times 10^{4}$ cm*-3*, velocity components $v_{x,{\rm R}}=v_{y,{\rm R}}=10$ cm s*-1* (note $\mathbf{v}_{\rm L}\gg\mathbf{v}_{\rm R}$), $B_{y,{\rm R}}=2\times 10^{-4}$ G, and $T_{\rm R}=10$ K. Both sides have $B_{x,{\rm L}}=B_{x,{\rm R}}=10^{-4}$ G constant in time, $v_{z,{\rm L}}=v_{z,{\rm R}}=0$, and $B_{z,{\rm L}}=B_{z,{\rm R}}=0$. The interface between left and right states is placed at $L=0.3\,L_{\rm box}$. In all tests, we let the system evolve for $t=10^{4}$ yr.

The reference model includes the full calculation of ambipolar diffusion, time-dependent chemistry with the recombination rate coefficient from Eq. (26), and electron-sticking following Bai (2011). We set the dust-to-gas mass ratio to $\mathcal{D}=\rho_{d}/\rho=10^{-2}$ and the grain size distribution to $\varphi(a)\propto a^{p}$ with $p=-3.5$ and size range $a_{\rm min}=10^{-7}$ cm to $a_{\rm max}=10^{-5}$ cm. The dust is given a bulk density of $\rho_{0}=3$ g cm*-3*. The initial ionisation fraction is set to $f_{i}=n_{\rm e^{-}}/n_{\rm H_{2}}=10^{-7}$, but this has no influence on the evolution except for tests without chemistry (i.e. nochem1e$x$). See Appendix G for additional details on the chemical initial conditions.

Since the gas is dominated by molecular hydrogen, we assume a molecular gas with constant $\gamma=7/5$ and constant mean molecular weight $\mu=2$. Cosmic rays propagate from left to right with initial effective column density $N_{{\rm eff},1}=5\times 10^{21}$ cm*-2*, and we assume high cosmic ray spectrum, i.e., Eq. (37) with $b_{1}=5.34\times 10^{-5}$ s*-1* and $b_{2}=-0.384$ (Sect. 5).

9.2 General behaviour of the models

We evolve the shock for $t=10^{4}$ yr and, in Figs. 6 and 7, we report131313An animation of the evolution of reference and ideal models is available at https://vimeo.com/286491689; the evolution without thermal processes (reference_noth, ideal_noth) is available at https://vimeo.com/290131160. the results for reference, ideal, reference_noth and ideal_noth models.

The solution to the ideal_noth test, which is the simplest physical model we explore, exhibits five features (from left to right): a fast shock, a slow rarefaction, a contact discontinuity, a slow shock, followed by another fast shock. While the leftmost and rightmost fast shocks are clear, the slow rarefaction, located near $\sim 1.8\times 10^{17}$ cm at 10 kyr, is very low amplitude and only just visible in $v_{y}$. The contact discontinuity and slow shock are, meanwhile, adjacent to each other with the transition between the two occurring at $\sim 2\times 10^{17}$ cm (at 10 kyr). The solution as a whole moves to the right at $\sim 3.5$ km s*-1*. The $z$-components of the magnetic field and velocity ($B_{z}$ and $v_{z}$) remain equal to zero since both are initially zero and there are no $z$-derivatives in the problem. Thus, the magnetic field vector does not rotate for these particular initial conditions, even though the code is capable of this.

The ideal model, which includes cooling (Sect. 3), demonstrates a considerably different solution. The shock structure is modified significantly by the cooling of the hot gas in the left initial state ($T_{\rm L}=10^{3}$ K $\rightarrow$ $\sim 25$ K) and the shock-heated gas in the intermediate region between left and right fast shocks (cf. the ideal_noth model). The wave speeds and the extent of the region between the leftmost fast shock and the contact discontinuity/slow shock are subsequently reduced. Even the self-similarity of the ideal_noth solution is broken.

The addition of ambipolar diffusion smears out gradients in the magnetic field ($B_{y}$ in this case) and subsequently heats the affected regions. This is particularly visible for the region downstream from the rightmost fast shock in the reference_noth model (Fig. 7). Note, however, that the diffusion of the magnetic field in the vicinity of the leftmost fast shock only becomes substantial once cooling is included (i.e. the reference model); this is a result of recombination, which is more efficient at lower temperatures (Fig. 2), and reduces the ionisation fraction. That said, even including ambipolar heating, magnetic diffusion does not drastically modify the structure and evolution of the solution; cooling has a much more significant effect.

Before turning to a more detailed examination of how the microphysics affects the shock solution, we first describe our primary means of presenting and comparing the different tests. Figure 8 reports the density in the different models at $t=10^{4}$ yr for the region $x=10^{17}$ – $2.4\times 10^{17}$ cm. The absolute value of the gas density is plotted in the four top panels, where the black solid line indicates the reference model, while the shaded grey area denotes the envelope values/extrema of the density across all models for visual reference. The lower panels show the relative density variation $r_{\rho}=(\rho_{\rm model}-\rho_{\rm reference})/\rho_{\rm reference}$ for each model, sorted from the largest $|r_{\rho}|$ to the smallest.

The upper bound of the grey envelope in Fig. 8 is mainly set by the ideal MHD test, which has the largest positive $r_{\rho}$. Conversely, the nochem1e7 model, with constant ionisation fraction $f_{i}=10^{-7}$, gives the largest negative $r_{\rho}$ values and is mostly responsible for the lower bound of the envelope. In the region of the slow shock ($x\sim 1.9\times 10^{17}$ cm), however, it is instead the ideal model which sets the lower bound of the envelope and the nochem1e7 model which sets the upper bound.

Analogously, Figs. 9, 10, and 11 report, respectively, on the modulus of the magnetic field ($\mathbf{B}$), the temperature ($T$), and the resistivity coefficient ($\eta_{AD}$) as well as the relative differences with respect to the reference model.

9.3 The effect of modified chemistry

Chemistry plays a key role in the evolution of the shock because the abundances of the ions that control the Hall parameter $\beta_{i,n}$, and hence the resistivity coefficient $\eta_{\rm AD}$ (see Sect. 6.1), are determined by the chemistry. Indeed, $\beta_{i,n}$ and $\eta_{\rm AD}$ provide the main interplay between chemistry and hydrodynamics (see Fig. 1). To understand which chemical processes are most influential to the shock evolution, we change the recombination efficiency (model arec), the electron sticking (models stick01 and stick1), or turn off the chemistry (nochem1e4, nochem1e5, nochem1e6, and nochem1e7). As already discussed, parameters that reduce the ionisation fraction should shift models towards stronger ambipolar diffusion.

9.3.1 Cation-electron recombination rate

In our model, cations only recombine with electrons via $k_{\rm rec}$ (Eq. 26). Positively-charged grains and electrons can “recombine”, but this is an aspect of the grain chemistry and sticking coefficient, which is discussed below. When we adopt an alternative $k_{\rm rec}$ (Eq. 50; the arec model), we note that the solution becomes “less ideal” relative to the reference model (e.g. Fig. 8). The recombination rate in the arec model is much more effective than the one used in the reference model. Thus, the arec model has a lower abundance of free electrons, which diminishes the global ionisation fraction and enhances the magnetic diffusion (Fig. 11).

Given our reduced network, wherein cations are represented by Mg+, Eq. (26) is the most appropriate recombination rate to use. Although the arec model demonstrates that one must be careful when choosing an effective recombination rate for reduced networks, Figs. 8 – 11 instead show that it is not the most important effect in determining the evolution and structure of the shock.

9.3.2 Electron-grain sticking coefficient

The electron-grain sticking coefficient, $S(T)$, dictates the likelihood that electrons attach to grains after a collision. This not only affects the abundance of free electrons, but also the fraction of negatively-charged grains. Since the grain-cation sticking coefficient is typically greater than the electron-cation rate (see Fig. 2), increasing the fraction of negatively charged grains then increases the probability that Mg+ recombines with grains.

Figure 12 shows the contributions of Mg+, e-, and g- to $\eta_{\rm AD}$ in the reference, stick01 ($S(T)=0.1$), and stick1 ($S(T)=1$) models. Note that, since the fitting function for $S(T)$ (Sect. 4.3) is of order 0.1 (see Fig. 6 in Bai 2011), the stick01 model is very similar to the reference model. An electron-grain sticking coefficient of $S(T)=1$, however, results in a decreased abundance of e- and subsequent increased abundance of g-. As can be seen in Fig. 13, the increase of g- sticking partners reduces the abundance of Mg+, but because electron-grain sticking is more prevalent than cation-grain sticking, there is still a net increase in g- grains. As Fig. 12 shows, negatively-charged grains are the most important contributor to $\eta_{\rm AD}$. Thus, increasing the sticking coefficient decreases the magnetic resistivity and moves the system towards the ideal case (Figs. 8 – 11).

It is clear that the sticking coefficient plays an important role in determining the strength of the non-ideal terms, but we point out that its modelling is subject to significant uncertainties, such as in the depth of the potential well $D$ (see Sect. 4.3). For example, in the cases reported in Bai (2011), and for the temperature range of the present study, the sticking coefficient has values in the range $0.1\lesssim S(T)\lesssim 0.7$ for $1\leq D\leq 3$ eV.

9.3.3 Constant ionisation fraction

The ionisation fraction is determined by the chemistry. Thus, if we turn off chemistry altogether and instead force a constant ionisation fraction, we naturally find quite different results relative to the reference model.

In models nochem1e4, nochem1e5, nochem1e6, and nochem1e7 the ionisation fraction is set to a constant $f_{i}=10^{-4}$ to $10^{-7}$, respectively. In these cases, the chemical initial conditions remain unaltered during the evolution, which means $n_{\rm g^{\pm}}=0$, $n_{\rm e^{-}}=n_{\rm Mg^{+}}=f_{i}n_{\rm H_{2}}$ (see Appendix G); grains remain neutral, and electrons and Mg+ never recombine.

The evolution and structure of the shock is therefore controlled by $f_{i}$ and, when the ionisation fraction is large (e.g. nochem1e4), the results approach the ideal MHD limit. Conversely, when it is low (e.g. nochem1e7), the magnetic field is strongly diffused and nearly passive (Fig. 9). In fact, the models with constant ionisation fractions of $10^{-4}$ and $10^{-7}$ produce the largest deviations from the reference model and set the bounds of the gray envelopes in Figs. Fig. 8 – Fig. 11.

Evidently, a consistently calculated ionisation fraction is critical for obtaining a physically realistic shock structure and evolution (e.g. Flower et al., 1985). While it is true that the nochem1e6 model results are relatively similar to the reference model, this is only known because we first did the calculation consistently. As such, we caution others from using constant ionisation fractions when calculating non-ideal coefficients, unless it has already been established with a full calculation that this is a good approximation.

9.4 The effect of the cosmic ray parameters

In the present set up, since we do not include any external radiation, cosmic rays are the main driver of ionisation, and their effect is therefore of paramount importance. This is similar to the conditions at high column depths in a quiescent molecular cloud with little ongoing star formation.

9.4.1 Initial effective column density

The first parameter we modify is the initial effective column density $N_{\rm eff,1}$, i.e, the assumed column density that the cosmic rays have travelled through before entering the simulation box at $x=0$ (see Fig. 5). The default value is $N_{\rm eff,1}=5\times 10^{21}$ cm*-2*, which is comparable to the accumulated column density from $x=0$ to $x\simeq 2\times 10^{17}$ cm (i.e. the position of the shock front at $t=10^{4}$ yr); given $n_{\rm tot}\simeq 10^{4}$ cm*-3* (the average density of the shocked gas), $N_{\rm shock}=n_{\rm tot}\cdot x\simeq 2\times 10^{21}$ cm*-2* (Eq. 36). If we instead choose $N_{\rm eff,1}=5\times 10^{20}$ cm*-2* (model N5e20), $N_{\rm shock}$ now dominates over $N_{\rm eff,1}$. Conversely, if we choose $N_{\rm eff,1}=5\times 10^{22}$ cm*-2* (model N5e22), the opposite is true. The relative importance of $N_{\rm eff,1}$ with respect to $N_{\rm shock}$ explains the behaviour of the N5e20 and N5e22 models in Figs. 8–11. In the first case, $N_{\rm shock}$ dominates, so the variation of column density along the shock is relevant, even more than it is for the default value ($N_{\rm eff,1}=5\times 10^{21}$ cm*-2*; reference model). In the second case, $N_{\rm eff,1}$ dominates and $N_{\rm shock}$ is nearly ignorable. This is also clear from Fig. 14, where the ionisation rate in the N5e22 model is almost constant throughout the simulation box. Indeed, in Figs. 8–11 and Fig. 14, the N5e22 model is most similar to the cr16 model, which employs a constant cosmic ray ionisation rate of $\zeta=10^{-16}$ s*-1* (also see below).

9.4.2 Constant ionisation and a “low” spectrum

Next, we evolved shock models with a set of constant ionisation rates, using $\zeta=10^{-18},10^{-17},10^{-16}$, and $10^{-15}$ s*-1*, and labelled with cr$x$ indicating a fixed $\zeta=10^{-x}$ s*-1*. In general, since a larger $\zeta$ generates a larger number of free electrons, and thus a higher ionisation rate, the closer to the ideal model the results will be. Analogously, when we use the “low” spectrum fit (model crlow), i.e., Eq. (37) with $b_{1}=1.327\times 10^{-12}$ s*-1* and $b_{2}=-0.211$, we observe greater non-ideal behaviour relative to the reference model. We also observe that the crlow model lies between cr16 and cr17 (Fig. 14).

Evidently, the treatment of cosmic ray attenuation and subsequent ionisation rate plays a role in the shock structure and evolution. In particular, because cosmic rays control the ionisation fraction in our set up, they play an important role in determining the resistivity coefficients. Meanwhile, from Fig. 14, it is clear that, when including cosmic ray attenuation consistently, it is difficult to obtain an ionisation rate comparable to the canonical, constant rate of $10^{-17}$ s*-1* (Spitzer & Tomasko, 1968), even if there is a substantial attenuating column between a region of interest and the source of cosmic rays. Even if one adopts a “low” energy CR spectrum (see Sect. 5.1), the resultant ionisation rate is still a factor of 2-3 above the canonical ionisation rate.

9.5 The effect of different dust parameters

9.5.1 Dust-to-gas mass ratio

As we already know from Sect. 9.3.2, charged grains play a key role in determining the non-ideal behaviour of the shock. Thus, if we decrease the dust-to-gas ratio from $\mathcal{D}=10^{-2}$ to $10^{-4}$ (model d2g1e4), we expect that the contribution to $\eta_{\rm AD}$ from (negative) grains will decrease. A reduction in the amount of dust will also naturally decrease the probability that Mg+ and e- stick to grains. These effects are shown in Fig. 12, where it can be seen that decreasing the dust-to-gas ratio produces a large decrease in $\eta_{\rm AD}$. In this case, $\eta_{\rm AD}$ is dominated by Mg+ while the contribution from g- grains is negligible. Consequently, the shock structure is similar to the ideal case.

Next, we consider a discontinuous dust-to-gas ratio with a very low value ($\mathcal{D}=10^{-6}$) in the initially upstream region (i.e. $x<0.3L_{\rm box}$), but a typical value ($\mathcal{D}=10^{-2}$) in the initially downstream region. This model is intended to mimic the propagation of a dust-free shock into a dense, cold, and dust-rich cloud. The upstream region of the d2g_step model ($x\lesssim 1.7\times 10^{17}$ cm; Fig. 12) is very similar to the same region in the d2g1e4 model. In contrast, the downstream region ($x\gtrsim 1.9\times 10^{17}$ cm) is similar to the reference model, where $\mathcal{D}=10^{-2}$. This behaviour suggests that the advection of dust (or lack thereof) can strongly affect $\eta_{\rm AD}$ and the shock structure. However, the partial overlap between d2g1e4 and reference models with respect to d2g_step indicates that the net effect is less trivial than the sum of two dust-to-gas ratios.

9.5.2 Grain size distribution

To understand the effect of the grain size distribution $\varphi\propto a^{p}$ on the shock evolution, we now modify141414We do not vary $a_{\rm max}$, since with a power-law distribution the smaller grains comprise the largest total surface. $a_{\rm min}$, $p$, and the bulk density $\rho_{0}$. The impact of these parameters on $\eta_{\rm AD}$ is reported in Fig. 12. First, in amin1e6, we adopt a larger minimum grain size ($a_{\rm min}=10^{-6}$ cm) relative to the reference model ($a_{\rm min}=10^{-7}$ cm). With the removal of smaller grains, similar to d2g1e4, $\eta_{\rm AD}$ decreases and is now dominated by Mg+.

Analogously, in pexp25, if we adopt a shallower power-law exponent ($p=-2.5$; the reference model has $p=-3.5$), the dust mass becomes more heavily distributed towards larger grains, with a net effect on $\eta_{\rm AD}$ that is similar, but slightly less evident, than in amin1e6. Evidently, small grains dominate $\eta_{\rm AD}$ and, in both cases, their removal results in a solution that tends toward the ideal case (albeit, not as strongly as in d2g1e4).

If we now adopt a steeper power-law exponent, $p=-5$ (pexp50), the resulting shift of dust mass to smaller sizes decreases the abundances of e- and Mg+ while enhancing the fraction of negatively charged grains (see, e.g., stick1; Sect. 9.3.2). This is visible in Fig. 13, where the relative abundance of Mg+ decreases relative to the reference, and becomes comparable to the fraction of g- grains. The effect of the steeper power-law is thus a smaller overall value of $\eta_{\rm AD}$ (similar to stick1), but still higher than the pexp25 and amin1e6 models.

Finally, by increasing the bulk grain density from 3 to 10 g cm*-3* (model bulk10), we find that $\eta_{\rm AD}$ is somewhat reduced (relative to reference and similar to pexp50), but dominated by Mg+ instead of g- grains (Fig. 12). Since we are keeping the total dust mass constant, increasing the bulk density effectively decreases the number of grains, which results in less recombinations of cations and electrons with charged grains, and therefore larger abundances of these species (Fig. 13).

9.5.3 Grain-grain reactions

In the nogg model, we examine the effect of grain-grain reactions by removing them from the network. As can be seen in Figs. 8–11, 12, the effects on the shock structure and evolution are negligible. In the current context, grain-grain reactions are only marginally involved; the chemistry is dominated by the production of Mg+ and e- from H2 ionisation via cosmic rays, followed by the interaction of electrons and cations with grains to form g- and g+, which then eventually recombine with Mg+ to form H2. Similarly, as was also shown in M16, grains with charge $Z=\pm 2$ do not play a key role in the chemistry.

9.6 The effect of changing the thermal processes

To better understand the temperature structure of the shock (Fig. 10), we compute the cooling time by combining Eqs. (13) and (14) and taking the time derivative (see Appendix H for details):

[TABLE]

We integrate this equation numerically151515Using the odeint solver from the scipy package. assuming typical values $n_{\rm tot}=10^{4}$ cm*-3*, $\zeta=2.5\times 10^{-16}$ s*-1*, and initial temperature $10^{3}$ K. The results are shown in Fig. 15.

We note that the cooling time $\tau_{\rm cool}$ follows $\Lambda\propto T^{3.3}$ (Eq. 23) and, given that we evolve the system for $t=10^{4}$ yr, it is not surprising that the gas decreases to $T<100$ K. The gas does not have time to cool further, however, and thus does not reach the equilibrium temperature where $\Lambda=\Gamma_{CR}$ (which corresponds to a vertical line in Fig. 15, i.e., $\tau_{\rm cool}\to\infty$). If we now reduce the cooling function by a factor of 10 ($0.1\times\Lambda$ in Fig. 15), the cooling time increases, and the temperature in the shock and upstream region increases (see also Fig. 10). Analogously, if we increase the cooling by a factor of 10 ($10\times\Lambda$ in Fig. 15), the gas naturally reaches lower temperatures ($T<10$ K) in a shorter time.

This trend is confirmed in Fig. 10, where, for the cool01 model, the temperature is higher relative to the reference everywhere but in the quiescent downstream region. Indeed, the cool01 model produces the highest temperatures, in general, of all the models examined. For cool10, meanwhile, the opposite is true and, because the cooling time scale is now shorter, the quiescent downstream region even cools somewhat with respect to the initial conditions.

If we instead examine the effect of the cosmic ray heating rate $\Gamma_{\rm CR}$ on the cooling time, we find that the equilibrium temperature correlates with the ionisation rate $\zeta$. While it does not affect the cooling time scale where $\Lambda$ dominates, as can be seen in Fig. 15, an increased cosmic ray ionisation rate of $\zeta=10^{-15}$ s*-1* does raise the minimum temperature slightly, while decreasing the value to $\zeta=10^{-18}$ s*-1* decreases the equilibrium temperature to $T<4$ K.

That said, when the cosmic ray heating term is turned off entirely (model noheat), the impact on the temperature profile (see Fig. 10) is negligible, in contrast to what one would expect from Fig. 15. The temperature in models cr18 and cr15 (with constant cosmic ray ionisation rates of $\zeta=10^{-18}$ s*-1* and $10^{-15}$ s*-1*, respectively) does, meanwhile, deviate from the reference. In fact, the effect of a constant cosmic ray ionisation rate on the temperature is not direct via cosmic ray heating but rather indirectly through the chemistry and MHD heating processes.

The temperature behaviour in the different models suggest that it is primarily the ambipolar diffusion heating that prevents the temperature from reaching the levels observed in the ideal MHD model, which indeed provide the lowest temperatures observed in any of the models considered here. To demonstrate this, in Fig. 16, we plot the the ambipolar diffusion heating rate, which is given by the spatial derivative of the third term of Eq. (9). Since $\partial_{x}B_{x}=0$ and $B_{z}=0$, the heating rate can be written as

[TABLE]

Comparing the ambipolar heating rate (Fig. 16) with the temperature (Fig. 10), it is clear that the trends observed in the temperature of the different models can be explained by heating due to ambipolar diffusion.

10 Conclusions

We have developed an open-source161616https://bitbucket.org/tgrassi/lemongrab/ 1D, time-implicit, MHD code that includes ambipolar diffusion, chemistry, dust, and consistent cosmic-ray propagation. The code has been employed to explore the evolution of an oblique magnetic shock in order to understand the effects of the different microphysical parameters on the results. We have shown that, even in a simple application, microphysics plays a crucial role, and that the uncertainties are manifold.

Our main findings can be summarised as follows:

• •

In the absence of external radiation, cosmic rays plays a key role by controlling the ionisation level of the gas (via the chemistry) and hence determining the amplitude of the non-ideal MHD effects. In this study, we find that the widely-used cosmic-ray ionisation rate of $\zeta=10^{-17}$ s*-1* results in up to $\sim 90\%$ relative error in the density and $\sim 60\%$ in the magnetic field strength after $10^{4}$ yr relative to a self-consistent treatment of propagation.

• •

Dust is responsible for a large fraction of the neutral-ion momentum exchange. Reducing the dust-to-gas mass ratio from $\mathcal{D}=10^{-2}$ to $\mathcal{D}=10^{-4}$ results in a $\sim 280\%$ relative change in the density and a $\sim 60\%$ relative change in the magnetic field.

• •

The gas density is strongly affected by the parameters of the grain size distribution. Increasing the lower size limit of the distribution to $10^{-6}$ cm from $10^{-7}$ cm produces a change of $\sim 200\%$ in the density.

• •

Chemistry is also paramount. When turned off, the ionisation fraction is arbitrary (and constant in time), which affects the resistivity. A low ionisation fraction of $f_{i}=10^{-7}$ results in very strong ambipolar diffusion while, in contrast, values of $f_{i}=10^{-5}$ and $10^{-4}$ produce results that are very similar to the ideal MHD case.

• •

Reducing the cosmic-ray ionisation rate increases the temperature in certain regions, not because the direct cosmic-ray heating decreases, but because the ambipolar diffusion heating increases in these regions as a result of the lower ionisation. Analogously, higher ionisation rates (i.e. more highly ionised gas) show lower temperatures. This effect could be limited to the present set up, but due to its potential consequences, is worth exploring in more complex 3D simulations of pre-stellar cores and other environments.

• •

We find almost no change in the results when direct cosmic-ray heating is turned off or when grain-grain reactions are removed.

We remark that, given the complexity of the processes discussed in this paper and the interactions between them, our findings are particularly valid within the current set-up, and should not be arbitrarily generalised. For example, models with higher dimensions, complex chemistry, more realistic cooling functions, or other additional physics could change the importance of the different physical processes. To understand the interplay between the physical processes, the chemical model employed in this study was intentionally kept simple; determining the effects of a network with thousands of reactions, where many rates coefficients are uncertain, is beyond the aims of this paper. Nevertheless, the tests presented here show that self-consistent microphysics cannot be ignored in the context of non-ideal MHD simulations, and the choice of processes and parameters (mainly chemistry, cosmic rays, and dust) significantly affects the evolution of the dynamics.

Acknowledgement

We acknowledge Y. Fujii, O. Gressel, C. McNally, and R. Xu for useful comments and discussions. MP acknowledges funding from the European Unions Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 664931. The research leading to these results has received funding from the Danish Council for Independent Research through a Sapere Aude Starting Grant to TH. The Centre for Star and Planet Formation is funded by the Danish National Research Foundation (DNRF97). This research was supported by the DFG cluster of excellence “Origin and Structure of the Universe” (http://www.universe-cluster.de/). This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Ref no. FOR 2634/1 ER685/11-1. We would also like to thank the referee T. Hartquist for bringing additional relevant references to our attention.

Appendix A Algorithm testing

To validate the code implementation, we present here results for two standard MHD benchmarks. First, a Brio-Wu shock tube (Brio & Wu, 1988) with ideal MHD and, second, a C-shock tube (Masson et al., 2012) with MHD and ambipolar diffusion.

A.1 Brio-Wu MHD shock

To test the (ideal) MHD, we performed a Brio-Wu shock test (Brio & Wu, 1988) with $\mathbf{U}=\mathbf{0}$, except for $(\rho,P,B_{y})=(1,1,\sqrt{4\pi})$ where $x\leq 0.5$ and $(\rho,P,B_{y})=(0.2,0.2,-\sqrt{4\pi})$ otherwise, $B_{x}=\sqrt{4\pi}$, $\gamma=5/3$, $x\in[0,1]$. We used 1024 grid points. Note that for Ramses we employed rational units, i.e. $B_{\rm ramses}=B/\sqrt{4\pi}$. In Fig. 17, at $t=0.12$, we present our results. We also compared our HLL solver with the Ramses HLL and HLLD solvers, respectively with no slope limiter (slope_type=0) and MinMod limiter (slope_type=1) (Teyssier, 2002; Fromang et al., 2006; Teyssier et al., 2006). As expected, the solver without slope limiter is slightly more diffusive. We note that our results are indistinguishable from those obtained with the Ramses solver without slope limiter.

A.2 Non-ideal MHD C-shock

To determine if our ambipolar diffusion implementation is functioning properly, we benchmarked our code against the non-ideal MHD, non-isothermal C-shock presented in Masson et al. (2012, Sect. 2.4.2). The initial conditions are $\mathbf{U}=\mathbf{0}$, except for $(\rho,v_{x},v_{y},P,B_{y})=(0.5,5,0,0.125,\sqrt{2})$ where $x\leq 0.5$ and $(\rho,v_{x},v_{y},P,B_{y})=(0.988,2.5303,1.1415,1.4075,3.4327)$ otherwise, $B_{x}=\sqrt{2}$, $\gamma=5/3$, $\gamma_{\rm AD}=75$, $\rho_{\rm ion}=1$, $x\in[0,1]$. We used 400 grid points. Note that Masson et al. (2012) use rational units, hence $B_{\rm Masson}=B/\sqrt{4\pi}$. We evolve the system until the shock reaches a stationary state. In Fig. 18, we compare our results at $t=1$ to the analytic solution of Masson et al. (2012). As can be seen, other than a very slight difference in the position of the shock, we accurately reproduce the analytic solution.

Appendix B Grain size distribution properties

In this study, the grain size distribution is given by $\varphi(a)\propto a^{p}$ between sizes $a_{\rm min}$ and $a_{\rm max}$. The average grain mass is thus

[TABLE]

and from whence we compute the corresponding dust number density needed by the chemistry:

[TABLE]

where $\rho$ and $\rho_{\rm d}$ are the mass densities of gas and dust, respectively, $\mu_{g}$ is the dust mean molecular weight, and $\mathcal{D}$ is the dust-to-gas mass ratio.

To better compare our results with models that employ a constant grain size, we also derive the average surface area

[TABLE]

With $p=-3.5$, and $a_{\rm min}=10^{-7}$ to $a_{\rm max}=10^{-5}$ cm, this corresponds to an average grain size of $\sqrt{\langle a^{2}\rangle}=2.1\times 10^{-7}$ cm.

Appendix C Polynomial fitting functions for grain reactions

Here, we report the fitting functions for the reaction rates that involve dust grains discussed in Sect. 4.2. The grain distribution used here is $\varphi(a)\propto a^{p}$ with $p=-3.5$ and $a_{\rm min}=10^{-7}$ to $a_{\rm max}=10^{-5}$ cm. The fitting functions are of the form

[TABLE]

with $T$ in K, $k(T)$ in cm3 s*-1*, and coefficients $c_{i}$ given in Tab. 6. Fits are valid in the range $T=3$ K to $T=10^{4}$ K.

Appendix D Evaluating the charged grain–H2 momentum transfer rate

Eq. (47) evaluated for $p=-3.5$, $a_{\rm min}=10^{-7}$ cm, $a_{\rm max}=10^{-5}$ cm, $\delta=1.3$, and $\alpha_{\rm pol}=8.06\times 10^{-25}$ cm3 is equal to

[TABLE]

where $a_{1}=-2.7913\times 10^{-10}$, $a_{2}=1.7065\times 10^{-9}$, and $a_{3}=1.6369\times 10^{-9}$. All coefficients are in cgs units and the resulting rate is in ${\rm cm^{3}\,s^{-1}}$.

Appendix E Momentum transfer collisions

In Fig. 19, we report, for reference, the different processes described in Pinto & Galli (2008, henceforth P08b). Colours indicate the type of interaction and the corresponding equation number in P08b. Magnetic fields only interact with neutrals indirectly via charged colliders. The processes which are included in this work are enclosed by a dashed box. The “fit/other” label denotes fits to experiments or theoretical calculations (see Tab. 1 of P08b); “Langevin” and “Lgv” refer to the Langevin model (see also our Sect. 6.2.1); “H.Sphere” and “H.Sph” refer to the hard sphere approximation and means the rate equations are multiplied by $(1-S)$ where $S$ is the sticking coefficient; “Coulomb” is the standard Coulomb momentum transfer rate and the equations employed in P08b for these interactions is indicated. We refer the reader to P08b and references therein for additional details.

Appendix F Additional details for the Marchand et al. (2016) benchmark

In this Appendix, we include additional details to aid the interested reader in reproducing the results of M16.

The chemical network employed by M16 is reported in their Tab. A1 and includes the same reactions as in UN90. For the sake of completeness, however, the network is presented in Tab. 7. Here, we also report on their assumptions:

• •

All molecular ions except H${}_{3}^{+}$ (i.e. O${}_{2}^{+}$, HCO+, OH+, O2H+, CH${}_{2}^{+}$, and CO+) are represented by m+.

• •

Metallic, atomic cations (excluding oxygen and carbon) are indicated with M+.

• •

M and m are the neutral counterparts to M+ and m+, respectively.

• •

When O+ is produced, it immediately turns into OH+ (i.e. m+) by reacting with H2.

• •

Analogously, when H${}_{2}^{+}$ is produced (due to cosmic ray ionisation; see their Tab. 2), it turns immediately into H${}_{3}^{+}$.

• •

CO and C are the same species, since all neutral carbon is assumed to be in the form of CO.

• •

Neutral products are ignored, since their reservoir is assumed to be constant with time.

The network thus includes the following species: e-, O, O2, M, M+, H2, C, He, He+, m, m+, and H${}_{3}^{+}$. Neutral (g0) and charged dust (g±) is also included. Given the assumptions listed above, we obtain Tab. 7, the chemical network solved in M16. Note that, as written, mass is not conserved, but charge is.

The only difference between M16 and UN90 is given by the additional rates (A.1)-(A.3) in M16, which are taken from Eqs. (22)-(24) in Pneuman & Mitchell (1965), although the coefficients in those expressions are different from the ones reported by M16. Based on the code released by the authors171717Commit 2b23528., these rates should in fact be

[TABLE]

and are relevant when $T_{\rm gas}>10^{4}$ K and $n_{\rm tot}\gtrsim 10^{15}$ cm*-3*.

Dust evaporates as explained in Sect. 2.4.2 of M16 (see their Fig. 2); this behaviour can be reproduced by scaling the dust-to-gas mass ratio by a factor

[TABLE]

where

[TABLE]

with $\varepsilon=10^{-3}$, $T_{\rm mid}=(T_{\rm max}+T_{\rm min})/2$, and the sum is over grain species carbon (with parameters $w_{i}=0.85,\,T_{\rm min}=750\,{\rm K},\,T_{\rm max}=1100\,{\rm K})$, (MgFe)SiO4 $(w_{i}=0.144,\,T_{\rm min}=1200\,{\rm K},\,T_{\rm max}=1300\,{\rm K})$, and Al2O3 $(w_{i}=0.006,\,T_{\rm min}=1600\,{\rm K},\,T_{\rm max}=1700\,{\rm K})$.

We tested the validity of our assumptions using Krome (Grassi et al., 2014). The initial conditions are given in Tab. A2 of M16, and the dust is initially neutral with $n_{\rm d}=1.73\times 10^{-10}\,n_{\rm tot}$. In Krome, HCO*(+)* and Mg*(+)* are used as proxies for m*(+)* and M*(+), and we impose $\mathrm{d}n_{x}/\mathrm{d}t=0$ for the neutrals. We modelled the system until all species reach equilibrium. Empirically, we used $t_{\rm end}=10^{6}$ yr for models with $n_{\rm tot}<10^{10}$ cm-3* and $t_{\rm end}=1$ yr otherwise. We also imposed

[TABLE]

when $n_{\rm tot}\geq 10^{20}$ cm*-3* (i.e. when Eq. (F) dominates) in order to avoid an over-production of ions given the fact that we impose $\mathrm{d}n_{\rm H_{2}}/\mathrm{d}t=\mathrm{d}n_{\rm K}/\mathrm{d}t=\mathrm{d}n_{\rm Na}/\mathrm{d}t=0$.

The model results are identical to the results reported in M16.

Appendix G Initial abundances

We derive the initial abundances for each species from the total mass density ($\rho$), the dust-to-gas mass ratio ($\mathcal{D}$), and the ionisation fraction ($f_{i}=n_{\rm e^{-}}/n_{\rm H_{2}}$) as follows:

[TABLE]

which guarantees $\rho=\rho_{\rm H_{2}}+\rho_{\rm e^{-}}+\rho_{\rm Mg^{+}}+\rho_{\rm g(Z=0)}$ and global charge conservation.

Appendix H Cooling time

To compute the cooling time of a static volume of gas, we first derive its temperature evolution given an initial temperature $T(t=0)$, assuming a density $n_{\rm tot}\simeq n_{\rm H_{2}}$ and a cosmic ray ionisation rate $\zeta$. By equating Eq. (13) and Eq. (14), we find

[TABLE]

where $E^{*}$ is the sum of kinetic and magnetic energy, and we assume both to be constant in time for the current purpose. The time derivative of the temperature is thus

[TABLE]

Since we are considering a static volume of gas, i.e., all quantities are spatially constant, from Eq. (9) we obtain $\partial_{t}E=-\Lambda_{\rm chem}+\Gamma_{\rm CR}$, followed by

[TABLE]