1. Introduction
The interior atmosphere of Jupiter-like gas giant planets is typically characterized by two dynamically distinct regions: a neutral outer envelope following the laws of hydrodynamics (HD), and a hot ionized interior where the hydrogen transitions into a conducting metallic liquid state that produces interactions between the momentum and magnetic field, following the laws of magnetohydrodynamics (MHD) (Guillot Reference Guillot2005; Liu, Goldreich & Stevenson Reference Liu, Goldreich and Stevenson2008; Stanley & Glatzmaier Reference Stanley and Glatzmaier2010). While the location of this transition might change depending on the planet's mass and age, the existence of the two regions is expected to be robust, given typical pressures, temperatures and compositions of gas giant planets. Attempts at modelling each region individually have successfully reproduced and helped with the understanding of a lot of the major observations of Jupiter and Saturn to this date. This includes, among many other things, the formation, characteristics and dynamics of the jets and vortices (Rhines Reference Rhines1975; Busse Reference Busse1976; Dowling & Ingersoll Reference Dowling and Ingersoll1988, Reference Dowling and Ingersoll1989; Cho & Polvani Reference Cho and Polvani1996; Showman Reference Showman2007; Scott & Polvani Reference Scott and Polvani2007, Reference Scott and Polvani2008; Liu et al. Reference Liu, Goldreich and Stevenson2008; Lian & Showman Reference Lian and Showman2008, Reference Lian and Showman2010; Glatzmaier, Evonuk & Rogers Reference Glatzmaier, Evonuk and Rogers2009; Schneider & Liu Reference Schneider and Liu2009; Stanley & Glatzmaier Reference Stanley and Glatzmaier2010; Warneford & Dellar Reference Warneford and Dellar2014; O'Neill, Emanuel & Flierl Reference O'Neill, Emanuel and Flierl2015; Heimpel, Gastine & Wicht Reference Heimpel, Gastine and Wicht2016) as well as the generation and morphology of the magnetic field (Stanley & Glatzmaier Reference Stanley and Glatzmaier2010; Wicht & Tilgner Reference Wicht and Tilgner2010; Jones et al. Reference Jones, Boronski, Brun, Glatzmaier, Gastine, Miesch and Wicht2011; Gastine et al. Reference Gastine, Wicht, Duarte, Heimpel and Becker2014; Jones Reference Jones2014; Rogers & McElwaine Reference Rogers and McElwaine2017; Dietrich & Jones Reference Dietrich and Jones2018; Duarte, Wicht & Gastine Reference Duarte, Wicht and Gastine2018).
However, in reality, the two regions are not completely independent of each other. The transition from neutral to fully ionized is believed to happen continuously as a function of the radius (Bagenal, Dowling & McKinnon Reference Bagenal, Dowling and McKinnon2006; Liu et al. Reference Liu, Goldreich and Stevenson2008; Cao & Stevenson Reference Cao and Stevenson2017; Zaghoo Reference Zaghoo2018). This suggests a continuous transition from HD to MHD dynamics. Many of the modelling studies mentioned above have either ignored or parametrized the interaction of their modelling region with this transition region, where the dynamics begins to change. Some examples include the presence of a Rayleigh-like friction in general circulation models (GCMs) modelling the neutral region. This friction is coined ‘MHD drag’ and is supposed to account for the interaction between the jets and the metallic interior (Glatzmaier Reference Glatzmaier2008; Liu et al. Reference Liu, Goldreich and Stevenson2008; Schneider & Liu Reference Schneider and Liu2009; Perna, Menou & Rauscher Reference Perna, Menou and Rauscher2010). On the other hand, some modelling efforts of the interior MHD region have begun including the effects of variable conductivity as a function of radius, accounting for the steep drop-off of conductivity as one approaches the neutral region (Jones et al. Reference Jones, Boronski, Brun, Glatzmaier, Gastine, Miesch and Wicht2011; Gastine et al. Reference Gastine, Wicht, Duarte, Heimpel and Becker2014; Jones Reference Jones2014; Dietrich & Jones Reference Dietrich and Jones2018). Despite the relative success of these methods, some call into question these techniques (Glatzmaier Reference Glatzmaier2008; Chai, Jansen & Vallis Reference Chai, Jansen and Vallis2016), and many studies from both communities explicitly express interest in further understanding the coupling between the HD and MHD regions (Gastine et al. Reference Gastine, Wicht, Duarte, Heimpel and Becker2014; Jones Reference Jones2014; Chai et al. Reference Chai, Jansen and Vallis2016; Heimpel et al. Reference Heimpel, Gastine and Wicht2016).
A better understanding of the transition region could also be important for understanding hot Jupiters, gas giant planets orbiting close to other stars. The outer atmospheres of hot Jupiters are expected to be ionized, due to both high temperatures and incident radiation from the nearby host star (Batygin & Stevenson Reference Batygin and Stevenson2010; Perna et al. Reference Perna, Menou and Rauscher2010; Menou Reference Menou2012; Koskinen et al. Reference Koskinen, Yelle, Lavvas and Cho2014; Koll & Komacek Reference Koll and Komacek2018), and a few studies have already implemented GCMs that include an electrically conducting atmosphere and a magnetic field (Batygin, Stanley & Stevenson Reference Batygin, Stanley and Stevenson2013; Rogers & Showman Reference Rogers and Showman2014; Rogers & McElwaine Reference Rogers and McElwaine2017). The resulting circulations depend significantly on the interaction of the flows with the magnetic field, which has implications for the interpretation of observed hot spots on hot Jupiters.
All the previous numerical work described above has been carried out using either regular HD, which models the dynamics of an electrically neutral fluid, or MHD, which models the dynamics of a fully ionized electrically conducting fluid, incorporating the interaction between the fluid and the magnetic field. These are called single-fluid or single-species models because they model only one type of molecule (either neutral or ionized). However, it is likely that this continuous transition occurs via partial ionization (Zaghoo Reference Zaghoo2018), implying a coexistence of ionized and neutral molecules in this region, both following their own respective mean dynamics but occupying the same fluid volume and interacting via collisions.
The main goal of this work is to improve our understanding of the partially ionized turbulent dynamics occurring in the transition region. A more rigorous understanding of the plasma physics and dynamical regimes there can shed light on the commonly used assumptions. To that effect, in § 2 we introduce and explore partially ionized magnetohydrodynamics (PIMHD), a two-fluid model – one neutral and one ionized – coupled by a collision term proportional to the difference in velocities. Unlike the single-species (fully ionized or fully neutral) models, the coexistence of two species introduces a new frictional dissipation of energy and a source of heating due to collisions between the differentially moving species. In § 3 we motivate and describe the approach to our study using numerical simulations of two-dimensional (2-D) turbulence, whose results are presented and discussed in § 4. Finally, in § 5 we summarize the PIMHD numerical experiments and give tentative parameter values for Jupiter, making some connections to our motivating discussion in the current section.
2. Partially ionized magnetohydrodynamics
2.1. PIMHD system
In this study, we investigate incompressible PIMHD with uniform species densities, ignoring the complications of compression, stratification and buoyancy. Incompressibility is expected to hold true in the deep atmospheres of gas giant planets, although this is possibly less accurate for the outer regions of hot Jupiter atmospheres where partial ionization is also relevant. The assumption of uniform species densities is harder to justify, and we admit that it would certainly play a role in real geophysical applications at large scales. Despite this, we make these assumptions because they simplify analysis and allow better comparison to previously established turbulence results. The PIMHD system then becomes the following:
\begin{gather} \rho_n \left(\frac{\partial}{\partial t} + \boldsymbol{v}_n \boldsymbol{\cdot} \boldsymbol{\nabla}\right) \boldsymbol{v}_n ={-} \boldsymbol{\nabla} p_n -\rho_i \rho_n \alpha ( \boldsymbol{v}_n- \boldsymbol{v}_i) + \mu_n \boldsymbol{\nabla}^2 \boldsymbol{v}_n + \boldsymbol{F}_n, \end{gather}
\begin{gather} \rho_i \left(\frac{\partial}{\partial t} + \boldsymbol{v}_i \boldsymbol{\cdot} \boldsymbol{\nabla}\right) \boldsymbol{v}_i ={-} \boldsymbol{\nabla} p_i + \rho_i \rho_n \alpha ( \boldsymbol{v}_n- \boldsymbol{v}_i) + \boldsymbol{J} \times \boldsymbol{B} + \mu_i \boldsymbol{\nabla}^2 \boldsymbol{v}_i + \boldsymbol{F}_i, \end{gather}
\begin{gather} \frac{\partial \boldsymbol{B}}{\partial t} = \boldsymbol{\nabla} \times (\boldsymbol{v}_i \times \boldsymbol{B}) + \eta \boldsymbol{\nabla}^2 \boldsymbol{B} + \boldsymbol{F}_B, \end{gather}
\begin{gather} \boldsymbol{J} = \frac{1}{\mu_0} \boldsymbol{\nabla}\times\boldsymbol{B}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{B} = 0, \quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}_n = 0, \quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}_i = 0, \quad \rho_{tot} = \rho_i + \rho_n, \end{gather}
where the subscripts represent an ionized (‘i’) or neutral (‘n’) component, potentially representing dissociated hydrogen ions and recombined atoms, respectively (Guillot Reference Guillot2005; French et al. Reference French, Becker, Lorenzen, Nettelmann, Bethkenhagen, Wicht and Redmer2012). For each species, 
$\boldsymbol {v}$ is the velocity, 
$\rho$ is the density, 
$p$ is the pressure, 
$\mu$ is the dynamic viscosity, 
$\boldsymbol {B}$ is the magnetic field, 
$\mu _0$ is the vacuum permeability, 
$\eta$ is the magnetic diffusivity, and 
$\boldsymbol {F}$ is a generic force that may include body or gravitational forces and other forms of dissipation, for example. These two species do not form two different layers – they occupy the same space, i.e. at each point there are two velocities corresponding to 
$\boldsymbol {v}_i$ and 
$\boldsymbol {v}_n$.
The PIMHD system can be derived from the Boltzmann equations for singly charged ions, electrons and a single neutral species (Draine Reference Draine1986; Meier Reference Meier2011; Meier & Shumlak Reference Meier and Shumlak2012). The derivation process is similar to that of the MHD system from an electron–ion plasma, the difference being the presence of extra collision and reaction terms between neutrals, ions and electrons. One combines the momentum equations of each ionized species with Maxwell's equations, ignoring any static charge sources (quasi-neutral approximation) and light waves (the electric field is set by Ohm's law). Further simplifications in the dynamics mainly come from ignoring the electron inertia and electron pressure, thus reducing the equations of motion to that of the magnetic field and of the centre of mass between the ions and electrons. Making these approximations in MHD requires assuming that the electron-to-ion mass ratio is very small, as well as assuming that we are looking at length scales much larger than the ion or electron skin depths. We want to emphasize here that we are taking the MHD approximation for the ion species, which ignores many two-fluid effects commonly considered in astrophysical plasmas (Ballester et al. Reference Ballester, Alexeev, Collados, Downes, Pfaff, Gilbert, Khodachenko, Khomenko, Shaikhislamov and Soler2018). We believe that the MHD regime is valid in the system we are attempting to study, and we discuss the breakdown of some of these assumptions in § 2.2 and appendix A.
Partially ionized models have been used in previous studies of the Earth's thermosphere/ionosphere, the Sun's chromosphere (Khodachenko et al. Reference Khodachenko, Arber, Rucker and Hanslmeier2004; Zaqarashvili, Khodachenko & Rucker Reference Zaqarashvili, Khodachenko and Rucker2011; Khomenko & Collados Reference Khomenko and Collados2012; Leake et al. Reference Leake, DeVore, Thayer, Burns, Crowley, Gilbert, Huba, Krall, Linton and Lukin2014; Martínez-Sykora et al. Reference Martínez-Sykora, De Pontieu, Hansteen, Rouppe van der Voort, Carlsson and Pereira2017; Song Reference Song2017), magnetic reconnection (Lazarian, Vishniac & Cho Reference Lazarian, Vishniac and Cho2004; Smith & Sakai Reference Smith and Sakai2008; Malyshkin & Zweibel Reference Malyshkin and Zweibel2011; Leake et al. Reference Leake, Lukin, Linton and Meier2012), protoplanetary disks (Balbus Reference Balbus2009), as well as molecular clouds and the interstellar medium (Draine Reference Draine1980; Nakano & Umebayashi Reference Nakano and Umebayashi1986; Falle Reference Falle2003; Oishi & Mac Low Reference Oishi and Mac Low2006; O'Sullivan & Downes Reference O'Sullivan and Downes2007; Tilley & Balsara Reference Tilley and Balsara2010; Meyer et al. Reference Meyer, Balsara, Burkhart and Lazarian2014; Xu, Yan & Lazarian Reference Xu, Yan and Lazarian2016; Xu & Lazarian Reference Xu and Lazarian2017). As far as we are aware, there have not been any studies looking at partially ionized turbulence in the planetary atmosphere setting we are investigating here. Some of the main differences between the context of previous work and our deep atmosphere context is that the former typically deals with ionization fractions much smaller than one, as well as compressibility effects. We have also not found any study that does a systematic parameter space study of the turbulent PIMHD system that will be discussed in § 3.
The two fluids are coupled via collisions, represented by the second term on the right-hand sides of (2.1a) and (2.1b). The coefficient 
$\alpha$ measures the strength of the coupling; it is approximately proportional to the collision cross-section of the two species and their thermal velocities, which in turn depends on the square-root of the temperature under equilibrium assumptions (Draine Reference Draine1986; Meier Reference Meier2011; Leake, Lukin & Linton Reference Leake, Lukin and Linton2013). For our purposes, 
$\alpha$ will be a parameter that we vary, although we will discuss possible realistic values of 
$\alpha$ in § 5. Looking at the energy equation (and ignoring other forms of dissipation) we see the effects of collisions on the energy of the individual species (
$s \in \{i,n\}$) and as a whole:
\begin{gather} \frac{\textrm{d} E_s}{\textrm{d}t} ={-} \rho_i \rho_n \alpha \langle | \boldsymbol{v}_s|^2\rangle + \rho_i \rho_n \alpha \langle \boldsymbol{v}_n \boldsymbol{\cdot} \boldsymbol{v}_i \rangle, \end{gather}
\begin{gather} \frac{\textrm{d} E}{\textrm{d}t} ={-} \rho_i \rho_n \alpha \langle | \boldsymbol{v}_i - \boldsymbol{v}_n|^2\rangle, \end{gather}
where 
$\langle \cdot \rangle$ implies a domain integral, 
$E_n = KE_n = \rho _n \langle |\boldsymbol {v}_n|^2 \rangle /2$, 
$E_i = KE_i + E_B = \rho _i \langle |\boldsymbol {v}_i|^2 \rangle /2 + \langle |\boldsymbol {B}|^2 \rangle / (2 \mu _0)$ and 
$E = E_i + E_n$. Collisions conserve momentum, and the sign-indefinite term in (2.2a) tells us that the two species may exchange some energy via the collisions. Looking at (2.2b) we see that total energy is lost from these interactions, in a process we are calling ‘collisional heating’ (CH) (Vasyliunas & Song Reference Vasyliunas and Song2005). In the absence of collisions and other dissipation terms, the two species are uncoupled and behave as HD and MHD independently. Note that CH is something not accounted for in one-fluid models, and could therefore prove problematic if it is shown to be significant in these planetary systems.
Another new and important parameter of the PIMHD system is the ionization fraction, which we will denote by 
$\chi \equiv \rho _i/\rho _{tot}$. Since 
$\rho _{tot} = \rho _i + \rho _n$, we see that 
$(1-\chi ) = \rho _n/\rho _{tot}$. The ionization fraction plays a role in the dynamics by influencing the acceleration and Reynolds number (which measures the relative strength of the advection term to the diffusion term) of each species, as well as the strength of the collision term. We note here that (2.1a)–(2.1d) are not valid for ionization fractions that are strictly 0 or 1, as the fluid description breaks down as those limits are approached. In the very extreme limits, we expect that certain assumptions made about time-scale separations between molecular motion and mean motion are no longer valid as the densities become low and the mean free path for self-collisions becomes too large (e.g. leading to the breakdown of the thermal equilibrium assumption) (Draine Reference Draine1986). Furthermore, before this occurs, as one approaches 
$\chi \rightarrow 0$, Ohm's law must be altered, as will be discussed in § 2.2 and appendix A. In any case, since we expect transition regions in gas giant planets to span all values of ionization fraction from 0 to 1, we will not focus on extremes of ionization fraction in this work.
Ionization fraction also modifies the magnetic diffusivity 
$\eta$. In MHD the origin of magnetic diffusivity comes from collisions between ions and electrons. However, in a partially ionized system, we also have collisions between neutrals and electrons. Incorporating this in the expression for magnetic diffusivity gives us
\begin{equation} \eta = \left(1+r\left(\frac{1-\chi}{\chi}\right)\right)\eta_{MHD}, \end{equation}
where 
$r$ is the ratio of cross-sections of ion–electron collisions to neutral–electron collisions, and 
$\eta _{MHD}$ is the magnetic diffusivity for regular MHD. Typically we would expect 
$r \ll 1$ (Leake et al. Reference Leake, Lukin, Linton and Meier2012), implying a sudden increase in magnetic diffusivity for small values of 
$\chi$ (low ionization fraction). Since we will not be dealing with extreme values of 
$\chi$ in this work, this effect will not be important here.
2.2. Limiting cases
Before moving on to the numerical experiments in § 3, we find it useful to explore the limiting cases of the PIMHD system (while staying within the bounds of our assumptions that make the system valid, as discussed in the previous section). We will look at the extreme limits of 
$\alpha$ and comment briefly on the ionization fraction limits.
The relative strength of the collision and advection terms in (2.1a) and (2.1b) is determined by the ratio of the eddy turnover time to the time scale of collisions. We define the eddy turnover time in the usual way as 
$\tau _{eddy} \equiv L/U$, where 
$L$ is a typical length scale of the system, and 
$U$ is a typical velocity. There are two time scales for collisions: 
$\tau _{coll,i} \equiv (\rho _n \alpha )^{-1}$ in the ion equation and 
$\tau _{coll,n} \equiv (\rho _i \alpha )^{-1}$ in the neutral equation. These are typically denoted in the literature as collision frequencies 
$\nu _{in} = \rho _n \alpha$ and 
$\nu _{ni} = \rho _i \alpha$. Since at the moment we are dealing with both 
$\chi \sim \textit{O}(1)$ and 
$(1-\chi ) \sim \textit{O}(1)$, it is convenient to define 
$\tau _{coll} \equiv (\rho _{tot}\alpha )^{-1}$ so that 
$\tau _{coll,i} = \tau _{coll}/(1-\chi )$ and 
$\tau _{coll,n} = \tau _{coll}/\chi$. This allows us to define our second main non-dimensional parameter (after 
$\chi$),
\begin{equation} \tilde{\alpha} \equiv \frac{\tau_{eddy}}{\tau_{coll}} = \frac{L \rho_{tot} \alpha}{U}, \end{equation}
which will determine the strength of the collision term compared to the advection term in the PIMHD system and will be a measure of how coupled the two fluids are. The limits in the cases below are really being applied to 
$\tilde {\alpha }$.
Part of our goal is to predict the collisional heating, defined to be
\begin{equation} CH \equiv \rho_i \rho_n \alpha \langle | \boldsymbol{v}_i - \boldsymbol{v}_n|^2\rangle. \end{equation}
The collisional heating is not only of interest for its implications to the astrophysical systems mentioned in § 1, but also because it is a measure of how coupled the two fluids are. Equation (2.5) and its equivalent wavenumber spectrum (to be defined) will be of interest throughout the rest of this work. We expect two extreme regimes: (1) 
$\tilde {\alpha } \ll 1$, where the ions and neutrals are not coupled and therefore follow their own separate dynamics, colliding and exchanging energy as they do so; and (2) 
$\tilde {\alpha } \gg 1$, where the ions and neutrals are extremely coupled, meaning 
$\boldsymbol {v}_i \approx \boldsymbol {v}_n$, thereby lowering 
$CH$. We will look at the two regimes separately and discuss some findings in each.
Let us begin with the limit of 
$\tilde {\alpha } \ll 1$. In the case where 
$\tilde {\alpha } = 0$ exactly, we recover two uncoupled fluids behaving as HD and MHD with no collisional heating. But suppose now that 
$0 < \tilde {\alpha } \ll 1$. In the case of isotropic three-dimensional (3-D) turbulence, since both HD and MHD turbulence cascade energy to smaller scales (Alexakis & Biferale Reference Alexakis and Biferale2018), we expect that the collisional heating will become negligible once 
$\tilde {\alpha } < 1$ and will eventually go to zero as we keep decreasing 
$\tilde {\alpha }$. The case of 2-D PIMHD turbulence is quite different due to the presence of an inverse cascade of energy in the neutral species, which causes energy to go to larger and larger scales. The energy builds until some dissipative force is able to balance it. For a finite-size domain of typical length 
$L_0$, as long as 
$\alpha > \mu _n/(\rho _i \rho _n L_0^2)$, this dissipative force is not the viscosity but the collisional heating. At steady state we expect collisions to balance the energy injected into the neutrals, which we call 
$I_n \equiv \langle \boldsymbol {v}_n \boldsymbol{\cdot} F_n \rangle$, and so 
$CH \approx I_n$ for 
$\tilde {\alpha } \ll 1$. This prediction becomes independent of 
$\alpha$ because, at steady state, all of the energy being injected at the forcing scale is expected to be dissipated away by collisional heating, and thus what changes for different values of 
$\alpha$ is not the dissipation rate, but the energy at the largest scales. Indeed, if we further assume that 
$|\boldsymbol {v}_n|\gg |\boldsymbol {v}_i|$, which should be the case for small values of 
$\tilde {\alpha }$ due to the inverse cascade of the neutrals but not the ions, then combining this result with (2.5) and the definition of 
$E_n$ we can say further that
\begin{equation} E \approx E_n \approx \frac{I_n}{2 \rho_i \alpha} \end{equation}
for 
$\tilde {\alpha } \ll 1$. The balance between 
$CH$ and 
$I_n$ has allowed us to approximately relate the energy of the neutral species with the energy injection rate, the ionization fraction and the collision coefficient.
If this limit of 
$\tilde {\alpha }$ were realized in the transition regions of gas giant planets, this could have possible implications for the saturation of the jets, whose formation is arguably attributed to the inverse cascade of kinetic energy in the presence of latitudinally varying rotation (Rhines Reference Rhines1975). Their saturation speeds, and therefore the effective Rhines scale, could depend on the value of 
$\alpha$. We should note here that we have assumed that the Rossby deformation radius is much larger than the domain size, but we do not expect our results to change for finite Rossby deformation radius since the arguments above still hold. Namely, the energy will still be dissipated at large scales (although possibly not the largest available scales) where viscous dissipation is negligible. Apart from a possible large-scale friction, (2.2a) tells us that collisions might be responsible for energy exchange between neutrals and ions. Indeed, when 
$|\boldsymbol {v}_n|\gg |\boldsymbol {v}_i|$, we might expect the second term on the right-hand side to dominate the friction-like term for the ion kinetic energy equation, thus leading to an injection of kinetic energy from the neutrals into the ions. Given that the ions would then cascade energy to the small scales, this could prove to be another route for the energy to be taken away from large scales and dissipated efficiently.
In the high collisional limit, 
$\tilde {\alpha } \gg 1$, we are also able to make some predictions. The following results are valid in both two and three dimensions, as they do not depend on any turbulent cascade. It is possible to do an asymptotic expansion of our variables in 
$\tilde {\alpha }^{-1}$, since this will be very small. Doing so leads us to conclude that, at 
$\textit{O}(\tilde {\alpha })$, 
$\boldsymbol {v}_i^{(0)} = \boldsymbol {v}_n^{(0)} \equiv \boldsymbol {u}$. Thus, to lowest order, the two fluids are completely coupled and 
$CH = 0$. Going to 
$\textit{O}(1)$ gives us, for the momentum equation of each species,
\begin{gather} \rho_{n} \left( \frac{\partial}{\partial t} + \boldsymbol{v}_n^{(0)} \boldsymbol{\cdot} \boldsymbol{\nabla}\right) \boldsymbol{v}_n^{(0)} ={-} \boldsymbol{\nabla} p_n - \rho_i \rho_n (\boldsymbol{v}_n^{(1)}- \boldsymbol{v}_i^{(1)}) + \mu_n \boldsymbol{\nabla}^2 \boldsymbol{v}_n^{(0)} + \boldsymbol{F}_n, \end{gather}
\begin{gather} \rho_{i} \left( \frac{\partial}{\partial t} + \boldsymbol{v}_i^{(0)} \boldsymbol{\cdot} \boldsymbol{\nabla}\right) \boldsymbol{v}_i^{(0)} ={-} \boldsymbol{\nabla} p_i + \rho_i \rho_n (\boldsymbol{v}_n^{(1)}- \boldsymbol{v}_i^{(1)}) + \boldsymbol{J} \times \boldsymbol{B} + \mu_i \boldsymbol{\nabla}^2 \boldsymbol{v}_i^{(0)} + \boldsymbol{F}_i, \end{gather}
The dynamics for 
$\boldsymbol {u}$ can be found by adding (2.7a) and (2.7b) and plugging in the fact that 
$\boldsymbol {v}_i^{(0)} = \boldsymbol {v}_n^{(0)} = \boldsymbol {u}$.
If we instead divide each equation by their respective density and subtract one from the other, we get rid of the left-hand sides and end up with an equation for 
$\boldsymbol {v}_n^{(1)}-\boldsymbol {v}_i^{(1)}$, which we can then use to get an expression for the next-order correction of the collisional heating 
$CH$. We end up with the following equations, which are correct down to 
$\textit{O}(1)$ in 
$\tilde {\alpha }^{-1}$:
\begin{gather} \rho_{tot} \left( \frac{\partial}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}\right) \boldsymbol{u} ={-} \boldsymbol{\nabla} (p_n+p_i) + \boldsymbol{J} \times \boldsymbol{B} + (\mu_n+\mu_i) \boldsymbol{\nabla}^2 \boldsymbol{u} + \boldsymbol{F}_n + \boldsymbol{F}_i, \end{gather}
\begin{gather} \frac{\partial \boldsymbol{B}}{\partial t} = \boldsymbol{\nabla} \times (\boldsymbol{u} \times \boldsymbol{B}) + \eta \boldsymbol{\nabla}^2 \boldsymbol{B} + \boldsymbol{F}_B, \end{gather}
\begin{gather} CH = \frac{\rho_n \rho_i}{\rho_{tot}^2} \frac{1}{\alpha} \left\langle \left|\frac{\boldsymbol{J}\times\boldsymbol{B}}{\rho_i} - \frac{\boldsymbol{\nabla} p_i}{\rho_i} + \frac{\boldsymbol{\nabla} p_n}{\rho_n} + \left(\frac{\mu_i}{\rho_i} - \frac{\mu_n}{\rho_n}\right)\nabla^2 \boldsymbol{u} + \frac{\boldsymbol{F}_i}{\rho_i} - \frac{\boldsymbol{F}_n}{\rho_n}\right|^2 \right\rangle. \end{gather}
Note that 
$CH \propto \alpha ^{-1}$, rather than the order-one correction one might expect, because the cross-terms in 
$|\boldsymbol {v}_n-\boldsymbol {v}_i|^2$ go away, leading to 
$|\boldsymbol {v}_n-\boldsymbol {v}_i|^2 = \alpha ^{-2}|\boldsymbol {v}_n^{(1)}-\boldsymbol {v}_i^{(1)}|^2$, which one combines with 
$CH \propto \alpha |\boldsymbol {v}_n-\boldsymbol {v}_i|^2$ to give an 
$\alpha ^{-1}$ dependence.
Looking at the dynamical equation for 
$\boldsymbol {u}$ reveals that it behaves like an MHD fluid, but with total densities, pressures, viscosities and forces. This suggests that, in the large-
$\tilde {\alpha }$ limit, the two fluids are coupled so that one-fluid models for the partially ionized region become valid. Since a turbulent fluid has a continuum of time scales, it is more appropriate to think of 
$\tilde {\alpha }$ as the collisional strength at a typical scale 
$L$ (which has a corresponding typical eddy turnover time and velocity), implying that the highly coupled limit could potentially only be valid up to a certain scale, depending on the actual value of 
$\tilde {\alpha }$ and choice of 
$L$. We will investigate the scale-by-scale properties of PIMHD in § 4.2. Equation (2.8c) gives us a prediction of the collisional heating based on order-one quantities. It tells us that collisions are caused by an imbalance in acceleration between the two species. For example, the ions feel the Lorentz force while the neutrals do not; therefore the magnetic field accelerates the ions in a different direction than the neutrals, which would then cause collisions. The results presented above, for both extremes of 
$\tilde {\alpha }$, will be tested and discussed in § 4.
Now we will briefly mention the limiting cases in ionization fraction, 
$\chi$, while maintaining the assumption of large 
$\tilde {\alpha }$. Owing to the separation of time scales that occurs when 
$\chi \ll 1$ or 
$(1-\chi )\ll 1$, it is numerically challenging to simulate these parameter limits, and we did not explore these regimes in our work. We therefore leave the details of the full discussion to appendix A and summarize the results here. In the fully ionized limit, the ions do not feel the collisions and make up most of the fluid, thus making the dominant dynamics single-fluid MHD, albeit with a modified pressure and body force in the highly collisional limit due to the fact that ions drag around neutrals. The low ionization limit is a bit more subtle. Certain assumptions in the derivation of MHD no longer hold, and thus the induction equation (2.1c) must be modified to include the Hall term (Pandey & Wardle Reference Pandey and Wardle2008), even at large scales. Furthermore, the neutrals, which dominate in this limit, still interact with the magnetic field indirectly via collisions with the ions, leading to what is called ‘ambipolar MHD’ in the high collisional limit, wherein the collisional effects act to enhance the magnetic diffusion. The limit of both low ionization and large collisional coupling, where ambipolar MHD is valid, are particularly relevant for many astrophysical applications, such as protoplanetary disks, molecular clouds and the interstellar medium (Draine Reference Draine1980; Oishi & Mac Low Reference Oishi and Mac Low2006; Balbus Reference Balbus2009; Tilley & Balsara Reference Tilley and Balsara2010; Meyer et al. Reference Meyer, Balsara, Burkhart and Lazarian2014; Xu et al. Reference Xu, Yan and Lazarian2016; Xu & Lazarian Reference Xu and Lazarian2017).
3. Methodology
Based on this introduction and discussion of the PIMHD system, we will now describe the numerical experiments used to test some of the predictions from § 2, as well as study the fully turbulent system scale by scale.
Our numerical experiments will comprise solely 2-D incompressible PIMHD turbulence. We acknowledge that a series of rotating 3-D simulations would be ideal. However, a large parameter sweep consisting of approximately 100 simulations at various ionization fractions, collision strengths and Reynolds numbers would be computationally demanding. We choose instead to explore this parameter space for the 2-D case first, with the expectation that it will provide guidance for future 3-D studies. There are two main reasons why we believe our choice is not restrictive and does not make this study irrelevant to its more realistic counterpart. The first relies on the fact that the high collisional results in § 2.2 are not dimension-dependent and thus should hold for both 2-D and 3-D turbulence. Secondly, those results which do depend on the dimensionality really only depend on the directions of the energy cascades, which we are respecting in our 2-D simulations, since 3-D rotating HD turbulence is expected to cascade energy to larger scales like 2-D HD turbulence, and 3-D MHD turbulence also cascades total energy to smaller scales.
Other simplifications will also be made for tractability of both the analysis and the numerics. In a realistic setting, assuming 
$\mu _s$ is not changing, the kinematic viscosity 
$\nu _s \equiv \mu _s/\rho _s$ will be a function of the ionization fraction and will thus affect the Reynolds number for each species. However, in this study we aim to isolate the effects of ionization fraction on the dynamics, and also wish to perform a large number of simulations. We therefore choose to keep the kinematic viscosity, and thus the Reynolds number, constant (and equal) for each species. For similar reasons we will fix 
$\eta$ so that the magnetic Prandtl number 
$Pr_m \equiv \nu _i/\eta = 1$ for all simulations. Although this is not likely to be true in realistic planetary settings (especially in the transition region where we expect the magnetic diffusivity to be large), we make this choice because the focus of this work is not to study the effects of 
$Pr_m$, which is purely an MHD parameter. This means we are ignoring the effects of the magnetic diffusivity's dependence on ionization fraction, seen in (2.3).
Equations (2.1a)–(2.1d) were solved in a doubly periodic domain with side length 
$2 {\rm \pi}$ using a modification of a 2-D MHD code written by Professor Pablo Mininni at the University of Buenos Aires, Argentina. The code was extended to include a neutral species and thus solve the PIMHD equations. It is a standard parallel pseudo-spectral code with a fourth-order Runge–Kutta scheme for time integration and a two-thirds dealiasing rule. More details on the parallelization can be found in Gómez, Mininni & Dmitruk (Reference Gómez, Mininni and Dmitruk2005). All runs started from random initial conditions, were continuously forced, and were carried out for long enough so that a statistically steady state was reached. All data were averaged at this state unless otherwise stated. A 
$512^2$ resolution was used for most of the experiments to explore the parameter space, with some 
$1024^2$ runs to ensure that the results are not resolution-dependent and to explore higher Reynolds numbers.
In an attempt to approach more realistic forcing mechanisms in geophysical flows, where convection or baroclinic instability might convert potential energy to kinetic energy, the forcing of each species is proportional to its density: 
$\boldsymbol {F}_s = \rho _s \boldsymbol {f}$, where 
$s \in \{i,n\}$. The forcing function 
$\boldsymbol {f}$ is identical for both species; and 
$\boldsymbol {f}$ is random, white-in-time and spectrally focused around wavenumber magnitude 
$k_f$, an input parameter. More specifically, at each time step, a wavenumber 
$\boldsymbol {k}_r$ of magnitude 
$k_f$ is chosen at random, and 
$\boldsymbol {\hat {f}}(\boldsymbol {k})$ (Fourier transform of 
$\boldsymbol {f}$) is set to zero everywhere except for at 
$\boldsymbol {k}_r$ where it has the magnitude 
$f_k/\sqrt {{\rm \Delta} t}$, with 
$f_k$ being another input parameter. This has the effect of setting the energy injection rate of each species to be 
$I_i = \rho _i f_k^2$ for ions and 
$I_n = \rho _n f_k^2$ for neutrals; see Chan, Mitra & Brandenburg (Reference Chan, Mitra and Brandenburg2012) for more details. Constant energy injection rate into the system is ideal for studying situations in which there is little to no large-scale dissipation and hence a large-scale condensate forms (Gallet & Young Reference Gallet and Young2013). The magnetic field was also forced using the function 
$\boldsymbol {f}$, but with a different random seed, and the magnetic energy injection rate was set to be 
$I_B = I_i/4$ for all runs. For most runs we set 
$k_f = 8$, so that both inverse and forward cascades could be resolved. However, 
$k_f = 4$ and 
$k_f = 32$ runs were carried out as well, to better resolve the forward and inverse cascade, respectively. In all runs, given different values of 
$k_f$, 
$I_s$ was chosen so that 
$u_f \equiv |\boldsymbol {v}_i(k_f)| \sim |\boldsymbol {v}_n(k_f)| \sim 1$ for both species.
In an attempt to better resolve inertial ranges while forcing at intermediate and larger wavenumbers, hyperviscosity, 
$(-1)^{p+1} \boldsymbol{\nabla}^{2p}$, was used in all runs and for all three fields, replacing the regular viscosity and regular magnetic diffusivity seen in (2.1a)–(2.1c). As long as the value of 
$p$ is not very large, hyperviscosity has been shown to have no significant effect on the turbulent properties of 3-D turbulence, and we expect the same to be the case for our work (Agrawal et al. Reference Agrawal, Alexakis, Brachet and Tuckerman2020). The value of 
$p$ was set to 2 in all runs except for those where 
$k_f = 32$, in which case 
$p = 4$. Furthermore, due to the inverse cascade of the square of the magnetic vector potential 
$|A|^2$ (Alexakis & Biferale Reference Alexakis and Biferale2018), where 
$\boldsymbol {\nabla } \times (A\hat {\boldsymbol {z}}) = \boldsymbol {B}$, we chose to include hypoviscosity in (2.1c) by adding 
$\eta _-\boldsymbol {\nabla} ^{-2} \boldsymbol {B}$. This acts to dissipate magnetic energy only at the largest scales and thus avoids the slow-forming condensate of 
$|A|^2$, making our simulations reach steady state faster. A condensate of the magnetic vector potential could possibly affect the dynamics of the ion species, but we consider this to be a purely MHD effect and thus not a focus of our work. The coefficient 
$\eta _-$ was chosen to ensure that the magnetic energy at the largest scales (
$|\boldsymbol {k}|=1$) was smaller than that at the next largest scale, thus avoiding the formation of a condensate.
In our simulations, we divided the momentum equations by 
$\rho _{tot}$ and absorbed it into the definition of our variables. Thus, in our simulations, 
$\boldsymbol {B} = \boldsymbol {B}/\sqrt {\mu _0 \rho _{tot}}$ and 
$\alpha = \rho _{tot} \alpha$ (making it equivalent to collision frequency), the latter being another input parameter. Doing this lets us directly employ the ionization fraction 
$\chi$ in the numerical integration of the equations. This means that 
$I_n = (1-\chi )f_k^2$ and 
$I_i = \chi f_k^2$. Using the four-fifths law, the typical velocity based on input parameters is 
$u_f = (\,f_k^2/k_f)^{1/3}$, which was maintained at 1 for all runs. Therefore, the eddy turnover time 
$\tau _{eddy} = (k_f f_k)^{-2/3}$. Furthermore, we define the numerical version of 
$\tilde {\alpha }$ based on input parameters:
\begin{equation} \tilde{\alpha} \equiv \frac{\alpha}{k_f u_f} = \frac{\alpha}{(k_f f_k)^{2/3}}. \end{equation}We can also define a parameter analogous to the Reynolds number in terms of our simulation input parameters:
\begin{equation} Re = \frac{f_k^2}{\nu k_f^{{-}2p+2/3}}, \end{equation}
where 
$p$ is the power of the hyperviscosity. Since we are keeping 
$\nu$ the same for both species, we do not distinguish between 
$Re_i$ or 
$Re_n$ and simply call it 
$Re$. Note that our use of hyperviscosity means that the parameter we call 
$Re$ is not exactly a Reynolds number. However, 
$Re$ can still be seen as the ratio of eddy to diffusive time scales, and thus measures the relative importance of advection compared to diffusion.
From now on, during any discussion of our numerical results, all references to these variables will be the numerical versions defined above. A summary of our runs can be seen in table 1.
Table 1. A summary of the runs performed for this work: 
$N$ is the resolution of the simulation and 
$p$ is the hyperviscosity exponent.
Since we are numerically integrating the two-fluid equations, when 
$\chi$ becomes very small we begin to encounter time-scale separation issues in the equations of motion, causing numerical difficulties (Falle Reference Falle2003; O'Sullivan & Downes Reference O'Sullivan and Downes2007). Furthermore, at extremely small values of 
$\chi$ we would be forced to alter the equations of motion to those seen in (A 1a) and (A 1b). For these reasons, and the fact that we are interested in spanning ionization fractions from zero to one, we chose to ignore extreme values of ionization fraction in our work. Set K8 was the most extensive set, with runs spanning 
$\tilde {\alpha } = [1.6\times 10^{-4}\ \textrm {to}\ 1.6\times 10^3]$. Since each set has a specific, fixed 
$k_f$ and 
$f_k$, 
$\tilde {\alpha }$ was modified by varying the numerical value of 
$\alpha$. For each value of 
$\tilde {\alpha }$, six runs were performed in which we varied 
$\chi$ from 0.1 to 0.99. The other sets were all run at 
$\chi = 0.5$ and focused mainly on the effect of the collisional strength on the dynamics.
In the forthcoming section, we will present and analyse the results of our simulations. We divide the analysis into two subsections: § 4.1, ‘global’ analysis, where we investigate the behaviour of volume-averaged statistics and compare to the predictions in § 2.2; and § 4.2, ‘spatial’ analysis, where we investigate the scale-by-scale effects of the collision strength on the two fluids.
4. Results
4.1. Global
In § 2.2 we saw that the collisional strength 
$\tilde {\alpha }$ sets the degree of coupling between the two fluids: in the small-
$\tilde {\alpha }$ limit we expect the two fluids to move independently, whereas, in the large-
$\tilde {\alpha }$ limit, they should be almost identical, following the MHD equations (2.8a) and (2.8b). In figure 1 we see three snapshots of the vorticity for each species, 
$\omega _s \equiv \boldsymbol {\nabla } \times \boldsymbol {v}_s$, representing, from left to right, 
$\tilde {\alpha } \ll 1$, 
$\tilde {\alpha } \sim 1$ and 
$\tilde {\alpha } \gg 1$, all from the K4 runs. The top row shows the neutral vorticity for each run, whereas the bottom row shows the ion vorticity for those runs. By visual inspection of figure 1(a) and (d), with 
$\tilde {\alpha } = 0.003$, one can indeed see that for small 
$\tilde {\alpha }$ the two species behave as they would in a completely uncoupled regime. The neutral vorticity shows clear signs of the 2-D HD inverse cascade of kinetic energy given by the two large-scale vortices, whereas the ion vorticity shows many filamented structures, typical of 2-D MHD turbulence. In this regime, we can approach the question of whether or not a one-fluid description is adequate for the proper determination of the dynamics by going into the centre-of-mass frame, typically done in other plasma settings (e.g. when deriving MHD itself). We define 
$\boldsymbol {V} \equiv \chi \boldsymbol {v}_i + (1-\chi ) \boldsymbol {v}_n$ and 
$\boldsymbol {D} = \boldsymbol {v}_i - \boldsymbol {v}_n$. The question now becomes: Is it possible only to account for the dynamics of 
$\boldsymbol {V}$ without knowing or integrating the dynamics of 
$\boldsymbol {D}$? Although it is not shown here, the simulations reveal that in the 
$\tilde {\alpha } \ll 1$ regime this is generally not possible – that is, the dynamics of 
$\boldsymbol {V}$ is partially determined by 
$\boldsymbol {D}$ and vice versa. However, as 
$\chi \rightarrow 0$ or 
$\chi \rightarrow 1$, we find that a one-fluid description (
$\boldsymbol {V}$ only) is sufficient, as one might expect. This was shown by comparing the relative magnitudes of 
$\boldsymbol {V}$ and 
$\boldsymbol {D}$ and noting when 
$|\boldsymbol {D}|\ll |\boldsymbol {V}|$.
Figure 1. Snapshots at steady state from K4 runs of the vorticity for each species, 
$\omega _s \equiv \boldsymbol {\nabla } \times \boldsymbol {v}_s$, representing, from left to right, 
$\tilde {\alpha } \ll 1$, 
$\tilde {\alpha } \sim 1$ and 
$\tilde {\alpha } \gg 1$. Panels (a)–(c) show the neutral vorticity for the three values of 
$\tilde {\alpha }$ whereas panels (d)–(f) show the ion vorticity for those runs.
As we increase 
$\tilde {\alpha }$ so that it is of order one, we see the neutral species begin to lose the large-scale vortices, which are dissipated away by collisional heating. In figure 1(b) and (e), at 
$\tilde {\alpha } = 3.15$, we no longer see obvious 2-D HD behaviour from the neutral species, but it also does not appear to be of a similar nature to the ion vorticity. Later scale-by-scale analysis will reveal that this regime is where the highest collisional heating is found and that most of the energy injected into the neutrals will be transferred to the ions or simply dissipated away. In figure 1(c) and (f), we see the case of 
$\tilde {\alpha } = 314.98$, and note immediately that the two fluids look identical from visual inspection. The two fluids are coupled and so the neutral species is behaving like an MHD fluid, as was predicted.
Although our predictions seem to agree qualitatively given the snapshots in figure 1, we now aim to confirm our results from § 2.2 in a quantitative way. In figure 2 we see two panels comparing global (time- and volume-averaged) quantities, where each data point is a single simulation. All runs performed in this study are included. Each panel aims to test the predictions made for each extreme of 
$\tilde {\alpha }$.
Figure 2. (a) Non-dimensional total energy versus the collision strength 
$\tilde {\alpha }$, rescaled by the energy injection rate of the neutrals over the ionization fraction. The red dashed line shows the prediction from (2.6), valid in the 
$\tilde {\alpha } \ll 1$ limit. (b) Non-dimensional collisional heating, rescaled by the square of the averaged Lorentz force, versus collision strength 
$\tilde {\alpha }$. The red dashed line shows the prediction from (2.8c), valid in the 
$\tilde {\alpha } \gg 1$ limit.
Figure 2(a) compares total energy with the collision strength rescaled by the energy injection rate of the neutral species over the ionization fraction, based on the right-hand side of (2.6), and with tildes representing suitable non-dimensionalization by a combination of 
$u_f$ and 
$k_f$. The red dashed line shows the prediction from (2.6), which is valid for small values of 
$\tilde {\alpha }$, and represents the balance between collisional heating and energy injection rate into the neutral species. The collapse of the data – for various values of ionization fraction, 
$Re$ and forcing wavenumber – on a single line, whose slope agrees with the predicted relationship over various orders of 
$\tilde {\alpha }$, confirms our claims.
Figure 2(b) focuses on the highly coupled regime, where the prediction for collisional heating was based on an expansion over 
$\tilde {\alpha }^{-1}$, revealing that 
$CH$ was given by the square of the difference in the forces acting on each species – given by (2.8c). This prediction is quite general but simplifies significantly for our simulations due to the fact that densities are uniform, 
$\nu _i = \nu _n$, and 
$\boldsymbol {F}_i = \boldsymbol {F}_n$. After non-dimensionalizing using 
$u_f$ and 
$k_f$, we are left with 
$\widetilde {CH} = (1-\chi )\chi ^{-1} \tilde {\alpha }^{-1} \widetilde {LF}$, where 
$\widetilde {LF} = \langle | \tilde {\boldsymbol {J}} \times \tilde {\boldsymbol {B}} |^2 \rangle$ and 
$\widetilde {(\cdot)}$ denotes the dimensionless version. Moving everything except 
$\tilde {\alpha }$ to the left-hand side, we get that, in the high collisional limit, the rescaled collisional heating (
$y$-axis of figure 2b) should be proportional to 
$\tilde {\alpha }^{-1}$, denoted by a red dashed line. We see that indeed the data collapse onto a single line, once again for various ionization fractions, 
$Re$ and forcing wavenumbers. The slope of the rescaled collisional heating seems to approach the theoretical one for large values of 
$\tilde {\alpha }$, confirming our predictions for the highly collisional limit.
4.2. Spatial
Our spatially averaged (‘global’) predictions for the extreme limits of 
$\tilde {\alpha }$, taken from § 2.2, have been confirmed. However, turbulence is an out-of-equilibrium, multi-scale process whose scale-by-scale analysis can reveal further interesting phenomena that are difficult to identify or predict otherwise. This is the purpose of the current subsection.
In our scale-by-scale analysis we will look at two general quantities – spectra, which tell how a quadratic quantity is distributed over scales, and fluxes, which tell us how that quadratic quantity is flowing through scales (Alexakis & Biferale Reference Alexakis and Biferale2018). The one-dimensional spectrum 
$KE_s(k)$ of the kinetic energy of a species 
$s$ is
\begin{equation} KE_s(k) = \frac{1}{2} \sum_{| \boldsymbol{k}|=k} | \widehat{\boldsymbol{v}_s}|^2( \boldsymbol{k}), \end{equation}
where 
$\widehat {(\,\cdot \,)}$ denotes the Fourier transform of 
$\boldsymbol {v}_s$. We will also be looking at the dimensionless collisional heating spectra,
\begin{equation} \widetilde{CH}(k) = \frac{(1-\chi) \chi \alpha}{u^3_f k_f} \sum_{| \boldsymbol{k}|=k} | \widehat{\boldsymbol{v}_i - \boldsymbol{v}_n}|^2( \boldsymbol{k}), \end{equation}and the dimensionless spectra of the square of the Lorentz force
\begin{equation} \widetilde{LF}(k) = \frac{1}{k_f^2 u_f^4} \sum_{| \boldsymbol{k}|=k} | \widehat{\boldsymbol{J} \times \boldsymbol{B}}|^2( \boldsymbol{k}). \end{equation} We have seen that 
$\tilde {\alpha }$ measures the relative strength between collisions and the nonlinear term at a typical scale 
$L$ with velocity 
$U$ (or, in the case of our simulations, 
$k^{-1}_f$ and 
$u_f$). Therefore, it measures the degree of coupling between the two fluids at 
$L$. However, due to the multi-time-scale nature of turbulence, different length scales couple at different values of 
$\tilde {\alpha }$. The global results from the last subsection showed us how to predict the average coupling of the fluid, given some information about the collision strength at the forcing scale. However, a scale-by-scale analysis allows us to see exactly how the two fluids gradually couple together among scales and to better understand non-extreme cases for 
$\tilde {\alpha }$. In figure 3 we see the steady-state average kinetic energy spectra of both species, for three different values of 
$\tilde {\alpha }$, as in figure 1. These spectra have ionization fraction 
$\chi = 0.5$ and were taken from the K8 set of runs, forced at intermediate wavenumbers so as to be able to resolve approximate inertial ranges in both large and smaller scales.
Figure 3. Kinetic energy spectra for the neutral species (orange, solid line) and ion species (green, dashed line) at (a) 
$\tilde {\alpha } = 0.002$, (b) 
$\tilde {\alpha } = 2.520$ and (c) 
$\tilde {\alpha } = 29.923$, all with ionization fraction 
$\chi =0.5$. As 
$\tilde {\alpha }$ increases, the two species become more and more coupled, starting from the largest scales which couple first.
In figure 3(a) we have the uncoupled limit, evident by the two spectra and their distinct shapes. The black dot-dashed lines show a 
$-5/3$ slope, the prediction for the large scales of 
$KE_n$ and the small scales of 
$KE_i$. The slope of the neutral kinetic energy is a bit steeper at large scales due to the presence of a condensate. Despite this, as well as a poorly resolved inertial range for the ion kinetic energy, we see reasonable agreement between expected and observed behaviour for each individual species. Figure 3(b) shows 
$\tilde {\alpha } = 2.52$, a moderately coupled case. In this figure we see that, indeed, the energies at scales down to the forcing scale seem to be lying on top of each other, implying a coupling at those scales. The collisions are strong enough at these scales to completely remove the inverse cascade of kinetic energy in the neutral species, either by dissipation or by transfer of energy to the ions (to be discussed further when we look at the fluxes). However, at scales smaller than the forcing there is a clear distinction between the two energies, implying a lack of coupling. It is this remaining ‘slippage’ between the ions and the neutrals in the small scales, along with the fact that 
$\tilde {\alpha } \sim \textit{O}(1)$, that maintains a large collisional heating, despite the small value of 
$|\boldsymbol {v}_i-\boldsymbol {v}_n|$ at large scales. In fact, collisional heating is largest at these values of 
$\tilde {\alpha }$, as seen in figure 2(b).
Finally, figure 3(c) shows a higher collisional case with 
$\tilde {\alpha } = 29.923$. In this figure we observe that scales smaller than the forcing scale are now beginning to couple; however, 
$\tilde {\alpha }$ is small enough so that the smallest scales are still not coupled. The two fluids look almost identical, like the same MHD fluid. At this point a natural question arises: At what scales do the fluids decouple and how does that depend on 
$\tilde {\alpha }$? We have attempted to numerically investigate the dependence of the decoupling wavenumber, which we are calling 
$k_{coll}$, on 
$\tilde {\alpha }$, but were not able to reach a definitive conclusion. Possible issues include: a small inertial range (limited by the numerical capabilities); a failure of time-scale assumptions common in homogeneous, isotropic turbulence which only are valid in a true inertial range (here, dissipation due to collisions might render that irrelevant); and more. We do, however, think that the results are worth showing, and so we have included this work in appendix B.
A quantity intimately related to the degree of coupling between the two species is the collisional heating, which is proportional to 
$|\boldsymbol {v}_i - \boldsymbol {v}_n|^2$. A scale-by-scale decomposition of 
$CH$ can tell us where kinetic energy dissipation (heating) due to collisions is most active, which might be of interest in astrophysical applications. In figure 4 we plot the non-dimensional spectra of 
$CH$, the black dashed curves, for the usual three cases of 
$\tilde {\alpha }$. These are taken from steady-state averages in the K32 set of runs, which were forced at small scales.
Figure 4. The spectra of dimensionless collisional heating 
$\widetilde {CH}$ (bold, black, dashed lines) and of the dimensionless right-hand side of (2.8c), proportional to the spectrum of the Lorenz force squared, 
$\widetilde {LF}$ (bold, light blue, solid lines). These two curves are expected to be equal when the two species are highly coupled. Shown are three limits of 
$\tilde {\alpha }$, namely (a) weak coupling, (b) moderate coupling and (c) strong coupling, allowing one to see the gradual coupling of the species and the change in form of the 
$CH$ spectrum.
Focusing only on the black dashed curves for now, we look first at figure 4(a), the low coupling case. Since the two species are not coupled, we can approximately say that, at scales larger than the forcing, the kinetic energy of the neutrals dominates and hence 
$CH \propto |\boldsymbol {v}_i-\boldsymbol {v}_n|^2 \approx |\boldsymbol {v}_n|^2 \propto KE_n$. This is confirmed by the 
$-5/3$ slope shown by the thin dot-dashed black line. Hence, the collisional heating is acting like a friction term for the neutral species, causing significant dissipation at the largest scales where the condensate lies. Although not shown here, apart from directly dissipating energy to heat, collisions also act to transfer energy from neutrals to ions at large scales (particularly for low ionization fractions). This possibility is evident by looking at (2.2a) – when 
$|\boldsymbol {v}_n|\gg |\boldsymbol {v}_i|$, the sign-indefinite term will dominate 
$|\boldsymbol {v}_i|^2$ for the ion species kinetic energy equation, thus providing the possibility to gain kinetic energy via collisions with the highly energetic neutrals at those scales. The species first begin to couple at the largest scales, as we saw when looking at the kinetic energy spectra. This coupling reduces 
$CH$, causing the length scale of maximal collisional heating to become smaller and smaller as we increase 
$\tilde {\alpha }$.
This is exemplified in figure 4(b), where we already see a significant change in the shape of the collisional heating spectrum, its peak around 
$k = 15$. This happens until 
$\tilde {\alpha } > 1$, after which the maximal heating is at the forcing scale and the shape of the 
$CH$ spectrum remains practically unchanged, with only its magnitude being reduced (proportional to 
$\tilde {\alpha }^{-1}$). The shape of the collisional heating spectrum at scales that are already coupled is set by the spectrum of 
$\langle |\boldsymbol {J} \times \boldsymbol {B}|^2 \rangle$, as predicted by (2.8c). The light blue, solid lines in figure 4 represent the spectra of the right-hand side of this equation, and indeed the two curves lie on top of each other at the scales that are coupled. Thus, understanding the spectrum of the square of the Lorentz force is key to understanding the scale-by-scale collisional heating. The small dark blue dashed line depicting a slope of 
$2.4$ (found by fitting that region of the spectrum) can be seen in figure 4(a) and (c), showing that the spectrum of the Lorentz force has not changed and is not affected by the collisions. One would still like to know what sets the spectrum of the Lorentz force (and hence of 
$CH$). A very clear power law with a positive exponent is seen in the figure, so one is tempted to speculate about its origins, particularly using equilibrium spectra arguments, as is done in 3-D homogeneous and isotropic turbulence. However, such equilibrium spectra are not necessarily universal if the forcing in the spectral shell around 
$k_f$ is dense, as is ours (Alexakis & Brachet Reference Alexakis and Brachet2019). A prediction of the spectral slope of the Lorentz force, although interesting and relevant for our applications, is beyond the scope of this paper.
The spectra have revealed the length scales at which coupling between the species occurs, as well as allowed us to observe at what scales collisional heating dominates and confirm our predictions about the spectral shape of 
$CH$. While these spectra have informed us of the distribution of various quantities among length scales, also of interest in turbulence is how certain conserved quantities, such as energy, move across scales. This is measured by the (spectral) flux. In this study, we will only look at two components of energy flux, what we call the kinetic component of energy flux for a species 
$s$, 
$\Pi _{U_s}$, and the magnetic component of energy flux for the ions, denoted by 
$\Pi _{B}$. The former is defined to be
\begin{equation} \Pi_{U_s}(k) = \langle \boldsymbol{v}^{{<}k}_s \boldsymbol{\cdot} (\boldsymbol{v}_s \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}_s) \rangle, \end{equation}
where 
$\boldsymbol {v}^{<k}_s$ stands for a filtering of the velocity 
$\boldsymbol {v}_s$ in Fourier space so that only the wavenumbers with modulus smaller than 
$k$ are kept. The flux 
$\Pi _{U_s}$ expresses the rate at which kinetic energy of a species 
$s$ is flowing out of scales larger than 
$\ell = 2 {\rm \pi}/ k$ due to nonlinear interactions only in velocity. Therefore, if energy is going from large to small scales, the energy flux will be positive, and vice versa. The magnetic component of energy flux for the ion species is defined to be
\begin{equation} \Pi_{B}(k) ={-} \langle \boldsymbol{v}_i^{{<}k} \boldsymbol{\cdot} (\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{B} ) \rangle + \langle \boldsymbol{B}^{{<}k} \boldsymbol{\cdot} ( \boldsymbol{v}_i \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{B} - \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}_i ) \rangle. \end{equation} As a reminder, in 2-D HD turbulence, energy flows to larger scales whereas 2-D MHD turbulence has been shown to cascade total energy forward to smaller scales. Thus, in an uncoupled system with 
$\tilde {\alpha } = 0$, we expect the neutral kinetic energy flux 
$\Pi _{U_n}$ to be negative at scales larger than the forcing and zero at scales smaller than the forcing. On the other hand, 
$\Pi _{U_i} + \Pi _{B}$ ought to be zero at large scales and positive at scales smaller than the forcing. Although 
$KE_i$ and 
$E_B$ are not individually conserved, the nonlinear interactions that cause the respective fluxes 
$\Pi _{U_i}$ and 
$\Pi _{B}$ each conserve energy and can thus also be analysed.
In figure 5 we look at all four fluxes, 
$\Pi _{U_n}$, 
$\Pi _{U_i}$, 
$\Pi _{B}$ and 
$\Pi _{U_i}+\Pi _{B}$, normalized by the total energy injection rate 
$I_{tot} = I_n + I_i + I_B$. We look at three cases of 
$\tilde {\alpha }$, the first two taken from the K8 set of runs, and the large 
$\tilde {\alpha }$ limit from the K8R run (all with 
$\chi = 0.5$). These fluxes are averaged over the steady state. The 
$\tilde {\alpha } \ll 1$ case, figure 5(a), confirms the uncoupled predictions. The neutral kinetic energy flux (orange solid line) is negative and roughly equal to the neutral injection rate: 
$I_n/ I_{tot} = 4/9$ (recall that 
$I_n = I_i = 4 I_B$). The total ion energy flux (black solid line) is positive and seen to be close to 
$(I_i+I_B)/I_{tot} = 5/9$. Since 
$\tilde {\alpha } \neq 0$, we do see a hint of some large-scale energy dissipation in the neutrals, given by the slight negative slope of the flux at large wavenumbers.
Figure 5. Kinetic component of energy flux for the neutral species, 
$\Pi _{U_n}$ (orange, solid line), and ion species, 
$\Pi _{U_i}$ (green, dashed line), as well as the magnetic component of energy flux for the ions, 
$\Pi _B$ (blue, dot-dashed line). Total ion energy flux is denoted by the solid, black line, and grand total energy flux 
$(\Pi _{U_n}+\Pi _{U_i}+\Pi _{B})$ by the thin, solid, red line, seen only in panel (c). The three panels depict the three limiting cases for 
$\tilde {\alpha }$: (a) uncoupled, where we see the HD and MHD predictions for fluxes; (b) moderate coupling, which is highly dissipative; and (c) high coupling, in which the neutrals are fully coupled to the ions and we recover the prediction from § 2.2 that the system behaves as a single MHD fluid with 
$\boldsymbol {F} = \boldsymbol {F}_i + \boldsymbol {F}_n$.
Turning now to the case where 
$\tilde {\alpha } \approx 1$, figure 5(b) shows that the neutral kinetic energy flux is practically zero for all wavenumbers, meaning that collisions have almost completely stopped any nonlinear transfer of kinetic energy among the neutrals. At the forcing scale and larger, the two species are coupled and behave like MHD, as we have seen. Since MHD has no inverse cascade of kinetic energy, the neutrals have no inverse cascade at these scales. At scales smaller than the forcing, the species are not yet coupled, and since the neutral species does not have a forward cascade of energy, the flux remains zero. However, there is still constant injection rate of kinetic energy into the neutrals of magnitude 
$I_n$, and this energy must go somewhere. That energy is either dissipated away by collisions, or transferred to the ion species. Seeing as, at best, approximately half of the total energy injection rate is being transferred spectrally by the ions, it seems that in this case approximately half of the energy injected at the forcing scale is immediately dissipated by collisions, the rest being transferred away by the ion species. For the case of uncoupled 2-D MHD turbulence, one would expect this flux of energy by the ion species and magnetic field to be constant in the inertial range down to the dissipation scales. However, in figure 5(b) we see a gradual drop-off of energy flux. This non-constant flux is a sign of strong dissipation of energy in the would-be inertial range due to collisional heating, which is largest in the runs where 
$\tilde {\alpha } \sim \textit{O}(1)$, as we have seen.
Moving on to figure 5(c), we notice that the ion and neutral species are fully coupled, even at the small scales, as seen by the fact that 
$\Pi _{U_n} = \Pi _{U_i}$. Now that the two species are coupled throughout most of the scales, the collisional heating is very weak and does not dissipate practically any energy that is being injected at the forcing scale. Figure 5(c) indicate that all of the energy being injected by the forcing – including the neutral energy injection – is going to small scales. This is shown by the total energy flux from all fields, denoted by the thin, red, solid line, which we see is positive and very close to 1. This is consistent with our 
$\tilde {\alpha } \rightarrow \infty$ limit studied in § 2.2, where we showed that the system behaves as a single MHD fluid with 
$\boldsymbol {F} = \boldsymbol {F}_i + \boldsymbol {F}_n$, as seen in (2.8a). Sticking to 
$\boldsymbol {v}_i$ and 
$\boldsymbol {v}_n$, what seems to be happening is that the two species exchange energy between each other very efficiently due to collisions, but, since ions also exchange energy with the magnetic field, part of the energy being injected into the neutral species ends up as magnetic energy, which is transferred to smaller scales. Although it is not shown here, we ran a few decaying turbulence experiments with no forcing to see if this description holds true. Indeed, after initializing the two species with identical random initial conditions (and 
$KE_i = KE_n = E_B$ at 
$t = 0$), the runs with very large 
$\tilde {\alpha }$ showed larger initial magnetic field growth, implying an indirect transfer of energy from the neutrals.
5. Discussion and conclusions
5.1. Conclusions
In this work we have investigated the 2-D partially ionized magnetohydrodynamic (PIMHD) system, a two-fluid model used for studying plasmas where some fraction of the ions have recombined to form neutral molecules that interact with the other ions via collisions. Although ionization fraction certainly plays a role in the dynamics of such plasmas, we focused more on the role that collisions have on the dynamics. In the limit where collision time scales are long compared to dynamical time scales, we found that the two species were weakly coupled and behaved like their uncoupled fluid counterparts – hydrodynamics (HD) for neutral species and magnetohydrodynamics (MHD) for the ions. In this limit, collisions act like a frictional force for whichever species has larger kinetic energy at that scale. Hence, collisions took the role of a large-scale dissipation for the neutral species. Using this observation and an energy balance argument we were able to predict the amount of collisional heating in these runs. For low ionization fractions, it is possible that ions may gain some kinetic energy via collisions at these large scales. At intermediate collision strengths, where dynamic and collision time scales are similar, we found the dynamics to be quite dissipative. Approximately half of the energy injected into the system is dissipated immediately at the forcing scale, and another quarter of it is dissipated by collisions along the forward cascade range. In this case, the neutral species has little to no nonlinear energy transfers.
For collision time scales much shorter than dynamical time scales we found that the two species are coupled and act like a single MHD fluid whose density, pressure, viscosity and external forces are the sum of each from the two species. This was confirmed numerically, particularly in figure 5(c), which showed all of the energy injected at the forcing scale being transferred to smaller scales. This species coupling in turn reduces the collisional heating significantly, which we showed, and also confirmed numerically, is determined by the square of the difference of accelerations of each species. In our specific implementation, the nature of the Lorentz force meant that the collisional heating no longer dominates at the largest scales (as in the low collisional case) but is present at smaller scales. A scale-by-scale analysis of our runs allowed us to further understand how the coupling of the two species gradually occurs as the collision time scale decreases. Although we numerically investigated 2-D PIMHD, our results for the strong collisional case are expected to hold for the 3-D counterpart since dimensionality did not play a role in the derivation of the predictions made.
5.2. Connection to Jupiter
This work was motivated by the transition region between the ionized interior and neutral upper atmosphere of gas giant planets such as Jupiter, where partial ionization effects are expected to be present. The current treatment of the dynamics in the transition region have been single-fluid models, either of MHD with large resistivity, or HD with a drag term to parametrize any magnetic effects (‘MHD drag’). Given our characterization of the behaviour of such a two-fluid plasma, we will attempt to apply some of what we learned here to the transition region of gas giant planets. In order to do so we must first estimate the ratio of eddy turnover time to collision time scale, 
$\tilde {\alpha }$, for the case of Jupiter's transition region, whose densities, temperatures, pressures and typical velocities and length scales are better constrained than any other gas giant planet.
Following (2.4), we must choose a typical velocity 
$U$ and length scale 
$L$, as well as total density 
$\rho _{tot}$, and finally the collision strength 
$\alpha$. Typical velocities for the transition region can range from 
$U = 10^{-2}~\textrm {m}\,\textrm {s}^{-1}$ close to the edge of the ionized interior, at a radius of roughly 
$0.9$ Jupiter radii, to 
$U = 10^2~\textrm {m}\,\textrm {s}^{-1}$ at the surface (Kaspi et al. Reference Kaspi, Galanti, Hubbard, Stevenson, Bolton, Iess, Guillot, Bloxham, Connerney and Cao2018). We take the typical length scale in this region from Cao & Stevenson (Reference Cao and Stevenson2017), who estimated a convective length scale of 
$L = 6 \times 10^5~\textrm {m}$. The total density is found from simulations by French et al. (Reference French, Becker, Lorenzen, Nettelmann, Bethkenhagen, Wicht and Redmer2012), who give values of 
$\rho _{tot} = 10^3~\textrm {kg}\,\textrm {m}^{-3}$ at 
$0.9$ Jupiter radii to 
$\rho _{tot} = 2 \times 10^2~\textrm {kg}\,\textrm {m}^{-3}$ closer to the surface.
The final piece is the collision strength 
$\alpha$, whose value is much harder to constrain, given the fact that, to date, no plasma experiment has been able to measure collision cross-sections at such high densities. Therefore, any value used in our estimate will be questionable, but we feel it would be best to give some estimate rather than none. Following § 5.2.1 in Meier (Reference Meier2011) and § III.D.2 from Meier & Shumlak (Reference Meier and Shumlak2012), under the assumptions that ions and neutrals have the same temperature and that both are in thermal equilibrium, we get the following expression for 
$\alpha$:
\begin{equation} \alpha = \alpha_1 \sqrt{T} - \alpha_2\ln \left(205 \sqrt{T} \right)\sqrt{T}, \end{equation}
where 
$\alpha _1 = 1.34 \times 10^{11}~\textrm {m}^2\,\textrm {kg}^{-1}\,\textrm {s}^{-1}\,\textrm {K}^{-1/2}$ and 
$\alpha _2 = 8.78 \times 10^9~\textrm {m}^2\,\textrm {kg}^{-1}\,\textrm {s}^{-1}\,\textrm {K}^{-1/2}$. We should note that these numbers are taken from cross-section estimates for charge exchange reactions (e.g. an electron hopping from a neutral molecule to an ionized molecule), which are expected to be dominant even at the relatively low (sub-
$10\,000~\textrm {K}$) temperatures considered here (E. T. Meier, personal communication), rather than more standard collisions. However, the functional form of the collision term does not change. To get an estimate for 
$\alpha$ at 0.9 Jupiter radii and near the surface, we use the temperature profile from French et al. (Reference French, Becker, Lorenzen, Nettelmann, Bethkenhagen, Wicht and Redmer2012), giving approximately 
$5000~\textrm {K}$ and 
$1000~\textrm {K}$, respectively. This in turn results in a range of 
$\alpha$ between 
$3.4 \times 10^{12}~\textrm {m}^2\,\textrm {kg}^{-1}\,\textrm {s}^{-1}$ in the interior and 
$1.8 \times 10^{12}~\textrm {m}^2\,\textrm {kg}^{-1}\,\textrm {s}^{-1}$ closer to the surface. Previous studies on shock waves in the interstellar medium (Smith & Mac Low Reference Smith and Mac Low1997) and for hot Jupiters (Perna et al. Reference Perna, Menou and Rauscher2010), have obtained similar values. Owing to the temperature dependence, and the relatively small size of the transition region in Jupiter, the collision strength does not change much. In a planet such as Saturn, with a much larger transition region (Cao & Stevenson Reference Cao and Stevenson2017), this might not be the case.
Combining everything, we are able now to give an estimate for the range of 
$\tilde {\alpha }$: closer to 0.9 Jupiter radii we estimate 
$\tilde {\alpha } = 2 \times 10^{23}$, whereas towards the surface we estimate 
$\tilde {\alpha } = 2 \times 10^{18}$. Given the 
$\tilde {\alpha }^{-1}$ dependence of 
$CH$ from (2.8c), this would imply an incredibly small (and presumably negligible) amount of collisional heating in the transition region of Jupiter, as well as single-species MHD-like dynamics, with two-fluid effects being negligible. In fact, 
$\alpha$ is so large that the two species would be coupled (and therefore behave as an MHD fluid) even for extremely low ionization fractions – since 
$\tau _{eddy}\nu _{in} = L \rho _i \alpha U^{-1}$ would be very large. This therefore seems to validate the single-fluid MHD models used so far in dynamo studies which include the transition zone (Jones et al. Reference Jones, Boronski, Brun, Glatzmaier, Gastine, Miesch and Wicht2011; Gastine et al. Reference Gastine, Wicht, Duarte, Heimpel and Becker2014; Jones Reference Jones2014; Dietrich & Jones Reference Dietrich and Jones2018). We want to emphasize, however, that, although the single-species MHD description holds throughout most of the transition region, whether the fluid behaves as an MHD or HD fluid depends on whether the Lorentz force is significant or not. Since we also expect 
$\eta _{MHD}$ to depend on temperature, it is quite possible that in the transition region the single-species description, although technically an MHD fluid, does not feel the Lorentz force and thus behaves as an HD system. Thus, the transition from MHD to HD could happen continuously, although not because neutrals make up the majority of the fluid, but because the Lorentz force becomes weaker and weaker as the magnetic diffusivity significantly increases as a function of radius, leading to neutral-like behaviour of the MHD fluid (Tobias, Diamond & Hughes Reference Tobias, Diamond and Hughes2007; Seshasayanan, Benavides & Alexakis Reference Seshasayanan, Benavides and Alexakis2014; Seshasayanan & Alexakis Reference Seshasayanan and Alexakis2016). As this transition happens, we go into the regime of low-magnetic-Reynolds-number MHD turbulence, which is where the MHD drag prescription is applied.
Acknowledgements
The authors would like to thank K. J. Burns, J. S. Oishi, B. Gallet, D. D. B. Koll and E. T. Meier for many helpful discussions throughout the development of the project. This project was funded by a grant from the National Science Foundation (OCE-1459702).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Limiting cases of ionization fraction 
$\chi$
Note that no longer having both 
$\chi \sim \textit{O}(1)$ and 
$(1-\chi ) \sim \textit{O}(1)$ now means that the two species feel collisions to a different degree, since, as we saw in § 2.2, 
$\tau _{coll,i} = \tau _{coll}/(1-\chi )$ and 
$\tau _{coll,n} = \tau _{coll}/\chi$. In the fully ionized limit, 
$\chi \rightarrow 1$, the collisional time scale in the ion equation, 
$\tau _{coll,i}$, becomes very large, yet for neutrals 
$\tau _{coll,n}$ does not change. This results in the ions not feeling collisions in the low collisional limit, whereas neutrals are still affected. In the high collisional limit, the collision term is still at least of order one for both species; however, the highest order dynamics are slightly different. Away from the dissipation range, a similar expansion in 
$\tilde {\alpha }^{-1}$ results in 
$\boldsymbol {v}_n^{(0)} = \boldsymbol {v}_i^{(0)} + (\boldsymbol {F}_n -\boldsymbol {\nabla } p_n)/ (\rho _i \rho _n \alpha )$ so that the order-one dynamics is simply the MHD of 
$\boldsymbol {v}_i$ but with a modified pressure 
$P = p_i + p_n$ and force 
$\boldsymbol {F} = \boldsymbol {F}_i+\boldsymbol {F}_n$.
In the low ionization limit, 
$\chi \rightarrow 0$, it is now the neutrals that are less affected by collisions. This limit is a bit more involved because it encompasses the breakdown of the validity of other approximations in the derivation of the MHD equations, regardless of the value of 
$\tilde {\alpha }$. There are two parameters of interest to recall: the electron-to-ion mass ratio, 
$\beta \equiv m_e/m_i$, and the ratio of ion skin depth, 
$d_i \equiv \sqrt {c^2 m_i/(4 {\rm \pi}e^2 n_i)}$, to typical length 
$L$, which we call 
$\lambda$. Here 
$m_e$ and 
$m_i$ are the masses of the electron and ion, respectively, 
$c$ is the speed of light, 
$e$ is the electric charge of the electron and 
$n_i = \rho _i/m_i$. In the derivation of MHD, these two parameters are taken to be much smaller than unity. However, when deriving PIMHD, ratios of 
$\chi$, 
$\lambda$ and 
$\beta$ appear and so, when 
$\chi \ll 1$, terms that would normally be negligible in MHD are no longer small. It can be shown that, if 
$\chi \sim \lambda \sim \beta \ll 1$, then the induction equation (2.1c) must be modified to include the Hall term, so that the nonlinear term becomes 
$\boldsymbol {\nabla } \times (\boldsymbol {v}_i \times \boldsymbol {B} - d_i (\boldsymbol {J} \times \boldsymbol {B})/\rho _i)$ (Pandey & Wardle Reference Pandey and Wardle2008). Note that this is now true even at large scales, not just small scales, as is normally the case. In the low collisional limit, the neutrals do not feel the collisions whereas the ions do. On the other hand, when 
$\tilde {\alpha } \gg 1$, we have, to highest order, 
$\boldsymbol {v}_i^{(0)} = \boldsymbol {v}_n^{(0)} +(\boldsymbol {J} \times \boldsymbol {B} - \boldsymbol {\nabla } p_i + \boldsymbol {F}_i) / (\rho _i \rho _n \alpha )$, leading to what is typically called ‘ambipolar MHD’ in the astrophysics community:
\begin{gather} \rho_{n} \left( \frac{\partial}{\partial t} + \boldsymbol{v}_n \boldsymbol{\cdot} \boldsymbol{\nabla}\right) \boldsymbol{v}_n ={-} \boldsymbol{\nabla} (p_n+p_i) + \boldsymbol{J} \times \boldsymbol{B} + \boldsymbol{F}_i+ \boldsymbol{F}_n, \end{gather}
\begin{gather} \frac{\partial \boldsymbol{B}}{\partial t} = \boldsymbol{\nabla} \times \left( \boldsymbol{v}_n \times \boldsymbol{B} - \frac{d_i}{\rho_i}\boldsymbol{J} \times \boldsymbol{B} + \frac{1}{\rho_i \rho_n \alpha} ( \boldsymbol{J} \times \boldsymbol{B} - \boldsymbol{\nabla} p_i + \boldsymbol{F}_i) \times \boldsymbol{B} \right). \end{gather}Counter to intuition, despite the neutrals making up the majority of the density in the fluid, the neutrals indirectly interact with the magnetic field via their collisional interaction with the ions and the dynamics is still that of MHD. It is only when 
$\chi \ll \lambda \sim \beta \ll 1$ that the behaviour of the fluid reverts to regular HD without interaction with the magnetic field. We should note here that we have treated 
$\chi$ and 
$\eta _{MHD}$ from (2.3) independently and have effectively held the latter fixed while varying 
$\chi$. In our claim that it is possible to have MHD behaviour even for 
$\chi \ll 1$, we have assumed that the magnetic diffusivity is low enough so that the Lorentz force is still of order one. In a realistic situation, it is likely that ionization fraction and idealized magnetic diffusivity 
$\eta _{MHD}$ are both related to temperature and so, as 
$\chi \rightarrow 0$, we would also expect 
$\eta _{MHD}$ to increase significantly; therefore, it is possible that the Lorentz force is negligible when 
$\chi \sim \beta \sim \lambda \ll 1$, making the dynamics that of HD (Tobias et al. Reference Tobias, Diamond and Hughes2007; Seshasayanan et al. Reference Seshasayanan, Benavides and Alexakis2014; Seshasayanan & Alexakis Reference Seshasayanan and Alexakis2016). However, this is a separate effect than that of ionization fraction, and our claims above still hold in general.
Appendix B. Determining $k_{coll}$ versus $\tilde {\alpha }$
A typical argument would tell us that, when the eddy turnover time at a scale 
$k_{coll}^{-1}$, 
$\tau _{eddy}(k_{coll})$, is equal to the collision time, 
$\tau _{coll} = (\rho _{tot}\alpha )^{-1}$, then that length scale is coupled. This sort of argument has been used before for estimating decoupling scales for both MHD waves and nonlinear dynamos in a partially ionized system (Xu et al. Reference Xu, Yan and Lazarian2016; Xu & Lazarian Reference Xu and Lazarian2017). This time-scale argument relies on an important assumption: that the velocity scales in the nonlinear (or wave) term and in the collision term have similar orders of magnitude independent of 
$\tilde {\alpha }$, i.e. 
$k_{coll} U^2 \sim \rho _{tot} \alpha U$, making 
$U$ cancel out and leaving one with a ratio of time scales. This assumption is equivalent to stating that 
$|\boldsymbol {V}| \sim |\boldsymbol {D}|$ (defined in § 4.1). However, the results of § 2.2 tell us that 
$|\boldsymbol {D}| = |\boldsymbol {v}_i-\boldsymbol {v}_n| \propto \alpha ^{-1}$, whereas 
$|\boldsymbol {V}| \sim \textit{O}(1)$. Although this scaling with 
$\tilde {\alpha }$ is only true for scales that are already highly coupled, we only expect 
$|\boldsymbol {D}|\sim |\boldsymbol {V}|$ far from the coupling scale, thus making the scaling non-trivial. The difficulty arises in the fact that, by definition, at the coupling scale, we are not in any extreme limit of 
$\tilde {\alpha }$ and so we cannot make a general statement about how 
$|\boldsymbol {D}|$ scales with 
$\tilde {\alpha }$. Furthermore, although it is not shown, we note that the spectra for 
$|\boldsymbol {V}|$ and 
$|\boldsymbol {D}|$ do not have the same power-law exponents, implying different 
$\tau _{eddy}$ tendencies with 
$k$, further complicating any length-scale analysis for finding 
$k_{coll}$. We therefore revert to simply observing the decoupling wavenumber 
$k_{coll}$ versus 
$\tilde {\alpha }$ in our simulations. We define 
$k_{coll}$ such that 
$|\boldsymbol {D}|(k_{coll})/|\boldsymbol {V}|(k_{coll}) = \delta$, where we have the freedom to choose 
$\delta$ as long as it is small. This is equivalent to saying that the fluid is coupled when 
$|\boldsymbol {v}_i - \boldsymbol {v}_n| \ll |\chi \boldsymbol {v}_i + (1-\chi )\boldsymbol {v}_n|$, with the freedom to define just how much smaller.
In figure 6 we see 
$k_{coll}$ versus 
$\tilde {\alpha }$ for two cases: (a) 
$k_{coll}<k_f$ and (b) 
$k_{coll}>k_f$. For each case, a 
$\delta$ was chosen so that the maximum number of runs could have 
$k_{coll}$ within the inertial range. Only these runs were kept in the figure. The former case is taken from various runs in the K32 set, which is forced at small scales to be able to better resolve the inertial range at scales larger than the forcing. The latter is taken from various runs in the K4 set, which is forced at larger scales to be able to better resolve the inertial range at scales smaller than the forcing. Notice that different values of 
$\delta$ are used for each case, but that 
$\delta <1$ for both cases. Owing to the limited scale separation possible in numerical experiments, this restricted the possible values of 
$\delta$ and made it difficult to test the robustness of the calculations of 
$k_{coll}$ as one varies 
$\delta$. The dashed lines represent best-fit power laws, whose exponents are seen in the legends of each panel. Although the power laws are very close to integer values, this may be a coincidence and we do not claim any physical significance, although we also have no reason to reject any. We see that the scaling of 
$k_{coll}$ with 
$\tilde {\alpha }$ depends strongly on which inertial range one is looking at, a reflection of the difference of behaviour in the two inertial ranges. Our results tell us that the scales larger than the forcing couple more slowly than those at scales smaller than the forcing.
Figure 6. Decoupling wavenumber 
$k_{coll}$ versus collision strength 
$\tilde {\alpha }$ for (a) 
$k_{coll}<k_f$ and (b) 
$k_{coll}>k_f$. We define 
$k_{coll}$ such that 
$|\boldsymbol {D}|(k_{coll})/|\boldsymbol {V}|(k_{coll}) = \delta$. The dashed lines represent the best-fit power law, whose exponent is shown in the legend.
