© ESO 2021

1. Introduction

Magnetic fields are present in almost every astrophysical system, playing a fundamental role in its dynamics. Nevertheless, the issue of their origin remains an open question. There is consensus at present that the relatively strong magnetic fields measured across a range of scenarios, such as the solar surface, the interstellar medium, or galactic disks, cannot be primordial in origin but, rather, are the result of the amplification of a much weaker seed field.

Typically, the amplification process is described within the framework of ideal magnetohydrodynamics (MHD). Magnetic fields become stronger by means of their interaction with the flow of the plasma. Rotational and convective motions within the plasma produce a transformation of kinetic energy into magnetic energy, a phenomenon that is referred to as a MHD-dynamo (Larmor 1919; Cowling 1933; Moffatt 1978). However, the induction equation of ideal MHD describes the temporal evolution of a pre-existing magnetic field. That is to say, it does not allow the creation of new magnetic field when starting from a completely field-free initial condition. Therefore, additional non-ideal processes must be taken into account to explain the origin of magnetic fields.

Biermann (1950) showed that spatial variations of the density and pressure of electrons produce small charge imbalances that may give rise to equally small currents and magnetic fields. The condition for those currents to appear is that the electric field resulting from the charge imbalance must not be curl free, which implies that the gradients of pressure and density must be misaligned. This mechanism was proposed in the first place as an explanation for the origin of the magnetic fields of stars (see also Mestel & Roxburgh 1962). Later, it was applied to other astrophysical contexts such as the interestellar medium (Lazarian 1992; Kulsrud et al. 1997), the intergalactic medium (Subramanian et al. 1994; Gnedin et al. 2000; Kulsrud & Zweibel 2008), the early universe (Widrow et al. 2012), or the Sun’s local dynamo (Khomenko et al. 2017). Descriptions of how a weak magnetic field grows due to the effect of small-scale dynamo can be found in the works of Batchelor (1950) and Schekochihin et al. (2002), for instance.

The Kelvin-Helmholtz instability (KHI) is a common phenomenon in fluid and plasma dynamics. It occurs when perturbations appear at the interface of media moving with different velocities. There is a transfer of energy from the fluid flow to the perturbations so that their amplitudes start to increase. Initially, they follow a phase of exponential growth with a rate that is proportional to the shear flow velocity and inversely proportional to their length scale, as shown by the linear analysis of Chandrasekhar (1961). Then, during the so-called non-linear stage, the interface becomes greatly distorted due to the formation of large-scale vortexes and as the energy of the flow cascades to smaller scales (so that smaller and smaller vortexes are generated) a turbulent regime is reached (Krasny 1988). The process of inverse cascade in which smaller vortexes merge into larger scale vortexes may be also present in the KHI (Matsumoto & Seki 2010).

This instability has been observed and studied in a wide variety of astrophysical environments. For instance, it can be perceived in Earth’s aurora (Hallinan & Davis 1970), stellar jets and outflows (Keppens et al. 1999), supernova remnants (Wang & Chevalier 2001), protoplanetary disks (Gómez & Ostriker 2005), molecular clouds (Shadmehri & Downes 2007), coronal mass ejections (Foullon et al. 2011), and solar spicules (Antolin et al. 2018). Its relevance for the research presented in this work resides in the fact that the motions that develop during its evolution allow for the operation of the Biermann battery (Modestov et al. 2014) and dynamo (Zhang et al. 2009) mechanisms.

Most of the works previously mentioned in the context of the Biermann battery have focussed on the study of fully ionised plasmas, whereas those that have taken into account the effect of partial ionisation have resorted to single-fluid models. The range of applicability of these models is restricted to scenarios in which there is a strong coupling between the charged and neutral components of the plasma. To the best of our knowledge, no previous work has approached the investigation of the Biermann mechanism in partially ionised plasmas by using multi-fluid models, which are expected to have a broader range of applicability. On the other hand, the multi-fluid framework has already been used for the study of the KHI in partially ionised plasmas, as shown in Watson et al. (2004), Soler et al. (2012), Martínez-Gómez et al. (2015) or Hillier (2019). However, these authors did not include the battery term in their models.

Therefore, the goal of this article is to try to fill the above-mentioned gap and improve the understanding of the magnetogenesis processes in partially ionised plasmas. For that purpose, we use a two-fluid model that treats the neutral component of the plasma as a separate fluid from the charged one (which includes ions and electrons) but allows them to interact with each other by means of collisions. Here, we focus on fluid models with isotropic and scalar pressures. For a review of the so-called collisionless fluid models with anisotropic pressures (temperatures), see Hunana et al. (2019a,b). Our model uses a generalised Ohm’s law (Khomenko et al. 2014) and takes into consideration several collisional terms, such as the momentum and heat transfer due to elastic collisions and charge-exchange collisions in the momentum and energy equations, and a collisional term in the induction equation. Due to the complexity of the dynamics of the system we are interested in, we find it necessary to base our research on numerical simulations. Hence, the equations of our model are solved using the code MANCHA-2F (Popescu Braileanu et al. 2019a), although we do present an analysis of the equations to aid in understanding the numerical results. In addition, we performed a basic comparison between our numerical results and the predictions from single-fluid models.

For the time being, we restrict our simulations to the 2D case and use plasma parameters that correspond to the solar chromosphere. However, we expect that our results are general enough to be applied to other partially ionised plasmas and easily extrapolated to other astrophysical scenarios.

The present paper is organised as follows. In Sect. 2, we show the equations of the two-fluid model and describe the numerical setup for the simulation of the KHI. In Sect. 3, we present the results of several series of simulations. We describe the evolution of the instability and the growth of the magnetic field. Then, we study the dependence on the ionisation degree of the plasma, on the different collisions terms included in the model and on the numerical dissipation. In Sect. 4, we analyse the equations of our model to explain the results obtained from the simulations. Finally, in Sect. 5, we present a summary of the performed research and discuss some lines that can be followed to improve it.

2. Numerical model

In this work, we use a two-fluid model to describe the properties and the temporal evolution of a partially ionised plasma. The plasma is assumed to be composed of three different kinds of particles: ions, electrons and neutrals, which are denoted by the subscripts ‘i’, ‘e’, and ‘n’, respectively. However, a strong coupling between ions and electrons is assumed and they are taken into account together as a charged (or ionised) fluid denoted by the subscript ‘c’. The other fluid is composed of neutral particles only. The derivation of this two-fluid model from a more general multi-fluid description can be found in Khomenko et al. (2014).

The equations that govern the evolution of the density, velocity, and energy of each fluid and of the magnetic field are shown in Sect. 2.1. These equations are solved numerically using the MANCHA-2F code, which is a two-fluid extension of the original single-fluid MANCHA3D code (Khomenko & Collados 2006; Felipe et al. 2010; González-Morales et al. 2018). The details of the numerical implementation of the two-fluid model can be found in Popescu Braileanu et al. (2019a). The specific setup of initial and boundary conditions used in the present work is explained in Sect. 2.2.

2.1. Equations

The evolution of the density of each one of the two fluids is described by its respective continuity equation:

$$\begin{array}{c}\hfill \frac{\partial {\rho }_{\mathrm{s}}}{\partial t}+\mathrm{\nabla }·\left({\rho }_{\mathrm{s}}{\mathit{V}}_{\mathrm{s}}\right)=0,\end{array}$$(1)

where s ∈ {c, n}, and ρ_c = ρ_i + ρ_e. No source terms appear on the right-hand side because we do not consider the processes of ionisation and recombination.

Neglecting the electrons’ inertia due to their small mass, m_e ≪ m_i, the velocity of the charged fluid is given by the velocity of the ions, that is V_c = V_i. Neglecting also the effect of gravity, the velocities V_c and V_n are computed from the following momentum equations:

$$\begin{array}{c}\hfill \frac{\partial \left({\rho }_{\mathrm{c}}{\mathit{V}}_{\mathrm{c}}\right)}{\partial t}+\mathrm{\nabla }·({\rho }_{\mathrm{c}}{\mathit{V}}_{\mathrm{c}}{\mathit{V}}_{\mathrm{c}}+P_{\mathrm{c}}\mathbb{I})=\mathit{J}\times \mathit{B}+{\mathit{R}}_{\mathrm{cn}}\end{array}$$(2)

and

$$\begin{array}{c}\hfill \frac{\partial \left({\rho }_{\mathrm{n}}{\mathit{V}}_{\mathrm{n}}\right)}{\partial t}+\mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathit{V}}_{\mathrm{n}}{\mathit{V}}_{\mathrm{n}}+P_{\mathrm{n}}\mathbb{I})={\mathit{R}}_{\mathrm{nc}}.\end{array}$$(3)

The variable P_c represents the sum of the pressures of ions and electrons, P_n is the pressure of neutrals, 𝕀 is the identity tensor, J is the current density, and B is the magnetic field. Neglecting the displacement currents in Ampére’s law, the variables J and B are related by:

$$\begin{array}{c}\hfill \mathit{J}=\frac{\mathrm{\nabla }\times \mathit{B}}{{\mathrm{\mu }}_0},\end{array}$$(4)

where μ₀ is the magnetic permeability.

The terms R_cn and R_nc represent the momentum transfer due to collisions between two fluids moving with different velocities and are defined as

$$\begin{array}{c}\hfill {\mathit{R}}_{\mathrm{cn}}={\alpha }_{\mathrm{eff}}{\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}({\mathit{V}}_{\mathrm{n}}-{\mathit{V}}_{\mathrm{c}})=-{\mathit{R}}_{\mathrm{nc}}.\end{array}$$(5)

The collisional parameter α_eff takes into account the elastic collisions between ions and neutrals, between electrons and neutrals, and the charge-exchange collisions between ions and neutrals. Thus,

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{eff}}=\alpha +{\alpha }_{\mathrm{cx}},\end{array}$$(6)

where α is the contribution of the elastic collisions and α_cx is the contribution of the charge-exchange collisions.

The elastic collisional parameter is given by

$$\begin{array}{c}\hfill \alpha =\frac{{\rho }_{\mathrm{e}}{\nu }_{\mathrm{en}}+{\rho }_{\mathrm{i}}{\nu }_{\mathrm{in}}}{{\rho }_{\mathrm{n}}{\rho }_{\mathrm{c}}},\end{array}$$(7)

where ν_st is the frequency of collisions of particles ‘s’ with particles ‘t’. Following the works of Braginskii (1965) and Draine (1986), these frequencies are computed as

$$\begin{array}{c}\hfill {\nu }_{\mathrm{st}}=\frac{n_{\mathrm{t}}{m}_{\mathrm{t}}}{m_{\mathrm{s}}+m_{\mathrm{t}}}\sqrt{\frac{8k_{\mathrm{B}}}{\pi }(\frac{T_{\mathrm{s}}}{m_{\mathrm{s}}}+\frac{T_{\mathrm{t}}}{m_{\mathrm{t}}})}{\sigma }_{\mathrm{st}},\end{array}$$(8)

where m_s is the mass of a particle “s”, T_s is the temperature, n_t is the number density of species “t”, n_t = ρ_t/m_t, k_B is Boltzmann constant and σ_st is the collisional cross-section. Then, the parameter α can be computed through the following expression

$$\begin{array}{c}\hfill \alpha =\frac{m_{\mathrm{in}}}{m_{\mathrm{n}}^2}\sqrt{\frac{8k_{\mathrm{B}}}{\pi }\frac{T_{\mathrm{cn}}}{m_{\mathrm{in}}}}{\sigma }_{\mathrm{in}}+\frac{m_{\mathrm{en}}}{m_{\mathrm{n}}^2}\sqrt{\frac{8k_{\mathrm{B}}}{\pi }\frac{T_{\mathit{cn}}}{m_{\mathrm{en}}}}{\sigma }_{\mathrm{en}},\end{array}$$(9)

in which T_cn = (T_c+T_n)/2 is the average temperature and m_st = m_sm_t/(m_s+m_t) is the reduced mass of particles ‘s’ and ‘t’. For the sake of simplicity, in the present study we assume that the plasma is composed of hydrogen only. Therefore, m_i ≈ m_n is the proton mass, m_in = 0.5m_i and m_en ≈ m_e. We choose the following values of the collisional cross sections (Leake et al. 2013): σ_in = 1.16 × 10⁻¹⁸ m² and σ_en = 10⁻¹⁹ m². Then, α can be expressed as

$$\begin{array}{c}\hfill \alpha =\frac{1}{m_{\mathrm{n}}^2}\sqrt{\frac{8k_{\mathrm{B}}}{\pi }{T}_{\mathrm{cn}}}[\sqrt{0.5m_{\mathrm{i}}}{\sigma }_{\mathrm{in}}+\sqrt{m_{\mathrm{e}}}{\sigma }_{\mathrm{en}}],\end{array}$$(10)

and since m_e ≪ m_i and σ_en < σ_in, the value of this parameter is mainly determined by the ion-neutral interaction, with a much smaller contribution coming from the electron-neutral collisions.

According to the definition of collision frequencies given by Eq. (8), the consevation of momentum is satisfied:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{s}}{\nu }_{\mathrm{st}}={\rho }_{\mathrm{t}}{\nu }_{\mathrm{ts}}.\end{array}$$(11)

In addition, the expression for the collisional parameter α given by Eq. (7) allows us to define effective collision frequencies between charges and neutrals as

$$\begin{array}{c}\hfill {\nu }_{\mathrm{cn}}=\alpha {\rho }_{\mathrm{n}}\text{and}{\nu }_{\mathrm{nc}}=\alpha {\rho }_{\mathrm{c}}.\end{array}$$(12)

The contribution of the charge-exchange collisions is computed as (see e.g. Pauls et al. 1995; Meier & Shumlak 2012):

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{cx}}=\frac{1}{m_{\mathrm{n}}}(V_0^{\mathrm{cx}}+V_{\mathrm{cn}}^{\mathrm{cx}}+V_{\mathrm{nc}}^{\mathrm{cx}}){\sigma }_{\mathrm{cx}},\end{array}$$(13)

where the velocities $V_0^{\mathrm{cx}}$ and $V_{\mathrm{st}}^{\mathrm{cx}}$ are defined as

$$\begin{array}{c}\hfill V_0^{\mathrm{cx}}=\sqrt{\frac{4}{\pi }{v}_{\mathrm{Tc}}^2+\frac{4}{\pi }{v}_{\mathrm{Tn}}^2+v_{\mathrm{D}}^2},\end{array}$$(14)

$$\begin{array}{c}\hfill V_{\mathrm{st}}^{\mathrm{cx}}=v_{\mathrm{Ts}}^2{[4(\frac{4}{\pi }{v}_{\mathrm{Tt}}^2+v_{\mathrm{D}}^2)+\frac{9\pi }{4}{v}_{\mathrm{T}\mathrm{s}}^2]}^{-0.5}\end{array}$$(15)

and the cross section σ_cx is given by

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{cx}}=1.12\times {10}^{-18}-7.15\times {10}^{-20}ln\left(V_0^{\mathrm{cx}}\right).\end{array}$$(16)

The thermal speed of each species ‘s’ is given by

$$\begin{array}{c}\hfill v_{\mathrm{Ts}}=\sqrt{\frac{2k_{\mathrm{B}}{T}_{\mathrm{s}}}{m_{\mathrm{s}}}},\end{array}$$(17)

and

$$\begin{array}{c}\hfill v_{\mathrm{D}}=|{\mathit{V}}_{\mathrm{c}}-{\mathit{V}}_{\mathrm{n}}|\end{array}$$(18)

is the velocity drift between charged and neutral particles.

The evolution of the magnetic field is given by the induction equation, namely

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=-\mathrm{\nabla }\times \mathit{E},\end{array}$$(19)

where the electric field E is given by a generalised Ohm’s law. Derivations of Ohm’s law in different multi-fluid descriptions can be found in, for instance, Zaqarashvili et al. (2011), Martínez-Gómez et al. (2016) and Ballester et al. (2018). Here, we follow the approach used by Bittencourt (1986) and Khomenko et al. (2014), with the difference that we neglect the effect of resistivity and Hall’s term. Therefore, the Ohm’s law used in the present work is given by

$$\begin{array}{c}\hfill \mathit{E}=-{\mathit{V}}_{\mathbf{c}}\times \mathit{B}-{\eta }_{\mathrm{H}}\mathrm{\nabla }{P}_{\mathrm{e}}+{\eta }_{\mathrm{D}}({\mathit{V}}_{\mathrm{n}}-{\mathit{V}}_{\mathrm{c}}),\end{array}$$(20)

where

$$\begin{array}{c}\hfill {\eta }_{\mathrm{H}}=\frac{1}{e{n}_{\mathrm{e}}},\end{array}$$(21)

$$\begin{array}{c}\hfill {\eta }_{\mathrm{D}}=\frac{m_{\mathrm{e}}({\nu }_{\mathrm{en}}-{\nu }_{\mathrm{in}})}{e},\end{array}$$(22)

and e is the elemental charge.

Combining Eqs. (19) and (20) we get the following expression for the induction equation:

$$\begin{array}{cc}\hfill \frac{\partial \mathit{B}}{\partial t}& =\mathrm{\nabla }\times ({\mathit{V}}_{\mathrm{c}}\times \mathit{B})+\mathrm{\nabla }\times \left[{\eta }_{\mathrm{H}}\mathrm{\nabla }{P}_{\mathrm{e}}\right]\hfill \\ \hfill & -\mathrm{\nabla }\times \left[{\eta }_{\mathrm{D}}({\mathit{V}}_{\mathrm{n}}-{\mathit{V}}_{\mathrm{c}})\right]={\mathit{A}}_{\mathrm{Ind}}+{\mathit{B}}_{\mathrm{Ind}}+{\mathit{C}}_{\mathrm{Ind}},\hfill \end{array}$$(23)

where A_ind, B_Ind, and C_Ind are the advective, Biermann battery and collisional terms, respectively. The collisional term is usually neglected because in the presence of a background magnetic field it is typically very small in comparison with the advection term. However, in this work we are interested in a situation in which there is no initial magnetic field. Thus, we keep this term to check whether in this case it may have a relevant influence on the evolution of the system.

It is interesting to consider the scenario with no initial magnetic field and no collisions between the charged and the neutral fluids. In this case, the only term that appears on the right-hand side of the induction equation is the Biermann battery term. Using the definition given by Eq. (21) for the coefficient η_H, the induction equation has the following form:

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=\mathrm{\nabla }\times \left[\frac{\mathrm{\nabla }{P}_{\mathrm{e}}}{e{n}_{\mathrm{e}}}\right]=-\frac{\mathrm{\nabla }{n}_{\mathrm{e}}\times \mathrm{\nabla }{P}_{\mathrm{e}}}{e{n}_{\mathrm{e}}^2}.\end{array}$$(24)

This formula tells us that even when there is not an initial magnetic field seed, a magnetic field is generated if the spatial variations in the density and pressure of electrons are misaligned. It must be noted here that this effect will not be present if the flow is barotropic, that is, if the pressure is a function of the density only, P_e = P_e(n_e) (Kulsrud et al. 1997).

The equations for the evolution of the internal energy of each fluid, e_c and e_n, respectively, are

$$\begin{array}{c}\hfill \frac{\partial e_{\mathrm{c}}}{\partial t}+\mathrm{\nabla }·\left(e_{\mathrm{c}}{\mathit{V}}_{\mathrm{c}}\right)+P_{\mathrm{c}}\mathrm{\nabla }·{\mathit{V}}_{\mathrm{c}}=\mathit{J}·{\mathit{E}}_{\mathrm{diff}}+Q_{\mathrm{cn}}\end{array}$$(25)

and

$$\begin{array}{c}\hfill \frac{\partial e_{\mathrm{n}}}{\partial t}+\mathrm{\nabla }·\left(e_{\mathrm{n}}{\mathit{V}}_{\mathrm{n}}\right)+P_{\mathrm{n}}\mathrm{\nabla }·{\mathit{V}}_{\mathrm{n}}=Q_{\mathrm{nc}},\end{array}$$(26)

where

$$\begin{array}{c}\hfill {\mathit{E}}_{\mathrm{diff}}=\mathit{E}+{\mathit{V}}_{\mathrm{c}}\times \mathit{B}\end{array}$$(27)

and Q_cn and Q_nc are the heat transfer terms associated with collisions.

The internal energies are related to the pressure values by

$$\begin{array}{c}\hfill e_{\mathrm{c}}=\frac{P_{\mathrm{c}}}{\gamma -1}\mathrm{and}{e}_{\mathrm{n}}=\frac{P_{\mathrm{n}}}{\gamma -1},\end{array}$$(28)

and, since m_i ≈ m_n, the energy transfer terms due to collisions (see e.g. Schunk 1977; Draine 1986) are given by

$$\begin{array}{c}\hfill Q_{\mathrm{cn}}={\alpha }_{\mathrm{eff}}{\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}[\frac{k_{\mathrm{B}}}{(\gamma -1)m_{\mathrm{n}}}(T_{\mathrm{n}}-T_{\mathrm{c}})+\frac{v_{\mathrm{D}}^2}{2}]\end{array}$$(29)

and

$$\begin{array}{c}\hfill Q_{\mathrm{nc}}={\alpha }_{\mathrm{eff}}{\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}[\frac{k_{\mathrm{B}}}{(\gamma -1)m_{\mathrm{n}}}(T_{\mathrm{c}}-T_{\mathrm{n}})+\frac{v_{\mathrm{D}}^2}{2}].\end{array}$$(30)

We note here that the expressions used for the momentum and energy transfer due to collisions, that is, Eqs. (5), (29), and (30) are only valid when the velocity drifts between the two species are much smaller than the thermal speeds of the fluids. The general expressions include correction factors that depend on the ratios between those speeds (Schunk 1977; Draine 1986).

Finally, as equations of state we use the ideal gas law, which relates the pressure, number density, and temperature of a fluid in the following way:

$$\begin{array}{c}\hfill P_{\mathrm{s}}=n_{\mathrm{s}}{k}_{\mathrm{B}}{T}_{\mathrm{s}}.\end{array}$$(31)

Due to the assumptions of thermal equilibrium between ions and electrons (T_i = T_e = T_c) and of charge neutrality of the plasma (n_i = n_e), we have

$$\begin{array}{c}\hfill P_{\mathrm{c}}=n_{\mathrm{c}}{k}_{\mathrm{B}}{T}_{\mathrm{c}}=2n_{\mathrm{i}}{k}_{\mathrm{B}}{T}_{\mathrm{c}}=2P_{\mathrm{i}}=2P_{\mathrm{e}}.\end{array}$$(32)

From the set of equations presented in the previous paragraphs, it can be checked that the only physical dissipative mechanism taken into account is the friction due to collisions between charged and neutral particles. This model does not include other real dissipative processes such as viscosity or resistivity (which, in fact, are different kinds of collisional interactions).

However, the discretisation of these equations that is necessary for solving them numerically leads to the appearance of high-frequency numerical noise. To reduce this unwanted noise and assure the stability of the simulations, a numerical diffusivity term is added to the right-hand side of each of the evolution equations. Details on the implementation of these terms in the MANCHA3D and MANCHA-2F codes can be found in Vögler et al. (2005) and Felipe et al. (2010). In general, each numerical diffusivity term is composed of a shock-resolving term, a hyperdiffusivity part and a constant contribution. In this work, we only consider the latter.

The numerical diffusion terms act as analogues of the viscous forces and the Ohmic dissipation and have a non-negligible impact on the dynamics of the plasma. In Appendix A, we present a brief study of the influence of the magnitude of the numerical dissipation on the results of our simulations.

2.2. Numerical setup

We consider a physical rectangular domain of dimensions [−L₀/2,L₀/2] × [ − L₀, L₀], where L₀ is the length in the x-direction, represented by a numerical mesh of N₀ × 2N₀ points. Periodic boundary conditions are applied in both directions.

This region is filled with a partially ionised plasma whose density is uniform along the horizontal direction but varies in the vertical direction. More precisely, the vertical density profile corresponds to that of a denser slab embedded in a lighter medium, with smooth transitions between the internal and external layers. For the plasma to be in mechanical equilibrium, the pressure of both the charged and the neutral fluids must be uniform. This implies that the initial temperature varies with height in a way that is inversely proportional to the density.

The ionisation degree of the plasma is defined as χ_c = ρ_c/ρ_T, where ρ_T = ρ_c + ρ_n is the total density. In the same way, the density fraction of neutrals is given by χ_n = ρ_n/ρ_T = 1 − χ_c. In general, these parameters vary with the density, pressure, and temperature of the fluids. However, for the sake of simplicity, here we neglect this dependence and assume that χ_c and χ_n are uniform at the initial time.

Then, the vertical density profile of each fluid “s” is given by the following expression:

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{s}}(z,t=0)& =\frac{{\chi }_{\mathrm{s}}}{2}[2{\rho }_{\mathrm{T}1}+({\rho }_{\mathrm{T}2}-{\rho }_{\mathrm{T}1})tanh\left(\frac{z-z_{\mathrm{L}1}}{w_{\mathrm{L}}}\right)\hfill \\ \hfill & +({\rho }_{\mathrm{T}1}-{\rho }_{\mathrm{T}2})tanh\left(\frac{z-z_{\mathrm{L}2}}{w_{\mathrm{L}}}\right)],\hfill \end{array}$$(33)

where ρ_T1 is the total density of the top and bottom regions, ρ_T2 is the total density of the central part, z_L1 and z_L2 are the central positions of the transition layers, and w_L is the width of those layers. For our simulations, we have set the following values of these parameters: ρ_T2 = 2ρ_T1, z_L1 = −L₀/2, z_L2 = L₀/2, and w_L = L₀/25.

As already mentioned, the KHI appears when there is a shear flow at the interface of two media. Therefore, to simulate this instability we consider an initial velocity profile in which the fluids at the central heavier region are moving along the x-direction at a speed U_2, s, while the top and bottom regions move at a speed U_1, s, with a smooth variation between those two values in the transition layers.

Since the evolution of the KHI depends on the velocity drift between the different media, defined as ΔU ≡ |U_1, s − U_2, s|, we can move to a reference frame in which one of the media is at rest. Here, we have chosen the outer layers to be at rest, so U_1, s = 0, while the central region is moving at a speed U_2, s = U₀. Therefore, the initial flow along the x-direction is given by

$$\begin{array}{c}\hfill V_{x,\mathrm{s}}(z,t=0)=\frac{U_0}{2}[tanh\left(\frac{z-z_{\mathrm{L}1}}{w_{\mathrm{L}}}\right)-tanh\left(\frac{z-z_{\mathrm{L}2}}{w_{\mathrm{L}}}\right)].\end{array}$$(34)

Then, we apply a small perturbation to the z-component of the velocity, given by

$$\begin{array}{cc}\hfill V_{z,\mathrm{s}}(x,z,t=0)& =-V_{0}sin\left(\frac{2\pi }{L_0}x\right)exp(-\frac{|z-z_{\mathrm{L}1}|}{0.5w_{\mathrm{L}}})\hfill \\ \hfill & +V_{0}sin\left(\frac{2\pi }{L_0}x\right)exp(-\frac{|z-z_{\mathrm{L}2}|}{0.5w_{\mathrm{L}}}),\hfill \end{array}$$(35)

where the value of V₀ is much smaller than the maximum flow speed of the plasma in the x-direction, that is, V₀ ≪ U₀. This perturbation corresponds to an oscillation of wavelength λ = L₀, with a maximum amplitude at the position of the transition layers and exponentially tends to zero when going away from those positions.

Figure 1 shows vertical cuts of the initial conditions described in the paragraphs above. The top panel shows the normalised temperature, T/T₀ (where T₀ is the temperature of the outer regions), as a dashed green line and the normalised profile of density, Θ, as a solid red line. This parameter is computed as Θ = ρ_s(z)/(χ_sρ_T1). The vertical profile of the normalised velocity in the x-direction, V_x/U₀, is represented at the middle panel. The normalised initial perturbation in the z-component of the velocity, V_z/U₀ is shown at the bottom panel: the solid line corresponds to the position x = −L₀/4, while the dashed line corresponds to x = L₀/4.

Fig. 1. Vertical cuts of the initial conditions of the simulations. Top: normalised density, Θ (red solid line) and normalised temperature (green dashed line) at x = 0. Middle: normalised flow speed at x = 0. Bottom: normalised perturbation of the z-component of the velocity at x = −L0/4 (solid line) and x = L0/4 (dashed line).

Finally, since we are interested in studying how partial ionisation affects the generation of magnetic field by means of the Biermann battery effect, we neglect the presence of an initial background magnetic field. We note that with the setup described in the previous paragraphs, the pressure of each species ‘s’ is not only a function of the fluid density, but it also depends on the temperature, which is non-uniform. It does not fulfill P_s = P_s(ρ_s) but, rather, P_s = P_s(ρ_s, T_s). Therefore, this configuration allows for the action of the Biermann battery effect.

3. Results of the simulations

For the present study, we used values of the initial densities and temperatures typically found in the solar chromosphere. It has been shown, for instance, by Zaqarashvili et al. (2015), Antolin et al. (2018) and Zhelyazkov & Chandra (2019) that motions in chromospheric structures, such as spicules, may lead to the development of the KHI. Using those structures as references, the initial total density and temperature in our simulations are ρ_T1 = 10⁻¹⁰ kg m⁻³ and T₀ = 10⁴ K. Since the sound speed of each fluid ‘s’ is

$$\begin{array}{c}\hfill c_{\mathrm{S},\mathrm{s}}=\sqrt{\frac{\gamma P_{\mathrm{s}}}{{\rho }_{\mathrm{s}}}},\end{array}$$(36)

using Eqs. (31) and (32), the chosen temperature gives values of c_S, c ≈ 16.6 km s⁻¹ and c_S, n ≈ 11.7 km s⁻¹.

In addition, we have set the length-scale to L₀ = 10⁴ m. This length is small in comparison with the thickness or height of the spicules but it is an appropriate value for the investigation of the effects of the charged-neutral interaction.

The works of, for example, Leake et al. (2005), Zaqarashvili et al. (2011), Soler et al. (2013), Martínez-Gómez et al. (2018), and Popescu Braileanu et al. (2019b) have shown that phenomena occurring in a plasma are more affected by ion-neutral collisions when the time and length scales of those phenomena are comparable to the corresponding collisional scales. For example, they have shown that the damping of waves and the plasma heating due to ion-neutral friction are more efficient when the wave frequencies are of the same order than the collision frequencies. On the other hand, the effect of these collisions is negligible when the oscillation frequency is much smaller than the frequency of the ion-neutral interaction.

For the present work, we need to compare the length-scale, L₀, and the numerical resolution, Δx ≡ L₀/N₀, with the collisional mean free path, λ_coll. This parameter can be estimated from Eq. (8) in the following way:

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{coll}}=\frac{V_{\mathrm{ref}}}{{\nu }_{\mathrm{st}}},\end{array}$$(37)

where V_ref is a velocity representative of the dynamics of the system. Assuming the thermal speed given by Eq. (17) as this reference value, the mean free path can be computed as

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{coll}}\approx \frac{m_{\mathrm{s}}+m_{\mathrm{t}}}{n_{\mathrm{t}}{m}_{\mathrm{t}}{\sigma }_{\mathrm{st}}}.\end{array}$$(38)

With a cross-section of σ_in = 1.16 × 10⁻¹⁸ m² for the case of ion-neutral collisions (Leake et al. 2013), we get a mean-free path on the order of 10 to 100 m. This is two orders of magnitude smaller than the length-scale chosen for our simulations. Hence, the dynamics at the large scales should not be greatly affected by collisions. However, the opposite may be expected at the smaller scales since the resolution would be comparable to the mean free path. For instance, if N₀ = 1000, we have that Δx = 10 ≲ λ_coll.

Keeping the values of ρ_T1, T₀, and L₀ mentioned in the previous paragraphs, we performed several series of simulations to study the dependence of the results on various parameters and physical effects. In Sect. 3.1, we present a detailed description of the dynamics of a small set of simulations that show the general evolution of the instability. In the remaining parts of this section we present studies on how the generation of magnetic field is affected by the presence of collisions and by the numerical resolution (Sect. 3.2), and on the dependence on the ionisation degree (Sect. 3.3).

3.1. Reference simulations

We start with a detailed analysis of the case of a partially ionised plasma with an ionisation degree of χ_c = 0.5, that is, with the same values of density for the charged and the neutral fluid. The same flow speed is considered for the two fluids, with the heavier region moving at a speed U₀ = 5 km s⁻¹. This corresponds to Mach numbers of the charged and neutral fluids of

$$\begin{array}{c}\hfill M_{\mathrm{c}}=\frac{U_0}{c_{\mathrm{S},\mathrm{c}}}\approx 0.3\text{and}{M}_{\mathrm{n}}=\frac{U_0}{c_{\mathrm{S},\mathrm{n}}}\approx 0.43,\end{array}$$(39)

respectively. Therefore, the flows are subsonic.

The same perturbation is applied to the z-component of the velocity of both fluids, with an amplitude of V_z, 0 = 10 m s⁻¹. However, this does not mean that charges and neutrals will have exactly the same dynamics. The reason is that, since they start with the same temperature, their pressures are different. This leads to some dissimilarities in their evolution, as we explain later.

3.1.1. Densities

Figure 2 shows the results of a simulation with N₀ = 1000 points, in which the collisional interaction has not been taken into account (α_eff = 0). This means that there is no coupling between the charged and neutral fluids. Density snapshots at different times are presented, where τ ≡ t/t₀ and t₀ is the domain crossing time, defined as t₀ ≡ L₀/U₀. Due to the chosen initial setup, the evolution of the instability is symmetric with respect to z = 0. However, this figure only shows the results for the charged fluid in the upper half of the domain, and those for the neutral fluid in the lower half.

Fig. 2. Density snapshots of a simulation with χc = 0.5 with no coupling between the charged (orange top half) and neutral (blue bottom half) fluids. Time is measured in the units of τ ≡ t/t0, where t0 ≡ L0/U0 is the domain crossing time. Asterisks represent key points whose temporal evolution is analysed in Sect. 3.1.3. An animation of this figure is available online.

At first sight, it seems that the two fluids have exactly the same dynamics: as time advances, the interfaces that separate regions of different densities are distorted because of the transfer of energy between the horizontal and vertical flows. During the initial steps, this distortion has the shape of the sinusoidal perturbation applied to the z-component of the velocities, whose amplitude increases exponentially, according to the linear analysis of Chandrasekhar (1961). At later stages, it turns into a large vortex (Krasny 1988). Then, as energy cascades to smaller scales, secondary vortexes appear and, finally, a turbulent phase is reached (Matsumoto & Seki 2010).

However, if we examine Fig. 2 more closely, some small dissimilarities between the evolution of the charged and the neutral fluids can be found. For example, in the panels corresponding to τ = 6 and τ = 7, it can be seen that the main vortex is slightly more entangled in the case of the charged fluid. In addition, at later times, a greater number of small-scale structures, such as secondary vortexes, can be noticed on the upper half of the plot than in the lower half.

This divergent behaviour exhibited by the two fluids can be understood by taking into account the results of Soler et al. (2012), who studied the effect of compressibility on the KHI. They found that for the case of small density ratios (ρ_T2/ρ_T1), compressibility reduces the growth rate of the instability in comparison with the incompressible limit. This limit is obtained by assuming that the sound speed tends to infinity.

Since in our simulations, c_S, c > c_S, n, compressibility should have a stronger effect on the neutral than in the charged fluid, and the growth rates of the instability should be slightly larger for the latter. This difference on the growth rates explains why the charged fluid is slightly faster in reaching the non-linear and turbulent stages of the KHI.

The existence of the above-mentioned differences is more clearly revealed by defining the variable |Δρ|≡|ρ_c − ρ_n|. We show in the upper half of Fig. 3 snapshots of this variable that correspond to the snapshots shown in Fig. 2. During the linear phase, the differences are very small, as it can be checked by the fact that the panel for τ = 3 is almost empty (small patches of light grey can be found around z/L₀ = 0.5). Later, the differences become more evident and they are clearly concentrated at the edges of the main vortex. But as time advances and turbulence comes into play, they appear in a larger region of the domain.

Fig. 3. Differences in the evolution of density of charged and neutral fluids for the case of χc = χn = 0.5; |Δρ|≡|ρc − ρn|. Top: αeff = 0 (no collisions); bottom: αeff = α (charge-exchange collisions not taken into account). The black lines represent density iso-contours.

The bottom half of Fig. 3 shows the results for a simulation that includes elastic collisions between the charged and neutral fluids (α_eff = α). Compared to the previous case, it can be seen that the dynamics of both fluids are much more similar. Nevertheless, there are still some differences in the evolution of densities that reveal that the coupling between the two fluids is not perfect. These differences have a relative magnitude of less than 0.1, which is notably smaller than the ones for the case without collisions, and are mainly found in the small-scale structures in the cores of the vortexes. On the other hand, the two fluids are strongly coupled in the large scales. These results are in good agreement with the findings of Hillier (2019).

3.1.2. Velocity drifts

Another useful variable to quantify the effect of collisions on the dynamics of the plasma is the velocity drift, v_D. When this variable tends to zero, the two fluids are perfectly coupled and they have the same dynamics, which can be precisely described by single-fluid models. If v_D ≠ 0 (and collisions are considered), there is a transfer of momentum and energy between the two species, as shown by Eqs. (5), (29) and (30).

In Fig. 4, we present snapshots of the x-component of the velocity drift, ΔV_x ≡ V_x,c − V_x,n. The top panels of Fig. 4 show the results for the case without collisions, while the bottom panels correspond to the simulation that includes elastic collisions. As expected, the behaviour shown by the top panels is similar to the one presented in Fig. 3 for the density: during the first stages of the simulation, the main differences in velocity appear around the vortex edges but later they extend to the whole turbulent mixing region. The maximum velocity drifts can be found in the later stages of the evolution and in the smaller-scales, with values of the order of 2.5 km s⁻¹. Over the course of the whole simulation, there are zones in which the charged fluid is moving faster than the neutral one and other in which it is slower. But there is no clear trend on the spatial distribution of these areas.

Fig. 4. Velocity drifts, ΔVx ≡ Vx,c − Vx,n, for the same simulations as in Fig. 3. The black lines represent density iso-contours. We note that the units of the top colourbar are km s−1 and the units of the bottom colourbar are m s−1: the presence of the collisional coupling reduces the magnitude of the velocity drifts.

On the contrary, a trend can be distinguished in the bottom panels of Fig. 4, which represent the case with charged-neutral collisions. Looking to the panel for τ = 5, it can be seen that in the leading part of the main vortex the velocity drifts are negative, which means that the neutral fluid has a larger horizontal speed than the charged fluid. The opposite behaviour is found in the trailing part. This trend is maintained during the remaining of the simulation, even at the turbulent stage. Although we do not include the corresponding plot here, we find a different trend in the z-component of the velocity drifts: in this case, the neutral fluid is moving faster in the upper part of the main vortex and slower in the lower part.

By comparing the bottom panels of Figs. 2 and 4, it can be checked that the largest velocity drifts are found in the regions with the lowest densities. This can be explained by the expression of the momentum transfer. According to Eq. (5), the momentum transfer is proportional to the density of both fluids. Therefore, a lower collisional friction should be expected in regions with lower densities, which leads to a weaker coupling and larger values of the velocity drifts. As shown by the bottom colourbar of Fig. 4, these maximum values of v_D are of the order of 10 m s⁻¹, two orders of magnitude smaller than in the case without collisions. They are also much smaller than the thermal speeds of the two fluids, in agreement with the range of applicability of the momentum and energy transfer expressions given by Eqs. (5), (29), and (30) (Schunk 1977; Draine 1986).

3.1.3. Magnetic field and temperature at key points

Now we pay attention to the temporal evolution of some selected points of the numerical domain. The chosen points are marked as asterisks in Fig. 2. We denote the point at the internal denser region of the plasma as P1, represented as a light blue asterisk. Points P2 and P3, represented by red and green symbols, respectively, are initially located at the transition layers between the internal and external layers of the plasma. Finally, points P4 (black) and P5 (pink) are located at the external lighter regions.

Figure 5 shows the evolution of the absolute value of the magnetic field at points P1, P2, and P4. Qualitatively, the magnetic field at P2 and P4 follows a similar evolution, which can be clearly separated in two stages: exponential growth and saturation (with even a slight decrease). However, during the first stage, at P2 the magnetic field is several orders of magnitude larger than that at the other point. This behaviour is consistent with what is expected from the Biermann battery mechanism: at the transition layers there are larger gradients of pressure and density than in the other regions of the domain. Later, during the saturation stage, the magnetic field at P4 reaches the same order of magnitude than that at P2. This is a consequence of the non-linear evolution of the KHI, which creates a mixing layer that, in this case, contains the two selected points. Then, at the centre of the domain (P1), the magnetic field is negligible in comparison with the other points, although it has a different evolution: for the particular case of this simulation, it continues growing over time without reaching a saturation phase. We checked that in other configurations, with larger values of numerical diffusivity, the saturation stage is indeed reached by all the selected points.

Fig. 5. Absolute value of the magnetic field, |By|, at points P1 (blue dotted line), P2 (red solid line), and P4 (black dashed line) as a function of the normalised time, τ.

The magnetic field at points P3 and P5 has not been represented on Fig. 5 because we have checked that it has the same absolute value as for P2 and P4, respectively, although the actual value has the opposite sign.

The evolution of the temperature at the points P1, P2, and P4 is shown in Fig. 6 (in the same way as for Fig. 5, the results for P3 and P5 are not shown here because they are equivalent to those of P2 and P4, respectively). Initially, the plasma at P1 (the denser region) has a temperature of 5 × 10³ K, while the temperatures at P2 (transition layer) and P4 (lighter region) are ∼6600 K and 10⁴ K, respectively. Over the whole simulation the temperature at P1 is almost the same, with very small variations. That is not the case for P2 and P4. As the instability develops and the plasma is carried by the flow, those static points are crossed by plasma that corresponds to regions with different densities and temperatures, which explains the large fluctuations of temperature at P2. The temperature at P4 remains almost constant for a longer time than for P2 due to the time needed for the development of the mixing layer. Finally, once this mixing layer is created, both points have a similar temperature, which is higher than the initial temperature at P2 but lower than the initial temperature of the external region.

Fig. 6. Temperature as a function of time at points P1 (blue thin line), P2 (red thick line) and P4 (black line). Solid lines represent the temperature of the charged fluid while dashed lines represent the temperature of the neutral fluid. Left and right panels: cases with αeff = 0 and αeff = α.

In addition, the left panel of Fig. 6, which represents a case with no ion-neutral collisions, shows that appreciable differences in the temperatures of the two fluids are present during the turbulent stage. In contrast, when collisions are taken into account, these large temperature differences are removed, as shown by the right panel, in which the solid and dashed lines overlap.

3.2. Magnetic field: growth rate and effect of resolution

In this section, we focus on the evolution of the magnetic field by analysing its average value over the whole domain instead of its evolution at some given points. The average of the magnetic field is computed as follows:

$$\begin{array}{c}\hfill ⟨|\mathit{B}|⟩=\sum _{i=1}^{N_0}\sum _{j=1}^{2N_0}\frac{\sqrt{B_x{(i,j)}^2+B_y{(i,j)}^2+B_z{(i,j)}^2}}{2N_0^2}.\end{array}$$(40)

Figure 7 shows that the spatially averaged magnetic field initially follows a linear phase in which its magnitude increases exponentially with time: from τ = 0 to τ = 4, ⟨|B|⟩ rises up by four orders of magnitude. Then, it continues increasing, but at a lower rate, until it reaches a maximum around τ = 9 and starts decreasing due to the effect of the physical and numerical dissipation mechanisms. We can compute the normalised growth rate of the magnetic field during the linear phase using the following formula:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_B=\frac{ln⟨|\mathit{B}({\tau }_2)|⟩-ln⟨|\mathit{B}({\tau }_1)|⟩}{{\tau }_2-{\tau }_1}.\end{array}$$(41)

Fig. 7. Averaged magnetic field as a function of time and resolution for the cases with χc = 0.5 (left) and χc = 0.1 (right). Solid lines represent simulations in which the ion-neutral collisional interaction is taken into account (without charge-exchange). Dotted lines represent simulations without collisional coupling. Green, blue, purple, and red colours correspond to simulations with N0 = 250, 500, 1000, 2000, respectively. The vertical dashed lines approximately mark the separation between the linear and non-linear phases of the instability. The inset in each panel shows a zoom of the interval τ ∈ [5, 15].

This produces a result of Γ_B ≈ 2.02.

According to Chandrasekhar (1961), the linear growth rate of the KHI (which is expressed in units of s⁻¹) is given by

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{th}}=\frac{2\pi }{\lambda }\mathrm{\Delta }U\frac{\sqrt{{\rho }_{\mathrm{T}1}{\rho }_{\mathrm{T}2}}}{{\rho }_{\mathrm{T}1}+{\rho }_{\mathrm{T}2}},\end{array}$$(42)

which states that perturbations with longer wavelengths, λ, grow slower. In the case of our simulations, the growth rate of the initial perturbation (with λ = L₀) can be expressed in normalised units as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{th}}={\gamma }_{\mathrm{th}}{t}_0=\frac{{\gamma }_{\mathrm{th}}{L}_0}{U_0}=\frac{{\gamma }_{\mathrm{th}}\lambda }{\mathrm{\Delta }U}=2\pi \frac{\sqrt{{\rho }_{\mathrm{T}1}{\rho }_{\mathrm{T}2}}}{{\rho }_{\mathrm{T}1}+{\rho }_{\mathrm{T}2}}\approx 2.96,\end{array}$$(43)

which is larger than the one found in the simulations. This discrepancy can be explained due to the fact that the analytical linear growth rate, Γ_th, was computed assuming a sharp interface between the two media, while in our simulations there is a smooth variation of density and velocity. The smooth transition reduces the growth rate of the instability, as shown by, for instance, Miura & Pritchett (1982), Drazin & Reid (2004) or Berlok & Pfrommer (2019). Using a model in which the velocity varies linearly between the inner and outer layers in a scale of 2w_L, Drazin & Reid (2004) derived an analytic expression for the growth rate for the perturbations in velocity or magnetic field. With the notation employed in the present work, the mentioned expression is given by

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{DR}}=\text{Im}\left[\frac{\mathrm{\Delta }U}{4w_{\mathrm{L}}}\sqrt{{(\frac{4\pi }{\lambda }{w}_{\mathrm{L}}-1)}^2-exp(-\frac{8\pi }{\lambda }{w}_{\mathrm{L}})}\right].\end{array}$$(44)

Inserting the parameters of our simulations in the formula above, we get a result of γ_DR ≈ 1.08. Then, the growth rate for the magnetic field, given in normalised units, would be Γ_DR = γ_DRt₀ ≈ 2.16, which is closer to the value computed from the numerical simulation.

The different colours in Fig. 7 represent simulations with different values of the resolution parameter N₀: 250 (green), 500 (blue), 1000 (purple), and 2000 (red). It can be seen that the results of the four resolutions coincide during the linear phase of the instability. However, they start to diverge once the non-linear stage is reached. As a larger number of points is used and, hence, smaller scales can be better resolved, the generation of magnetic energy is greater and the peak value of magnetic field is reached at a later time. This is mainly due to the effect of the numerical dissipation, which depends on the grid size of the numerical domain. For more details, see Eq. (A.3), which shows how improving the resolution of the simulation decreases the numerical dissipation. Furthermore, it takes longer for the instability to cascade its energy to the dissipation scales in a simulation with a smaller grid size.

In addition, the dotted lines of Fig. 7 correspond to simulations in which the charged-neutral collisions were neglected, and the solid lines correspond to the case with only elastic collisions (no charge-exchange). According to these results, greater magnetic fields are generated when the charged-neutral collisions are taken into account. This effect is studied in more detail in the next section.

3.3. Effect of collisions and dependence on ionisation degree

The comparison of the left and right panels of Fig. 7, which correspond to plasmas with ionisation degrees of χ_c = 0.5 and χ_c = 0.1, respectively, shows that the effect of collisions on the generation of magnetic field varies with the ionisation degree. To investigate such a dependence, we performed four different series of simulations. The details of those simulations are summarised in Table 1. Within each series we have varied the ionisation degree from a minimum value of χ_c = 0.01 to a maximum of χ_c = 0.99. The rest of the parameters, such as the total density, temperatures, and velocities, are the same as in the previous section.

Figure 8 displays three panels that correspond to different stages of the evolution of the instability. At τ = 3 (left panel), the instability is still in the linear phase (as shown in Fig. 7). At the time τ = 6 (middle panel), the magnetic field continues to grow but it does not follow a linear dependence with time. The right panel represents the time τ = 9, in which the plasma is already in a turbulent stage (see the rightmost panel of Fig. 2), and corresponds approximately to the moment when the magnetic energy reaches its maximum value (saturation).

Fig. 8. Spatially averaged magnetic field as a function of ionisation degree at three different times of the simulations. Black asterisks show the results without collisions. Red diamonds correspond to simulations in which ion-neutral collisions are considered. Blue crosses represent the results when the collisional term of the induction equation is also taken into account and the green squares represent the simulations that also included the process of charge-exchange.

The magnitudes of the averaged magnetic field at each stage are quite different (as indicated by the scales in each panel) but the general dependence on the ionisation degree is the same. When no collisions between the charged and the neutral fluids are considered (series I), the resulting magnetic field does not vary with the ionisation degree. The reason is that, in the uncoupled scenario, the generation of magnetic energy is related only to the evolution of the KHI in the charged fluid. Apart from the value of the shear flow velocity and the wavelength of the initial perturbation, which are the same for every simulation in the series, the instability depends on the density ratio between the heavier and lighter media and the sound speed. In our simulations, neither of these parameters are functions of the ionisation degree. So, the same dynamics should be expected throughout this first series.

When the charged-neutral interaction comes into play, Fig. 8 shows a remarkable dependence on the ionisation degree. A greater magnetic field is obtained as the ionisation degree is decreased, that is, as the plasma becomes more weakly ionised. And the variation with respect to the uncoupled case is clearly not negligible. For a partially ionised plasma with χ_c = 0.5, the variation is of about 27%, while for a weakly ionised plasma, with χ_c = 0.01, it reaches up to ∼80%.

Comparing the results of the three series that included collisional terms, it can be seen that the effect of the collisional term of the induction equation is negligible. The results for the case that includes that term (series III) almost overlap with the results corresponding to the case with only elastic collisions (series II). On the other hand, the inclusion of charge-exchange collisions (series IV) slightly increases the averaged magnetic field, specially in the weakly ionised regime. The physical reasons behind the results described in the paragraphs above are explained in Sect. 4.

4. Discussion: Analysis of the equations

To better understand the results of the simulations described in the previous section, it is necessary to analyse in greater detail the equations included in the two-fluid model. In the present section, we pay special attention to the Biermann battery term of the induction equation and the collisional terms that appear in the momentum and the energy equations.

4.1. Induction equation

First, we analyse the case without the collisional interaction between charged and neutral particles. Since our simulations start without an initial seed, the very first steps of the evolution of the magnetic field are described by Eq. (24), which can also be written as

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=-\frac{m_{\mathrm{i}}\mathrm{\nabla }{\rho }_{\mathrm{c}}\times \mathrm{\nabla }{P}_{\mathrm{c}}}{2e{\rho }_{\mathrm{c}}^2}=-\frac{k_{\mathrm{B}}\mathrm{\nabla }{\rho }_{\mathrm{T}}\times \mathrm{\nabla }{T}_{\mathrm{c}}}{e{\rho }_{\mathrm{T}}}.\end{array}$$(45)

Here, we have used the facts that ρ_c = χ_cρ_T, P_c = 2ρ_c/m_ik_BT_c and χ_c is uniform. With this formula, we can explain the behaviour presented in Fig. 8 for the series without collisions (black asterisks). Under the conditions of that series, the evolution of the magnetic field is independent from the ionisation degree of the plasma.

In addition, if we take into account that in our model there are no variations along the y-direction, that is, ∂y = 0, from Eq. (45) we obtain the result of the battery term only generating magnetic field in the y-direction, thus, B_x = B_z = 0. The evolution of the y-component of the magnetic field is given by

$$\begin{array}{c}\hfill \frac{\partial B_y}{\partial t}=\frac{m_{\mathrm{i}}}{2e{\rho }_{\mathrm{c}}^2}[\left({\partial }_x{\rho }_{\mathrm{c}}\right)\left({\partial }_z{P}_{\mathrm{c}}\right)-\left({\partial }_z{\rho }_{\mathrm{c}}\right)\left({\partial }_x{P}_{\mathrm{c}}\right)]\end{array}$$(46)

where ∂_x and ∂_z represent the derivatives with to x and to z, respectively.

Once a seed magnetic field appears due to the action of the Biermann battery term, the advection term of the induction equation becomes different from zero. Taking also into account that V_y, c = 0, the evolution of the magnetic field due to that term is given by

$$\begin{array}{c}\hfill \frac{\partial B_y}{\partial t}=-{\partial }_x\left(V_{x,\mathrm{c}}{B}_y\right)-{\partial }_z\left(V_{z,\mathrm{c}}{B}_y\right).\end{array}$$(47)

The fact that no magnetic field is present in the x and z directions is important for the evolution of the instability. It is known that a magnetic field that is oriented along the direction of the shear flow velocity has a stabilizing effect (Chandrasekhar 1961). In this case, the generated magnetic field is perpendicular to that flow, so it has no effect on the KHI.

4.2. Vorticity

The vorticity of any species ‘s’ is defined as the curl of its velocity, that is, ω_s ≡ ∇ × V_s. By operating the continuity and momentum equations, Eqs. (1)–(3), we obtain equations of evolution for the vorticities of both fluids:

$$\begin{array}{cc}\hfill \frac{\partial {\mathit{\omega }}_{\mathrm{c}}}{\partial t}& =\mathrm{\nabla }\times ({\mathit{V}}_{\mathrm{c}}\times {\mathit{\omega }}_{\mathrm{c}})+\frac{\mathrm{\nabla }{\rho }_{\mathrm{c}}\times \mathrm{\nabla }{P}_{\mathrm{c}}}{{\rho }_{\mathrm{c}}^2}\hfill \\ \hfill & +\frac{1}{{\mathrm{\mu }}_0}\mathrm{\nabla }\times \left[\frac{(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}}{{\rho }_{\mathrm{c}}}\right]+\mathrm{\nabla }\times \left[{\nu }_{\mathrm{cn}}({\mathit{V}}_{\mathrm{n}}-{\mathit{V}}_{\mathrm{c}})\right]\hfill \end{array}$$(48)

and

$$\begin{array}{cc}\hfill \frac{\partial {\mathit{\omega }}_{\mathrm{n}}}{\partial t}& =\mathrm{\nabla }\times ({\mathit{V}}_{\mathrm{n}}\times {\mathit{\omega }}_{\mathrm{n}})+\frac{\mathrm{\nabla }{\rho }_{\mathrm{n}}\times \mathrm{\nabla }{P}_{\mathrm{n}}}{{\rho }_{\mathrm{n}}^2}\hfill \\ \hfill & +\mathrm{\nabla }\times \left[{\nu }_{\mathrm{nc}}({\mathit{V}}_{\mathrm{c}}-{\mathit{V}}_{\mathrm{n}})\right].\hfill \end{array}$$(49)

It can be seen that the vorticity equation for the charged fluid has a great resemblance with the induction equation. So, in the same way, we rewrite it as

$$\begin{array}{c}\hfill \frac{\partial {\mathit{\omega }}_{\mathrm{c}}}{\partial t}={\mathit{A}}_{\omega \mathrm{c}}+{\mathit{B}}_{\omega \mathrm{c}}+{\mathit{L}}_{\omega \mathrm{c}}+{\mathit{C}}_{\omega \mathrm{c}},\end{array}$$(50)

where A_ωc is the advective term, B_ωc is called the baroclinic term, L_ωc is related to the Lorentz force, and C_ωc refers to the collisional interaction between the charged and neutral fluids. The same procedure can be applied to the neutral vorticity equation, with the difference that it does not contain a term related to the magnetic field:

$$\begin{array}{c}\hfill \frac{\partial {\mathit{\omega }}_{\mathrm{n}}}{\partial t}={\mathit{A}}_{\omega \mathrm{n}}+{\mathit{B}}_{\omega \mathrm{n}}+{\mathit{C}}_{\omega \mathrm{n}}.\end{array}$$(51)

4.3. Comparison of induction and vorticity terms

By comparing the induction equation, Eq. (23), and the vorticity equation for the charged fluid, Eq. (48), it can be checked that several of their terms share the dependence on the magnetic field, densities, pressures, and velocities, but differ in terms of their coefficients. The following relations are found:

$$\begin{array}{c}\hfill \frac{|{\mathit{B}}_{\mathrm{Ind}}|}{|{\mathit{B}}_{\omega \mathrm{c}}|}=\frac{m_{\mathrm{i}}}{2e}\sim 5\times {10}^{-9}\mathrm{kg}{\mathrm{C}}^{-1},\end{array}$$(52)

and

$$\begin{array}{c}\hfill \frac{|{\mathit{C}}_{\mathrm{Ind}}|}{|{\mathit{C}}_{\omega \mathrm{c}}|}\sim \frac{m_{\mathrm{e}}}{e}\sim 5\times {10}^{-12}\mathrm{kg}{\mathrm{C}}^{-1}.\end{array}$$(53)

If we compare Eqs. (52) and (53), we get

$$\begin{array}{c}\hfill \frac{|{\mathit{C}}_{\mathrm{Ind}}|}{|{\mathit{C}}_{\omega \mathrm{c}}|}\ll \frac{|{\mathit{B}}_{\mathrm{Ind}}|}{|{\mathit{B}}_{\omega \mathrm{c}}|}.\end{array}$$(54)

The collisional term is much less important in the evolution of the magnetic field than in the evolution of vorticity in comparison with the battery and the baroclinic terms. This explains the results shown in Fig. 8 for the series of simulations III (blue crosses): the effect of the collisional term in the induction equation is negligible.

Now, performing a dimensional analysis of the baroclinic terms from the charged and neutral vorticity equations we get

$$\begin{array}{c}\hfill |{\mathit{B}}_{\omega \mathrm{c}}|\sim \frac{1}{\gamma }{\left(\frac{c_{\mathrm{S},\mathrm{c}}}{L_{\mathrm{ref}}}\right)}^2\text{and}|{\mathit{B}}_{\omega \mathrm{n}}|\sim \frac{1}{\gamma }{\left(\frac{c_{\mathrm{S},\mathrm{n}}}{L_{\mathrm{ref}}}\right)}^2,\end{array}$$(55)

where L_ref is a reference value of the length scales of the system. Furthermore, if we compare both terms and assume that the two fluids have the same temperature, we get the next relation:

$$\begin{array}{c}\hfill \frac{|{\mathit{B}}_{\omega \mathrm{c}}|}{|{\mathit{B}}_{\omega \mathrm{n}}|}\sim \frac{c_{\mathrm{S},\mathrm{c}}^2}{c_{\mathrm{S},\mathrm{n}}^2}=2,\end{array}$$(56)

which implies that the baroclinic term tends to produce a greater vorticity in the charged fluid than in the neutral one, due to the greater sound speed of the former. This relation, in combination with Eq. (52), is relevant for the explanation of the results from the simulations that we present in Sect. 4.4.

Figure 9 shows the average values of the baroclinic terms of the charged and neutral vorticity equations as functions of the ionisation degree at the times τ = 3 (left), τ = 6 (middle), and τ = 9 (right). These results can be compared with Fig. 8. We recall that |B_Ind|=|B_ωc|m_i/(2e); thus, the dependence depicted here for the charged baroclinic term is also applicable to the battery term.

Fig. 9. Spatially averaged values of the baroclinic terms as functions of the ionisation degree, χc, at three different times: τ = 3 (left), τ = 6 (middle), and τ = 9 (right). Solid and dashed lines correspond to the terms ⟨|Bωc|⟩ and ⟨|Bωn|⟩, respectively. The black, red and green colours correspond to the cases with αeff = 0, αeff = α, and αeff = α + αcx, respectively.

In the case where there is no coupling between the two fluids, Fig. 9 shows that the baroclinic terms do not change with the ionisation degree. This agrees with the prediction of Eq. (45) for the battery term. However, when collisions are considered, each baroclinic term has a different behaviour: although both display a similar dependence on the ionisation degree (that is, their magnitudes decrease when χ_c increases), B_ωc is larger than in the uncoupled case whereas B_ωn is smaller. The neutral baroclinic term (dashed lines) is less affected by the interaction with the charged fluid in the weakly ionised regime, while the charged baroclinic term (solid lines) is less affected by the interaction with the neutral fluid in the strongly ionised limit.

In comparing Figs. 8 and 9, it can be seen that the magnetic field follows the same trends as the charged fluid baroclinic term or, equivalently, the battery term. However, this term does not directly depend on any kind of collisional parameter. So, the influence of the collisional interaction on the generation of magnetic field must come indirectly from the dependence of the Biermann battery term on the density or on the pressure of the charged fluid. The continuity equation that describes the temporal evolution of the density, Eq. (1), does not involve any term related to collisions. On the other hand, the internal energy equations, Eqs. (25) and (26) do include source terms on their right-hand sides that represent the heat transfer due to the collisional interaction between the two fluids. We analyse those equations in the following section.

4.4. Pressure, temperature and energy transfer

Inserting in Eqs. (25) and (26) the relations between the internal energies and the pressures given by Eq. (28), we obtain the following expressions for the temporal evolution of the charged and neutral pressures:

$$\begin{array}{c}\hfill \frac{\partial P_{\mathrm{c}}}{\partial t}+{\mathit{V}}_{\mathrm{c}}·\mathrm{\nabla }{P}_{\mathrm{c}}+\gamma P_{\mathrm{c}}\mathrm{\nabla }·{\mathit{V}}_{\mathrm{c}}=(\gamma -1)[\mathit{J}·{\mathit{E}}_{\mathrm{diff}}+Q_{\mathrm{cn}}]\end{array}$$(57)

and

$$\begin{array}{c}\hfill \frac{\partial P_{\mathrm{n}}}{\partial t}+{\mathit{V}}_{\mathrm{n}}·\mathrm{\nabla }{P}_{\mathrm{n}}+\gamma P_{\mathrm{n}}\mathrm{\nabla }·{\mathit{V}}_{\mathrm{n}}=(\gamma -1)Q_{\mathrm{nc}},\end{array}$$(58)

with Q_cn and Q_nc given by Eqs. (29) and (30), respectively. We recall here that Q_cn ≠ Q_nc but they follow the relation

$$\begin{array}{c}\hfill Q_{\mathrm{cn}}+Q_{\mathrm{nc}}={\alpha }_{\mathrm{eff}}{\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}{v}_{\mathrm{D}}^2,\end{array}$$(59)

which is equivalent to the expression Q_ab + Q_ba = −R_ab ⋅ (V_a−V_b) from Braginskii (1965).

For the present analysis, it is useful to express the energy transfer terms as $Q_{\mathrm{cn}}=Q_{\mathrm{cn}}^T+Q_{\mathrm{cn}}^V$ and $Q_{\mathrm{nc}}=Q_{\mathrm{nc}}^T+Q_{\mathrm{nc}}^V$, where

$$\begin{array}{c}\hfill Q_{\mathrm{cn}}^T=\frac{{\alpha }_{\mathrm{eff}}{\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}{k}_{\mathrm{B}}}{(\gamma -1)m_{\mathrm{n}}}(T_{\mathrm{n}}-T_{\mathrm{c}})=-Q_{\mathrm{nc}}^T\end{array}$$(60)

are the terms depending on the difference of temperatures, and then we have

$$\begin{array}{c}\hfill Q_{\mathrm{cn}}^V=Q_{\mathrm{nc}}^V={\alpha }_{\mathrm{eff}}{\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}\frac{v_{\mathrm{D}}^2}{2}\end{array}$$(61)

as the terms related to the velocity drifts. Since $Q_{\mathrm{cn}}^V$ and $Q_{\mathrm{nc}}^V$ are always positive, their effect is to increase the pressures (or internal energies) of both fluids, but the thermal exchange terms increase the pressure of one fluid while decreasing the pressure of the other one.

Compared to the neutral pressure equation, the equation for the charges pressure includes Joule heating term that depends on the current density. However, due to the small values of the magnetic field present in the simulations, this term can be neglected in comparison with the collisional term, as shown in Appendix C.

Furthermore, according to the setup described in Sect. 2.2, the initial temperatures of both fluids are the same, T_c = T_n, and $Q_{\mathrm{cn}}^T=0=Q_{\mathrm{nc}}^T$ initially. To check whether the thermal exchange terms remain negligible during the evolution as well, we now shift our focus to the equations for the temperatures.

From Eqs. (57) and (58), with the assistance of the continuity equation, Eq. (1), and the equation of state, Eq. (31), the following expression is obtained for each fluid s ∈ {c, n}:

$$\begin{array}{c}\hfill \frac{\partial T_{\mathrm{s}}}{\partial t}+{\mathit{V}}_{\mathrm{s}}·\mathrm{\nabla }{T}_{\mathrm{s}}+(\gamma -1)T_{\mathrm{s}}\mathrm{\nabla }·{\mathit{V}}_{\mathrm{s}}={\left(\frac{\partial T_{\mathrm{s}}}{\partial t}\right)}_{\mathrm{coll}},\end{array}$$(62)

with

$$\begin{array}{c}\hfill {\left(\frac{\partial T_{\mathrm{c}}}{\partial t}\right)}_{\mathrm{coll}}=\frac{(\gamma -1)}{n_{\mathrm{c}}{k}_{\mathrm{B}}}[Q_{\mathrm{cn}}^T+Q_{\mathrm{cn}}^V]\end{array}$$(63)

and

$$\begin{array}{c}\hfill {\left(\frac{\partial T_{\mathrm{n}}}{\partial t}\right)}_{\mathrm{coll}}=\frac{(\gamma -1)}{n_{\mathrm{n}}{k}_{\mathrm{B}}}[-Q_{\mathrm{cn}}^T+Q_{\mathrm{cn}}^V].\end{array}$$(64)

When there are no collisions, the temperature of each fluid evolves independently. Given that the velocities of neutrals and charges are not the same (see Fig. 4), it is expected that differences in the temperatures will appear during the evolution.

The charged-neutral collisional interaction produces temporal variations of temperature that depend on the number density of each fluid, that is, they depend on the ionisation degree of the plasma. As a consequence of $Q_{\mathrm{cn}}^V$, the temperature of neutrals and charges will increase, but it will do in different amounts because the heating rate is inversely proportional to their corresponding number density. In contrast, $Q_{\mathrm{cn}}^T$ can increase or decrease the temperatures of neutrals and charges, depending on the sign of the difference of temperature between them. The equations for the temperature evolution under the influence of the solely $Q_{\mathrm{cn}}^T$ term show that the temperatures tend to equilibrate on the time scale given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{col}}=\frac{1}{{\nu }_{\mathrm{cn}}+{\nu }_{\mathrm{nc}}},\end{array}$$(65)

as already shown by previous studies on partially ionised plasmas. Therefore, the net effect of the charged-neutral collisions is to increase the temperature of the two fluids by means of the dissipation of kinetic energy and then to lead the plasma to a new thermal equilibrium by means of the thermal exchange (Leake et al. 2005; Zaqarashvili et al. 2011; Soler et al. 2013; Oliver et al. 2016; Martínez-Gómez et al. 2018).

Coming back to the baroclinic terms in the vorticity equations, in the situation where T_c(x, z)≈T_n(x, z) (and consequently, ∇T_c ≈ ∇T_n) holds, and also ∇ρ_c/ρ_c ≈ ∇ρ_n/ρ_n is assumed¹, the relation |B_ωc|∼2|B_ωn| has to be valid for any ionisation degree, see Eq. (56). According to Fig. 9, the neutral and charges baroclinic terms, |B_ωn| and |B_ωc|, do vary with ionisation degree χ_c. They both decrease when this parameter is increased. However, the ratio between both, |B_ωc|/|B_ωn|, is stable at approximately 2, which is in agreement with our order of magnitude estimates above. In the weakly ionised limit, χ_c → 0, the values of |B_ωn| for the cases with and without collisions are almost equal. In the opposite limit of high ionisation, χ_c → 1, the charged baroclinic terms, |B_ωc|, for the cases with and without collisions become equal.

The reason for the dependence of the strength of the magnetic field generated by the battery term on the ionisation degree can now be understood taking into account the results from Fig. 9 and the behaviour of the collisional terms in Eq. (62). We consider in the first place a plasma that is almost fully ionised, with n_c ≫ n_n. In that case, we find that

$$\begin{array}{c}\hfill \left|{\left(\frac{\partial T_{\mathrm{c}}}{\partial t}\right)}_{\mathrm{coll}}\right|\ll \left|{\left(\frac{\partial T_{\mathrm{n}}}{\partial t}\right)}_{\mathrm{coll}}\right|,\end{array}$$(66)

which means that the temperature of the charged fluid is almost unaffected by the interaction with the neutral fluid. The latter will be subject to much greater temporal variations in temperature. The baroclinic term of the charged fluid will remain almost the same as in the case without collisions. So, for the relation given by Eq. (56) to hold, |B_n| is reduced.

On the other hand, in the weakly ionised regime, that is n_c ≪ n_n, we get

$$\begin{array}{c}\hfill \left|{\left(\frac{\partial T_{\mathrm{c}}}{\partial t}\right)}_{\mathrm{coll}}\right|\gg \left|{\left(\frac{\partial T_{\mathrm{n}}}{\partial t}\right)}_{\mathrm{coll}}\right|,\end{array}$$(67)

so now the charged fluid has much greater variations of temperature. This leads to an increase in the magnitude of the charged baroclinic term in comparison with the case without collisions. Consequently, the battery term, which is proportional to |B_ωc|, also increases and produces stronger magnetic fields.

The behaviour of the collisional terms described above can be checked in Fig. 10, where we represent the average of their absolute values at different times of the simulations. By comparing Figs. 9 and 10, we can see that the baroclinic term for the charged fluid and the collisional term given by Eq. (63) (represented by the solid red lines in the respective figures) follow similar trends – increasing as the ionisation degree of the plasma decreases.

Fig. 10. Spatially averaged values of the collisional terms in the temperature equations as functions of the ionisation degree, at three different times: τ = 3 (left), τ = 6 (middle), and τ = 9 (right). Solid and dashed lines represent the terms for the charged and the neutral fluids, respectively, which are given by Eqs. (63) and (64).

4.5. Comparison with a single-fluid model

To study the origin of cosmic magnetic fields, Kulsrud et al. (1997) used a single-fluid model to simulate the generation of magnetic fields by the Biermann battery term in partially ionised plasmas. With the notation of the present work, the induction equation of that single-fluid model can be written as

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=\mathrm{\nabla }\times [\mathit{V}\times \mathit{B}]-\frac{m_{\mathrm{i}}}{e}\frac{1}{1+{\chi }_{\mathrm{c}}}\frac{\mathrm{\nabla }{\rho }_{\mathrm{T}}\times \mathrm{\nabla }{P}_{\mathrm{T}}}{{\rho }_{\mathrm{T}}^2}.\end{array}$$(68)

Here, the variables V and P_T are the weighted mean velocity of the partially ionised plasma and the total pressure, and are given by

$$\begin{array}{c}\hfill \mathit{V}=\frac{{\rho }_{\mathrm{c}}{\mathit{V}}_{\mathrm{c}}+{\rho }_{\mathrm{n}}{\mathit{V}}_{\mathrm{n}}}{{\rho }_{\mathrm{T}}}={\chi }_{\mathrm{c}}{\mathit{V}}_{\mathrm{c}}+(1-{\chi }_{\mathrm{c}}){\mathit{V}}_{\mathrm{n}}\end{array}$$(69)

and

$$\begin{array}{c}\hfill P_{\mathrm{T}}=P_{\mathrm{c}}+P_{\mathrm{n}}=P_{\mathrm{e}}+P_{\mathrm{i}}+P_{\mathrm{n}}.\end{array}$$(70)

The behaviour of the magnetic field as a function of the ionisation degree predicted by Eq. (68) is similar to that depicted in Fig. 8. In the fully ionised limit (χ_c = 1), Eq. (68) is equivalent to the first line of Eq. (23). In addition, both the single-fluid model and the two-fluid model predict that when the ionisation degree is decreased, greater magnetic fields are generated. Nevertheless, the single-fluid model predicts a greater increase in the magnetic field. In the weakly ionised regime, when χ_c → 0, the battery term in Eq. (68) is a factor ∼2 larger than in the fully ionised case. On the other hand, Fig. 8 displays a maximum ratio between the two limits of ∼1.8 when charge-exchange collisions are considered or ∼1.7 when they are not. These discrepancies can be understood in terms of the assumptions that each model makes about the coupling degree between the two species.

The single-fluid model assumes a strong coupling between the two species, which have the same temperature and almost negligible velocity drifts, v_D. In contrast, in the two-fluid model the strength of the coupling is not assumed a priori but computed in each time step. It is represented by the factor α_effρ_cρ_n, which greatly depends on the total density of the plasma and the ionisation degree, as clearly shown when expressed in the following way:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{eff}}{\rho }_{\mathrm{c}}{\rho }_{\mathrm{n}}={\alpha }_{\mathrm{eff}}{\chi }_{\mathrm{c}}(1-{\chi }_{\mathrm{c}}){\rho }_{\mathrm{T}}^2.\end{array}$$(71)

For a given value of the total density, the largest friction coefficient is obtained when χ_c = 0.5, while it obviously tends to zero in the fully ionised and the fully neutral limits. Therefore, in the weakly ionised regime the coupling between the two species may be not as strong as it is assumed by the single-fluid model.

To check how the coupling and, consequently, the generation of magnetic field depend on the total density, we performed additional series of simulations with three different values of ρ_T. The results of this study are presented in Fig. 11, in which the spatially averaged magnetic field is plotted as a function of the ionisation degree at the time τ = 9. The blue dotted, grey dashed-dotted, and red dashed lines represent the results for the cases with ρ_T = 10⁻¹¹,  10⁻¹⁰, and 10⁻⁹, respectively, in which α_eff = α. The symbols correspond to simulations with no charged-neutral collisions. The theoretical prediction of the single-fluid model is included as a solid black line and has been computed by multiplying the results with α_eff = 0 by a factor 2/(1 + χ_c). It can be seen that the results from the two-fluid simulations approach the behaviour predicted by the single-fluid model as the total density of the plasma increases.

Fig. 11. Spatially averaged magnetic field as a function of ionisation degree for simulations with N0 = 500 and different values of total density, ρT. Symbols and lines correspond to simulations with αeff = 0 and αeff = α, respectively. The blue dotted, grey dashed-dotted, and red dashed lines represent the cases with ρT = 10−11,  10−10 and 10−9 kg m−3, respectively. The black line corresponds to the prediction of the single-fluid model.

It is also remarkable that the results presented in Fig. 11 for the case with no collisions do not depend on the total density of the plasma. This is explained by the previously noted fact that the evolution of the KHI does not depend on the total density but, rather, on the ratio of densities of the different media.

5. Summary and future works

In this paper, we investigated the effect that the partial ionisation in plasmas has on the Biermann battery mechanism. We used the two-fluid model implemented in the numerical code MANCHA-2F (Popescu Braileanu et al. 2019a) to simulate the evolution of the Kelvin-Helmholtz instability in plasmas with different ionisation degrees. From the initial conditions with no magnetic field, we studied how a weak magnetic field is created by means of the battery mechanism and how it evolves as the instability develops.

First, we compared the evolution of the charged and the neutral fluids in the cases without and with coupling due to elastic collisions, for a particular value of the ionisation degree. We found that when no collisions are considered, the KHI develops slightly faster in the charged fluid than in the neutral one, although the initial conditions of density, velocity, and temperature were identical for both of them. The explanation for this behaviour is that the charged and neutral fluids have different pressures (the charged fluid takes into account the pressure of electrons). As a consequence, the sound speed of the charged fluid is greater. With the density ratio between the internal and external regions of the plasma slabs chosen in our simulations, a greater sound speed leads to a slightly higher growth rate of the instability (Soler et al. 2012). In the case with collisions, the two fluids follow the same evolution at the large spatial scales but differences still appear at the smallest scales, which is in good agreement with the results of Hillier (2019).

Then we found that the collisional interaction between the charged and neutral fluids enhances the generation of magnetic field by means of the Biermann battery mechanism (see Figs. 7 and 8). We investigated how the generation and amplification of the magnetic field depends on the ionisation degree of the plasma and what role different collisional mechanisms play in this processes.

On the one hand, we found that if there is no charged-neutral coupling, the generation of magnetic does not depend on the ionisation degree. The reason behind this is that the evolution of the KHI does not depend on the total density of the charged fluid but on the density ratio between the different regions of the plasma, and this ratio has been kept the same for all the simulations.

On the other hand, when collisions are taken into account, we found that the generated magnetic field increases as the ionisation degree of the plasma decreases. The main reason for the increase in the magnetic field with decreasing ionisation degree comes from the elastic collisional terms included in the momentum and energy equations. The induction equation also contains a term related to elastic collisions but its contribution is negligible. Finally, we also considered the process of charge-exchange collisions and we find that it slightly enhances the generated magnetic field, compared to the case with only elastic collisions.

In order to find an explanation for the behaviour described above, we carried out an analysis of the equations used in the two-fluid model. The Biermann battery term does not have any explicit dependence on the ionisation degree of the plasma or any parameter related to collisions. However, the simulations presented in this work show that the magnetic field generated by this mechanism does depend on the ionisation degree when the interaction between charges and neutrals is taken into account. Our analysis has shown that this dependence comes from the fact that the battery term is a function of the gradients of pressure or temperature and these variables are greatly affected by the presence of collisions.

The charged-neutral interaction by means of collisions exerts two effects on the temperatures of the plasma. On the one hand, it increases the temperature of the two fluids by dissipating kinetic energy and transforming it into internal energy. The resulting increase of temperature is not the same for both fluids but it also depends on their number density, that is, on the ionisation degree. On the other hand, if the temperatures of the two fluids are different, collisions tend to reduce that difference until a thermal equilibrium is reached, with T_c ∼ T_n. Under this equilibrium, the baroclinic terms from the vorticity equations tend to fulfill the relation |B_ωc|∼2|B_ωn|, in contrast to the relation |B_ωc|∼|B_ωn| found when there is no coupling (as shown by Fig. 9). Collisions produce smaller variations of temperature in the fluid with a greater density. Therefore, the magnitude of the baroclinic term of that fluid will have approximately the same value as in the case with no coupling. Then, in order for the relation |B_ωc|∼2|B_ωn| to be fulfilled, the magnitude of the other baroclinic term must be modified accordingly. This is why |B_ωc| increases as the ionisation degree decreases. In addition, since the Biermann battery term is proportional to the baroclinic term of the charged fluid, the generation of magnetic field is also enhanced as the ionisation degree is reduced.

Finally, we compared our numerical results with the predictions from the single-fluid model used by Kulsrud et al. (1997). The single-fluid and the two-fluid model agree in the qualitative aspect: they both show that the generation of magnetic field by the Biermann battery mechanism is enhanced as the ionisation degree of the plasma is reduced. However, they differ in the magnitude of that enhancement. The single-fluid predicts a greater increase, which does not depend on the total density of the plasma. On the contrary, the results of the two-fluid model do depend on the total density because they depend on the coupling degree between the two species. This parameter is a function of the ionisation degree and the total density. As the total density is increased, the two-fluid results show a better agreement with the single-fluid predictions. The results presented in Fig. 11 show that the trend of an increase in the magnetic field as the ionisation degree decreases holds in the range χ_c ∈ [10⁻³, 1]. However, according to Eq. (71), the coupling between the two fluids should tend to zero as the ionisation degree tends to zero. And without coupling there should not be enhancement of the Biermann mechanism. For future works, it would be interesting to extend the present study to even lower values of ionisation degree and check whether the mentioned trend still holds.

For the simulations analysed in the present work, we have chosen values of density, temperatures, velocities, and length-scales that correspond to plasmas in the solar chromosphere. However, our results may easily be extrapolated and prove relevant in other astrophysical scenarios in which partially ionised plasmas are involved. The time and length scales of the evolution of the KHI and the strength of the magnetic field generated by the battery term would be different, but its qualitative behaviour regarding the interaction between the charged and neutral components of the plasma would be the same.

There are some results from our simulations that we briefly mentioned in this work but these areas have been left unexplored in greater depth, such as certain trends observed in the velocity drifts when the collisional coupling is considered (see Fig. 4). This issue should be analysed in more detail in future works.

Furthermore, here we have assumed that the initial ionisation degree is uniform while, in general, it should be a function of the temperature and density of the plasma. In addition, it would be interesting to study the dependence of the results on the initial conditions of temperature. A higher temperature corresponds to a larger value of the sound speed, which brings the plasma closer to the incompressible regime. With the setup chosen for our simulations, and according to the results of Soler et al. (2012), this would produce higher growth rates of the KHI. Therefore, stronger magnetic fields would be generated during the first stages of the instability. However, a higher growth rate also means that the dissipation scales are reached more quickly and the maximum magnetic field obtained may be smaller than in a case with a lower temperature.

Another important step would be to perform 3D simulations. The values of the magnetic field obtained in our 2D simulations are very small in comparison to those that are found in reality. It is expected that a 3D configuration would lead to a considerably greater amplification of the magnetic field since it would not be restricted to only one spatial direction – as in the 2D simulations of the present work. In addition, some relevant physical processes that have been neglected here, such as the real viscosity of the fluids or the thermal conduction, should be taken into account in future works and thus provide more realistic results.

Movies

Movie 1 associated with Fig. 2 Access here

Movie 2 associated with Fig. A.2 Access here

Acknowledgments

This work was supported by the European Research Council through the Consolidator Grant ERC-2017-CoG-771310-PI2FA. We thankfully acknowledge the technical expertise and assistance provided by the Spanish Supercomputing Network (Red Espanola de Supercomputacion), as well as the computer resources used: the LaPalma Supercomputer, located at the Instituto de Astrofisica de Canarias. We also thank the anonymous referee for constructive remarks and helpful suggestions.

References

Appendix A: Effect of numerical dissipation

In the MANCHA-2F code, the numerical diffusion terms are implemented in such a way that they replicate the effects of viscosity and resistivity. This means that, for instance, the following term is added to the right-hand side of the momentum equation, Eq. (2):

$$\begin{array}{c}\hfill {\left(\frac{\partial \left({\rho }_{\mathrm{c}}{\mathit{V}}_{\mathrm{c}}\right)}{\partial t}\right)}_{\mathrm{art}}=\mathrm{\nabla }·{\widehat{\tau }}_{\mathrm{c},\mathrm{art}},\end{array}$$(A.1)

where ${\widehat{\tau }}_{\mathrm{c},\mathrm{art}}$ is the numerical viscous stress tensor. The components of this tensor are computed as follows:

$$\begin{array}{c}\hfill {\widehat{\tau }}_{\mathrm{c},\mathrm{art},\mathrm{kl}}=\frac{1}{2}{\rho }_{\mathrm{c}}({\nu }_k\frac{\partial V_{\mathrm{c},\mathrm{l}}}{\partial x_k}+{\nu }_l\frac{\partial V_{\mathrm{c},\mathrm{k}}}{\partial x_l}),\end{array}$$(A.2)

with

$$\begin{array}{c}\hfill {\nu }_l=a_{\mathrm{diff}}{c}_{\mathrm{S},\mathrm{c}}\mathrm{\Delta }{x}_l.\end{array}$$(A.3)

Here, a_diff is a constant that controls the strength of the numerical diffusion, c_S, c is the sound speed of the charged fluid and Δx_l is the resolution of the numerical domain.

On the other hand, to take into account the effect of the physical viscosity, the right-hand side of the momentum equation for each species s should include the term

$$\begin{array}{c}\hfill {\left(\frac{\partial \left({\rho }_{\mathrm{s}}{\mathit{V}}_{\mathrm{s}}\right)}{\partial t}\right)}_{\mathrm{visc}}=\mathrm{\nabla }·{\widehat{\tau }}_{\mathrm{s}},\end{array}$$(A.4)

where

$$\begin{array}{c}\hfill {\tau }_{\mathrm{s},ij}={\xi }_{\mathrm{s}}(\frac{\partial V_{\mathrm{s},i}}{\partial x_j}+\frac{\partial V_{\mathrm{s},j}}{\partial x_i}-\frac{2}{3}{\delta }_{\mathit{ij}}\mathrm{\nabla }·{\mathit{V}}_{\mathrm{s}})\end{array}$$(A.5)

and the viscosity ξ_s is computed as (Braginskii 1965; Leake et al. 2013)

$$\begin{array}{c}\hfill {\xi }_{\mathrm{s}}=\frac{n_{\mathrm{s}}{k}_{\mathrm{B}}{T}_{\mathrm{s}}}{{\nu }_{\mathrm{ss}}}=\frac{\sqrt{\pi k_{\mathrm{B}}{T}_{\mathrm{s}}{m}_{\mathrm{s}}}}{2{\sigma }_{\mathrm{ss}}},\end{array}$$(A.6)

where the definition of the collision frequency ν_ss given by Eq. (8) has been used and σ_ss represents the cross-section of collisions between particles of the same species.

With the above definitions, we can compute the Reynolds numbers. Using the charged fluid as reference, the former dimensionless number is given by

$$\begin{array}{c}\hfill R_{\mathrm{e}}=\frac{\mathrm{\nabla }·\left({\rho }_{\mathrm{c}}{\mathit{V}}_{\mathrm{c}}{\mathit{V}}_{\mathrm{c}}\right)}{\mathrm{\nabla }·{\widehat{\tau }}_{\mathrm{c}}}\sim \frac{{\rho }_{\mathrm{c}}{V}_{\mathrm{c}}^2/L}{{\xi }_{\mathrm{c}}{V}_{\mathrm{c}}/L^2}=\frac{{\rho }_{\mathrm{c}}{V}_{\mathrm{c}}L}{{\xi }_{\mathrm{c}}},\end{array}$$(A.7)

where V_c and L are reference values of the velocity of the flow and the length-scale of the system, respectively. In addition, the density of the charged fluid can be expressed in terms of ρ_T = ρ_c + ρ_n, the total density of the plasma, as ρ_c = χ_cρ_T, where χ_c is the ionisation degree. Then, the Reynolds number can be rewritten as

$$\begin{array}{c}\hfill R_{\mathrm{e}}=\frac{{\chi }_{\mathrm{c}}{\rho }_{\mathrm{T}}{V}_{\mathrm{c}}L}{{\xi }_{\mathrm{c}}}.\end{array}$$(A.8)

Therefore, the Reynolds number varies as the ionisation degree of the plasma is varied.

For the purposes of this analysis, we can define a dimensionless number which plays an analogue role to Reynolds number. Here, we refer to that quantity as the numerical Reynolds number, which is given by

$$\begin{array}{c}\hfill R_{\mathrm{e},\mathrm{num}}=\frac{\mathrm{\nabla }·\left({\rho }_{\mathrm{c}}{\mathit{V}}_{\mathrm{c}}{\mathit{V}}_{\mathrm{c}}\right)}{\mathrm{\nabla }·{\widehat{\tau }}_{\mathrm{c},\mathrm{art}}}\sim \frac{V_{\mathrm{c}}L}{{\nu }_l}\end{array}$$(A.9)

and does not depend on the ionisation degree, in contrast with the real Reynolds number.

For our simulations, we decided to neglect the effect of the physical viscosity and we only take into account the numerical one. In this way, the effective Reynolds number is kept the same for every series of simulations and the comparisons between each other are more straightforward.

As described in Vögler et al. (2005) and Felipe et al. (2010), the strength of the numerical dissipation applied to each one of the temporal evolution equations presented in Sect. 2.1 is controlled by a different coefficient. However, in this work, we have chosen the same value for all the coefficients, given by a_diff.

In Fig. A.1 we present a comparison of the effect of numerical diffusivity for the cases of χ_c = 0.5 with α_eff = 0 (black lines) and α_eff = α (red lines), and N₀ = 1000. The spatially averaged magnetic field is plotted as a function of the normalised time. The results for the three values of a_diff coincide almost perfectly during the linear phase of the instability, while huge differences appear at later stages. For example, at the end of the simulation (τ = 15), the magnetic field for the case with a_diff = 2.5 × 10⁻⁴ is approximately one order of magnitude larger than for the case with the strongest diffusivity. The behaviour of the magnetic field displayed by this figure supports the statement that most of the dissipation occurs at the smallest length scales.

Fig. A.1. Mean magnetic field as a function of time and numerical diffusivity for the case with χc = 0.5. The solid line represents a simulation with adiff = 2.5 × 10−4 (in normalised units) while the dashed and the dotted lines correspond to simulations with adiff = 2.5 × 10−3 and adiff = 2.5 × 10−2, respectively. Black and red lines represent the cases with αeff = 0 and αeff = α.

The effect of a larger numerical diffusivity can also be checked in Fig. A.2, which displays density snapshots for the case with a_diff = 2.5 × 10⁻². The comparison of Figs. 2 and A.2 shows that in the former the edges of the vortexes are sharper while in the latter they are more blurred. Furthermore, the secondary vortexes present at τ = 7 in the case of low diffusivity do not appear when a much larger diffusivity is used. Finally, most of the small scale structure that can be seen in the right panels of Fig. 2 is removed in Fig. A.2.

Fig. A.2. Density snapshots of a simulation with χc = 0.5, αeff = α, and a high value of numerical diffusivity, adiff = 2.5 × 10−2. Same colour scheme as in Fig. 2 is used. An animation of this figure is available online.

Appendix B: Additional parameters of the plasma

Using the definition given by Eq. (A.6), the physical viscosities of the charged and neutral fluid are ξ_c ≈ 1.2 × 10⁻⁷ kg m⁻¹ s⁻¹ and ξ_n ≈ 1.7 × 10⁻⁵ kg m⁻¹ s⁻¹, respectively, where the following collisional cross sections have been used: σ_cc = σ_ii ≈ 1.2 × 10⁻¹⁶ m² and σ_nn = 7.73 × 10⁻¹⁹ m² (Meier & Shumlak 2012; Leake et al. 2013).

The thermal conductivities, κ_c and κ_n, are given by Braginskii (1965), Leake et al. (2013):

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{c}}=\frac{2k_{\mathrm{B}}}{{\sigma }_{\mathrm{cc}}}\sqrt{\frac{\pi k_{\mathrm{B}}{T}_{\mathrm{c}}}{m_{\mathrm{c}}}}\approx 3.8\times {10}^{-3}\mathrm{W}{\mathrm{m}}^{-1}{\mathrm{K}}^{-1},\end{array}$$(B.1)

and

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{n}}=\frac{2k_{\mathrm{B}}}{{\sigma }_{\mathrm{nn}}}\sqrt{\frac{\pi k_{\mathrm{B}}{T}_{\mathrm{n}}}{m_{\mathrm{n}}}}\approx 0.57\mathrm{W}{\mathrm{m}}^{-1}{\mathrm{K}}^{-1}.\end{array}$$(B.2)

Then, the Peclet numbers, which represent the ratio of advective transport rate to the diffusive transport rate due to thermal conduction, are given by

$$\begin{array}{c}\hfill P_{\mathrm{e},\mathrm{c}}=\frac{L_0{U}_0{\rho }_{\mathrm{c}}{c}_{\mathrm{p},\mathrm{c}}}{{\kappa }_{\mathrm{c}}}\approx 13600\end{array}$$(B.3)

and

$$\begin{array}{c}\hfill P_{\mathrm{e},\mathrm{n}}=\frac{L_0{U}_0{\rho }_{\mathrm{n}}{c}_{\mathrm{p},\mathrm{n}}}{{\kappa }_{\mathrm{c}}}\approx 100,\end{array}$$(B.4)

where c_c, s = γ/(γ − 1)k_B/m_s is the heat capacity of species s. The values of these dimensionless numbers reveal that the thermal conduction, especially that of the neutral fluid, may play a relevant role in the dynamics of the plasma. However, for the sake of simplicity and to focus only on the study of the charged-neutral collisional interaction, we do not take into account the effect of thermal conduction in the present research.

Appendix C: Comparison between Joule and collisional heating

Taking into account the definition for the variable E_diff, given by Eq. (27) and Ohm’s law, Eq. (20), we get that the Joule heating is given by

$$\begin{array}{c}\hfill \mathit{J}·{\mathit{E}}_{\mathrm{diff}}=\mathit{J}·[-{\eta }_{\mathrm{H}}\mathrm{\nabla }{P}_{\mathrm{e}}+{\eta }_{\mathrm{D}}({\mathit{V}}_{\mathrm{n}}-{\mathit{V}}_{\mathrm{c}})].\end{array}$$(C.1)

The results presented in the main text demonstrate that the collisional term of Ohm’s law can be neglected. So, the previous equation is reduced to

$$\begin{array}{c}\hfill \mathit{J}·{\mathit{E}}_{\mathrm{diff}}=-{\eta }_{\mathrm{H}}\mathit{J}·\mathrm{\nabla }{P}_{\mathrm{e}}=-\frac{1}{e{n}_{\mathrm{e}}}\frac{\mathrm{\nabla }\times \mathit{B}}{{\mathrm{\mu }}_0}·\mathrm{\nabla }{P}_{\mathrm{e}}.\end{array}$$(C.2)

Then, the ratio between the Joule heating term and the collisional term ($Q_{\mathrm{cn}}^T$) can be estimated as:

$$\begin{array}{c}\hfill \frac{|\mathit{J}·{\mathit{E}}_{\mathrm{diff}}|}{|Q_{\mathrm{cn}}^T|}\sim \frac{(\gamma -1)B_{\mathrm{ref}}}{{\left(\mathrm{\Delta }x\right)}^{2}e{\mathrm{\mu }}_0{n}_{\mathrm{c}}{\nu }_{\mathrm{cn}}}\sim {10}^{-4},\end{array}$$(C.3)

where the value of B_ref corresponds to the maximum magnetic field found in the simulations. This relation demonstrates that the term J ⋅ E_diff can be neglected in the present study.

All Tables

All Figures

	Fig. 1. Vertical cuts of the initial conditions of the simulations. Top: normalised density, Θ (red solid line) and normalised temperature (green dashed line) at x = 0. Middle: normalised flow speed at x = 0. Bottom: normalised perturbation of the z-component of the velocity at x = −L0/4 (solid line) and x = L0/4 (dashed line).
In the text

	Fig. 2. Density snapshots of a simulation with χc = 0.5 with no coupling between the charged (orange top half) and neutral (blue bottom half) fluids. Time is measured in the units of τ ≡ t/t0, where t0 ≡ L0/U0 is the domain crossing time. Asterisks represent key points whose temporal evolution is analysed in Sect. 3.1.3. An animation of this figure is available online.
In the text

	Fig. 3. Differences in the evolution of density of charged and neutral fluids for the case of χc = χn = 0.5; |Δρ|≡|ρc − ρn|. Top: αeff = 0 (no collisions); bottom: αeff = α (charge-exchange collisions not taken into account). The black lines represent density iso-contours.
In the text

	Fig. 4. Velocity drifts, ΔVx ≡ Vx,c − Vx,n, for the same simulations as in Fig. 3. The black lines represent density iso-contours. We note that the units of the top colourbar are km s−1 and the units of the bottom colourbar are m s−1: the presence of the collisional coupling reduces the magnitude of the velocity drifts.
In the text

	Fig. 5. Absolute value of the magnetic field, |By|, at points P1 (blue dotted line), P2 (red solid line), and P4 (black dashed line) as a function of the normalised time, τ.
In the text

	Fig. 6. Temperature as a function of time at points P1 (blue thin line), P2 (red thick line) and P4 (black line). Solid lines represent the temperature of the charged fluid while dashed lines represent the temperature of the neutral fluid. Left and right panels: cases with αeff = 0 and αeff = α.
In the text

Fig. 7. Averaged magnetic field as a function of time and resolution for the cases with χc = 0.5 (left) and χc = 0.1 (right). Solid lines represent simulations in which the ion-neutral collisional interaction is taken into account (without charge-exchange). Dotted lines represent simulations without collisional coupling. Green, blue, purple, and red colours correspond to simulations with N0 = 250, 500, 1000, 2000, respectively. The vertical dashed lines approximately mark the separation between the linear and non-linear phases of the instability. The inset in each panel shows a zoom of the interval τ ∈ [5, 15].

In the text

Fig. 8. Spatially averaged magnetic field as a function of ionisation degree at three different times of the simulations. Black asterisks show the results without collisions. Red diamonds correspond to simulations in which ion-neutral collisions are considered. Blue crosses represent the results when the collisional term of the induction equation is also taken into account and the green squares represent the simulations that also included the process of charge-exchange.

In the text

	Fig. 9. Spatially averaged values of the baroclinic terms as functions of the ionisation degree, χc, at three different times: τ = 3 (left), τ = 6 (middle), and τ = 9 (right). Solid and dashed lines correspond to the terms ⟨|Bωc|⟩ and ⟨|Bωn|⟩, respectively. The black, red and green colours correspond to the cases with αeff = 0, αeff = α, and αeff = α + αcx, respectively.
In the text

	Fig. 10. Spatially averaged values of the collisional terms in the temperature equations as functions of the ionisation degree, at three different times: τ = 3 (left), τ = 6 (middle), and τ = 9 (right). Solid and dashed lines represent the terms for the charged and the neutral fluids, respectively, which are given by Eqs. (63) and (64).
In the text

Fig. 11. Spatially averaged magnetic field as a function of ionisation degree for simulations with N0 = 500 and different values of total density, ρT. Symbols and lines correspond to simulations with αeff = 0 and αeff = α, respectively. The blue dotted, grey dashed-dotted, and red dashed lines represent the cases with ρT = 10−11,  10−10 and 10−9 kg m−3, respectively. The black line corresponds to the prediction of the single-fluid model.

In the text

	Fig. A.1. Mean magnetic field as a function of time and numerical diffusivity for the case with χc = 0.5. The solid line represents a simulation with adiff = 2.5 × 10−4 (in normalised units) while the dashed and the dotted lines correspond to simulations with adiff = 2.5 × 10−3 and adiff = 2.5 × 10−2, respectively. Black and red lines represent the cases with αeff = 0 and αeff = α.
In the text

	Fig. A.2. Density snapshots of a simulation with χc = 0.5, αeff = α, and a high value of numerical diffusivity, adiff = 2.5 × 10−2. Same colour scheme as in Fig. 2 is used. An animation of this figure is available online.
In the text