The formation and evolution of reconnectiondriven, slowmode shocks in a partially ionised plasma
^{1}
Department of Applied Mathematics and Theoretical PhysicsUniversity of
Cambridge,
Wiberforce Road,
Cambridge
CB3 0WA,
UK
email:
ah826@cam.ac.uk
^{2}
Kwasan and Hida Observatories, Kyoto University,
6068501
Kyoto,
Japan
Received: 29 January 2016
Accepted: 23 March 2016
The role of slowmode magnetohydrodynamic (MHD) shocks in magnetic reconnection is of great importance for energy conversion and transport, but in many astrophysical plasmas the plasma is not fully ionised. In this paper, we use numerical simulations to investigate the role of collisional coupling between a protonelectron, chargeneutral fluid and a neutral hydrogen fluid for the onedimensional (1D) Riemann problem initiated in a constant pressure and density background state by a discontinuity in the magnetic field. This system, in the MHD limit, is characterised by two waves. The first is a fastmode rarefaction wave that drives a flow towards a slowmode MHD shock wave. The system evolves through four stages: initiation, weak coupling, intermediate coupling, and a quasisteady state. The initial stages are characterised by an overpressured neutral region that expands with characteristics of a blast wave. In the later stages, the system tends towards a selfsimilar solution where the main drift velocity is concentrated in the thin region of the shock front. Because of the nature of the system, the neutral fluid is overpressured by the shock when compared to a purely hydrodynamic shock, which results in the neutral fluid expanding to form the shock precursor. Once it has formed, the thickness of the shock front is proportional to ξ_{ i}^{1.2} , which is a smaller exponent than would be naively expected from simple scaling arguments. One interesting result is that the shock front is a continuous transition of the physical variables of subsonic velocity upstream of the shock front (a cshock) to a sharp jump in the physical variables followed by a relaxation to the downstream values for supersonic upstream velocity (a jshock). The frictional heating that results from the velocity drift across the shock front can amount to ~2 per cent of the reference magnetic energy.
Key words: magnetic reconnection / magnetohydrodynamics (MHD) / shock waves
© ESO, 2016
1. Introduction
Magnetic reconnection is the change in the connectivity of a magnetic field and is responsible for violent energy release in many space and astrophysical systems. There are two key physical processes that are believed to be important in the energy conversion between magnetic energy and fluid (including thermal and kinetic) energy in magnetic reconnection: Joule heating and the work of the magnetic field on the fluid postreconnection. The latter of these mechanisms was proposed by Petschek (1964), and was shown to be highly effective at converting the magnetic energy to kinetic and internal energies. This is the case because this model was able to produce fast magnetic reconnection (roughly speaking the ratio of the inflow to outflow velocity is approximately 0.01–0.1 and only weakly dependent on the Lundquist number) as a result of the presence of standing slowmode shocks created by magnetic field relaxation. Priest & Forbes (1986) developed a general class of steadystate twodimensional (2D) reconnection models that included both the diffusiondominated reconnection region and the convective region that surrounds it. They found that both the Petschek (1964) and Sonnerup (1970) models were limiting cases of the general class of models they found. One key unifying feature of all the models was the presence of slow shocks associated with the relaxation of the magnetic field in the reconnection outflow region. Nevertheless, it should be noted that gradients and rotation in the magnetic field can produce rotational discontinuities (Petschek & Thorne 1967).
In this paper we are interested in the oblique form of the slowmode magnetohydrodynamic (MHD) shock known as the switchoff shock; the magnetic field is more perpendicular to the shock front downstream of the shock than upstream. This is associated with the removal of the component of the magnetic field parallel with the shock front, but the component perpendicular to the shock front does not change owing to the requirement for ∇·B. This results in a decrease in magnetic energy across the shock front that is balanced by an increase in both the thermal and kinetic energies of the fluid. The abrupt change in direction of the magnetic field results in the acceleration of the fluid with a strong component of the velocity parallel to the shock front. This can be viewed as analogous to the reconnection jet in the reconnection models cited in the previous paragraph.
Though there is still great debate about the role of slowmode shocks in creating fast magnetic reconnection, their importance in driving observed dynamics of the solar atmosphere should not be underestimated. Solar flares, which are huge energy releases in the solar atmosphere, are one of the observational signatures of magnetic reconnection in the solar corona (for a recent review see, for example, Shibata & Magara 2011). In the case of flare reconnection, thermal conduction also plays an important role for energy transport and results in isothermal slowmode shocks forming (e.g. Longcope & Bradshaw 2010; Takasao et al. 2015).
Another example that highlights the importance of reconnectiondriven, slowmode shocks was provided by Takasao et al. (2013), who performed a detailed 2D numerical investigation of reconnection associated with emerging magnetic flux in the solar atmosphere. The main target of the work was to understand the role of magnetic reconnection in jet formation as a result of reconnection low down in the solar atmosphere where dynamically the energy release can have only limited consequences. The results from this study showed that the slowmode MHD shock is crucial for the transport of energy released by the reconnection to areas of the solar atmosphere where they can be of greater dynamic importance. This example also highlights the great importance shocks can play in the transport of energy from reconnection to areas in the solar atmosphere
The connection between shock physics, magnetic fields, and partially ionised plasma is of great importance in many astrophysical systems. One key area of this study has been the study of shocks in the interstellar medium (ISM). Fastmode MHD shocks, where both the gas and magnetic field are compressed by the shock, in this weakly ionised medium involve a complex coupling between the magnetic field and the neutral gas. Low ionisation of the medium means that the shock front can become a continuous transition over a large region (known as a cshock; Draine 1980); this allows radiative cooling to more efficiently cool the gas as it shocks, significantly reducing the temperature of the postshock region and, as a result, significantly reducing the disassociation of molecules resulting from the shock (Chernoff 1987). However, certain solutions based on sufficiently large upstream velocities were found that, in spite of the continuous transition of the ionised fluid, the neutral gas would shock (known as a jshock). These solutions were normally associated with weak magnetic field strengths as this meant the magnetic field could not work to smooth out the shock front (Draine 1980). Chernoff et al. (1982) was able to use steadystate shock solutions to model molecular line intensities and predict the magnetic field strength of BNKL. For a review on this subject, see, for example, Draine & McKee (1993).
Other studies focussing on how aspects of magnetic reconnection other than its associated shocks are influenced by partial ionisation of the host plasma have lead to some very interesting results. The current sheet structure for a Harristype current sheet under the influence of partial ionisation, through the onefluid, ambipolar diffusion approximation, was found to form a power law with the magnetic field scaling as B ∝ x^{1 / 3} (Brandenburg & Zweibel 1994). Arber et al. (2009) extended this study to include a guide field in the current sheet showing that the current sheet develops into a J × B = 0 with a thickness (l_{new}) that can be estimated as l_{new} = l_{old}B_{g}/B_{ex}, where B_{g} and B_{ex} are the original guide and external magnetic fields, respectively (Singh et al. 2015). A similar problem related to the current sheet thickness of the KippenhahnSchlüter prominence model under ionneutral drift was investigated by Hillier et al. (2010), finding that the thickness of the current sheet under ambipolar diffusion would become l_{new} = l_{old}B_{hor}/B_{ver}, where B_{hor} and B_{ver} are the original horizontal and vertical magnetic fields.
The reconnection process has also been shown to be modified by the influence of partially ionised and weakly ionised plasma. Sakai & Smith (2008) and Sakai & Smith (2009) studied how reconnection in the penumbra of sunspots may be influenced by the low ionisation of the plasma finding that weak flows of the neutral hydrogen in the photosphere could greatly enhance the rate of reconnection in penumbral filaments. Singh et al. (2015) extended the fractal reconnection model of Shibata & Tanuma (2001) to include the physics of collisional coupling of a partially ionised plasma in the strongly coupled, intermediately coupled, and weakly coupled regimes through application of the growth rates for the tearing instability derived by Zweibel (1989). On application to the solar atmosphere they found that the fractal reconnection process would have to cascade down through all these levels to reach the kinetic scales that offer the possibility of fast magnetic reconnection. Leake et al. (2012) and Leake et al. (2013) investigated, through multifluid simulations, the coupling and decoupling of ion and neutral fluids in a reconnecting current sheet, finding that the reconnection inflow became decoupled, but the highvelocity reconnection jet was coupled. They also found that for Lundquist numbers smaller than for those required for the plasmoid instability (e.g. Huang & Bhattacharjee 2012) that the reconnection rate would become independent of Lundquist number. Looking at the possibility for the plasmoid instability to lead to Hallmediated reconnection in various astrophysical bodies, Vekstein & Kusano (2013) found that a prime site for the plasmoid instability leading to fast, bursty reconnection would be in the weakly ionised medium of protostellar disks. It has been shown that the development of the plasmoids in a reconnection region leads to the dynamic formation of a multitude of slow shocks associated with the formation, movement, and merger of the plasmoids (Tanuma et al. 2001; Mei et al. 2012; Shibayama et al. 2015). This finding highlights how the formation of plasmoids in a partially ionised plasma is intrinsically linked to the nature of shocks in that medium.
As yet, the nature of slowmode reconnection shocks in a partially ionised plasma are yet to be understood. In this paper we look at a simple onedimensional (1D) model that captures the necessary physics of the reconnection shock system but solves for the evolution of a neutral and ionised fluid that are coupled by collisions between the species.
2. Solving the equations for a partially ionised plasma
We investigate the dynamics of a neutral fluid and a chargeneutral, fully collisionally coupled ion electron plasma. The equations governing the neutral fluid are written as The equations governing the ionised fluid are written as (5)
We take both fluids to be ideal gases following the relations P_{n} = ρ_{n}R_{g}T_{n} and P_{p} = 2ρ_{p}R_{g}T_{p}, respectively. In this formulation α_{c}ρ_{α} = ν_{αβ}, which is the collision frequency of species α on species β where the subscripts α and β denote either n or p. We formulated these equations with extensive reference to previous numerical codes that were used to investigate partially ionised plasma (see Appendix A for more details.)
The key difference between this and single fluid models, i.e. those that approximate the ionneutral coupling through a modified Ohm’s law (e.g. Braginskii 1965), is that we solve each fluid separately, and couple them through the collisional coupling terms. These terms can be found as the source terms on the RHS of Eqs. (2), (3), (6), and (7). The code used to study this problem is the (PIP) code, which we developed to study the influence of partial ionisation on the dynamics of magnetised fluids. A basic description of the code, the full equations it can solve and tests of the collisional coupling terms are presented in Appendix A. For this study, as we are interested in a shock problem, we use an HLLD scheme (Miyoshi & Kusano 2005) because this was found to be stable down to very low ionisation fractions and plasma β values.
3. The model under consideration
In this study, we extend the model of slowshocks formed as a result of magnetic reconnection proposed by Petschek (1964). We examine the changes that occur when we look at plasmas of differing neutral fraction. The equations are normalised such that the density ρ_{tot} = 1, the collision frequency as determined by the bulk fluid density is ν = α_{c}(T_{0})ρ_{tot} = 1 = τ^{1} and the bulk Alfvén velocity the normalisation of the magnetic field is used such that . The value of α_{c}(T_{0}) is normalised to one, but α_{c}(T_{n},T_{p}) depends on the temperature as follows: (11)Physically this means that a bulk Alfvén wave approximately becomes coupled to both fluids after travelling through a unit length as defined by L_{norm} = V_{A}/ν. This normalisation means that the plasma β is defined as .
The initial conditions in normalised units are as follows: where ξ_{n} and ξ_{i} are the neutral and ion fractions, respectively. We are assuming that the two fluids are in thermal equilibrium at the start of the calculation. It should be noted that this angle of magnetic field and plasma β values used are known to produce a slowmode, switchoff shock and not a purely hydrodynamic shock (Takasao et al. 2015). The key parameters that we investigate are the ionisation fraction and plasma beta (ratio of gas to magnetic pressure) of the system, and so that these parameters are varied between simulations, the other parameters are kept the same through all the simulations. We take the adiabatic index of γ = 5 / 3. The 1D model under consideration here has been shown to give a good approximation of the dynamics of more complex reconnection simulations (Takasao et al. 2013).
As a result of the symmetry around x = 0 of the situation we are studying, we can take only onehalf the computational domain and set the boundary at x = 0 as a reflective boundary formulated in such a way that the magnetic field can penetrate the boundary. The boundary at x = x_{max} is an open boundary. As this is a 1.5D simulation the derivatives in the y direction are taken to be 0 and we do not include any velocity or magnetic field component in the z direction. The spatial resolution in the x direction of the simulations is Δx = 0.1ξ_{i0}/ξ_{i}, where ξ_{i0} = 0.1, unless otherwise stated.
3.1. The reference MHD solution
Though the MHD solution of this system has been well investigated, to provide sufficient context for the results we present a reference MHD solution. This is obtained by solving the system described above using Eqs. (5) to (10), where the source terms on the RHS of these equations are neglected (i.e. we use the ideal MHD equations and we examine only the ions).
Before we introduce the results, there is one key point that should be stressed. The equations used here are the ideal MHD equations, and these equations have the special property of having no inherent scale of reference. This is different from the case where gravity is involved, which naturally has the pressurescale height as its intrinsic length scale, or, as is relevant for this study, when partially ionised plasmas are considered as they have timescales associated with their collisional coupling. Therefore, when looking at an ideal MHD system, it is no surprise that it displays selfsimilarity, and as such we can present only one snapshot of the system here and know that it is representative of the system at all times. This does not apply to all MHD systems, but is a common feature of expanding shock systems (e.g. the hydrodynamic SedovTaylor point explosion Sedov 1959).
Fig. 1 Spatial distribution in the x direction of B_{y}a); gas pressure b); density c); v_{x}d); v_{y}e); and temperature f) for the ideal MHD simulation. There is a fastmode rarefaction wave at approximately x/t = 1 that drives a flow of material towards the shock front at x/t = 0.175. 

Open with DEXTER 
Figure 1 gives a snapshot of the distributions (going clockwise from the top left) of the B_{y} magnetic field (the B_{x} magnetic field is always constant to satisfy ∇·B = 0 and B_{x} = 0.3 in this case), gas pressure, density, temperature, v_{y} and v_{x}. The values at x/t = 1.3 are those of the initial conditions of the simulation. In this system, the relaxation of the magnetic field leads to the formation of two nonlinear waves: a fastmode rarefaction wave and a slowmode shock wave. The fastmode rarefaction wave is responsible for driving the inflow of material towards the shock front. It is characterised by a linear profile in the transition between the pre and postwave values for B_{y}, pressure, density, v_{x} and v_{y}.
Post shock there is a highvelocity jet in the y direction; this can be seen as analogous to a reconnection jet. The postshock velocity becomes v_{x} = 0 in the simulation frame, as does the postshock ycomponent of the magnetic field (B_{y}). The shock jump conditions for this particular MHD shock (in the shock reference frame) are (see, for example, Goedbloed & Poedts 2004) written as
where the subscripts U and D denote the upstream and downstream quantities. We also need to supplement this with a condition that the entropy must increase across the discontinuity, i.e. (24)Equation (19) tells us that the ratio of the two densities is given by the ratio of v_{xU}/v_{xD}, which in this case is approximately 1.9. By combining Eqs. (19), (20), and (22) we can derive the relation , where is the upstream Alfvén velocity (this was checked to be true in the simulations to an accuracy of 0.0001 per cent). As the upstream velocity in the y direction is likely to be significantly smaller than the Alfvén velocity, it can be said that approximately the postshock jet is travelling at the Alfvén velocity. Then by combining Eqs. (19) and (21) we can show that , which implies that the gas pressure has to increase to match both the dynamic and magnetic pressure drop, making it possible for pressure increases that are much larger than those possible for purely hydrodynamic shocks, as both of these terms can be significantly larger than p_{U} in this model for low plasma β.
Fig. 2 Spatial distribution in the x direction of B_{y}a); gas pressure b); density c); v_{x}d); v_{y}e); and temperature f) for the neutral (red) and ionised (blue) fluids at time t = 1τ. The green line indicates the total (pressure and density) or the difference (x and y velocities and temperature) for the two fluids. The dashed black line shows the reference MHD solution. The ionised plasma has taken on characteristics that are similar to the ideal MHD solution (but with faster wave speeds); the neutrals, however, are undergoing a violent expansion similar to that of a 1D point explosion. 

Open with DEXTER 
Fig. 3 Panel a): temporal evolution of the temperature; panel b): magnitude of the heating terms relating to collisional coupling at x = 0 for a timescale that follows the transition from the initiation to the weak coupling regimes (see Sect. 4.2). The normalisation of the power is . Because the temperatures cross at t ~ 3.5τ the thermal damping term become 0 at this point then reverse sign (i.e. the ionised fluid is losing heat to the neutral fluid before this time, but gaining heat from the neutral fluid after this point). 

Open with DEXTER 
4. Temporal evolution
In this section we take one specific set of parameters and study in detail the evolution of the formation of the shock and rarefaction waves through time, focussing on times that are smaller than the coupling time to those that are significantly longer. For this part of the study, we take ξ_{n} = 0.9 and β = 0.3.
4.1. Initiation
Fig. 4 Spatial distribution in the x direction of B_{y}a); gas pressure b); density c); v_{x}d); v_{y}e); and temperature f) for the neutral (red) and ionised (blue) fluids at time t = 10τ. The green line indicates the total (pressure and density) or the difference (x and y velocities, and temperature) for the two fluids. The dashed black line indicates the reference MHD solution. 

Open with DEXTER 
The initialisation of the dynamics is dominated by the evolution of the magnetised fluid; this is simply because of the initial conditions (i.e., only the magnetic field is out of equilibrium). When the system is initialised, as with the ideal MHD simulation in Sect. 3.1, a fastmode rarefaction wave and a slowmode shock form in the ionised fluid. For these initial conditions, the bulk Alfvén velocity of the fluid is normalised to 1, therefore the Alfven velocity with respect to the ionised fluid is larger, in this case . This results in a fastmode rarefaction wave travelling only in the ionised fluid away from the position of the original discontinuity in the magnetic field at approximately the ion Alfvén speed. This wave is accelerating the ionised fluid towards a slowmode MHD shock. Downstream of this shock, there is a jet of the ionised fluid that is travelling at a velocity of a few times the bulk Alfvén speed. However, as there is no force on the neutrals that can accelerate them in that direction until the coupling has taken effect, their velocity in the y direction is 0. This large velocity of the ionised fluid jet in the y direction, in conjunction with the initially 0 velocity of the neutral fluid in this direction, results in a large velocity difference v_{D} in the jet region.
The large v_{D} in the jet region results in the nonlinear terms of the coupling (those of order ) becoming much larger than those of order v_{D}. Therefore, the coupling of the fluids comes about via the nonlinear coupling terms in the energy equations. The lefthand panel of Fig. 3 gives the temperature as a function of time at x = 0 for both the ionised and neutral fluids. Initially the temperature of the ionised fluid increases drastically followed more slowly by the neutral fluid after about 3.5τ (where τ = ν^{1}) the two fluids have approximately the same temperature. The right panel shows the size of the heating (cooling) terms associated with frictional heating given by α_{c}(T_{n},T_{p})ρ_{n}ρ_{p}(v_{n} − v_{p})^{2} (solid line) and thermal damping term that drives thermal equilibrium given by 3α_{c}(T_{n},T_{p})ρ_{n}ρ_{p}(p_{n}/ρ_{n} − p_{p}/ 2ρ_{p}) (dashed line) and the two terms for the work performed on the fluid by the drift velocity given by  α_{c}(T_{n},T_{p})ρ_{n}ρ_{p}v_{n}·(v_{n} − v_{p})  (blue line) and  α_{c}(T_{n},T_{p})ρ_{n}ρ_{p}v_{p}·(v_{n} − v_{p})  (red line), respectively. The thermal damping term is equal in magnitude but opposite in sign for each fluid. This highlights how the frictional heating term, which represents a nonreversible increase in entropy, forms a significant proportion of the heating in the early stages.
The neutral fluid presents a very interesting response to its initial heating. The localised heating of this fluid results in an increase in pressure around x = 0, but it remains at its initial value elsewhere. This means we have a highly localised, hightemperature region, which, other than the external plasma not being cold, results in a SedovTaylor like explosive expansion of the inner hightemperature neutral layer. This can be seen from the neutral expansion, driven by a blast wave, shown in panel (d) of Fig. 2.
4.2. Weak coupling
Figure 4 shows the system early in its evolution at a time (10τ) that is an order of magnitude greater than that of the previous subsection. As the coupling between the ions and neutrals is beginning to take effect, we can see that structures that are present in both fluids are beginning to form. There are four components of the system: the rarefaction wave, neutral shock, ion shock, and postshock region. All the regions in which flows are present display velocity drifts and these drifts are of approximately the same magnitude as the flows themselves.
By this time, both an inflow in the preshock region and postshock jet are developing in the neutral fluid as a result of the coupling between the two fluids. There is a shock in the neutral fluid that is propagating away from the origin in advance of the ionised plasma. This is the continuation of the explosion in the neutral fluid described earlier. However, as the post neutralshock region is a prime site for frictional heating because of the large differences in velocity, it results in energy that is constantly added to this shock wave resulting in the peak in pressure that is some distance behind the shock front. The expansion has resulted in a density depletion at the origin.
4.3. Intermediate coupling
Fig. 5 Spatial distribution in the x direction of B_{y}a); gas pressure b); density c); v_{x}d); v_{y}e); and temperature f) for the neutral (red) and ionised (blue) fluids at time t = 100τ. The green line indicates the total (pressure and density) or the difference (x and y velocities, and temperature) for the two fluids. The dashed black line indicates the reference MHD solution. 

Open with DEXTER 
Figure 5 looks at the system at a time (100τ) that is approximately an order of magnitude later than that of the previous subsection, and as such the dynamics of the two fluids show greater coupling. The high pressure neutral region has sufficiently expanded so that the pressure at the front of this expansion has dropped to match that of the ambient pressure, and as such the shock in the neutrals has disappeared and is replaced with a smooth transition between the inflow region of neutral and ionised fluid and the jet region. It is in this smooth transition region, which represents a cshock, that the high levels of drift velocity are found (~0.15V_{A}). The two waves that we expect to form in this system, the slowmode shock and fastmode rarefaction wave, are visible, but are yet to form as completely distinct elements.
It can be seen that in the region of the cshock, the neutral fluid is undergoing the velocity transition before the ionised fluid. It is worth noting that this is a clear difference from the cshocks studied that are relevant to the ISM. This leads to the preheating in the shock front relative to the reference MHD solution.
4.4. Quasiselfsimilar state
Once the system has sufficiently evolved then it reaches a state that is approaching that of selfsimilarity. Figure 6 gives an example of this later stage at a time (2000τ) that is a factor of 20 greater than that used in the previous subsection. There are now two distinct waves in the system, the fastmode rarefaction wave and the slowmode shock. The pre and postshock values of the physical quantities are basically constant and large values of the drift velocity are only found in the shock front. As the drift velocity in the x direction is positive, we can tell that the neutral fluid undergoes the shock transition before the ionised fluid, but for the postshock jet the ionised fluid is accelerated first and then the neutral fluid is dragged with it. The nonlinear coupling terms related to frictional heating are most important in this region.
Fig. 6 Spatial distribution in the x direction of B_{y}a); gas pressure b); density c); v_{x}d); v_{y}e); and temperature f) for the neutral (red) and ionised (blue) fluids at time t = 2000τ. The green line indicates the total (pressure and density) or the difference (x and y velocities, and temperature) for the two fluids. The dashed black line indicates the reference MHD solution. The solution for the partially ionised case has reverted to that of the ideal MHD case outside of the wave fronts. 

Open with DEXTER 
It is worth considering what conditions are necessary to produce this steadystate shock. Eventually, the rarefaction wave and shock front decouple, and after this time the density, momentum, energy, and magnetic flux that leaves the rarefaction wave and as such composes the upstream conditions for the shock becomes almost constant in time (at least rapidly tends to such a state). This results in the shock jump conditions that need to be satisfied at each timestep becoming the same, which results, when travelling in the shock reference frame, in a steadystate shock. In the following, we provide an analysis of why this quasiselfsimilar solution develops and what determines the end state of this system.
Figure 7 shows the results from a simulation where β = 0.3 and ξ_{n} = 0.99. This calculation has been run for a longer time (using Δx = 10) specifically to analyse the evolution towards selfsimilarity. The xaxis of this figure is the length scale in the x direction divided by the time that the snapshot of the simulation was taken. In this way, a feature travelling at a constant velocity always appears at the same point in the figure at all times. Here we can see that the position of centroid of rarefaction wave and the start of the postshock region evolve in a selfsimilar fashion (i.e. these points travel at a constant velocity). The relative thickness of these two layers, however, does change; a detailed investigation of this is presented below.
Fig. 7 Quasiselfsimilar evolution of the shockrarefaction wave system. To highlight how the system is tending towards selfsimilarity, the xaxis is normalised by the time t of the snapshot of the simulation, therefore shows the velocity at which the waves are propagating in normalised units. As the time is increased, the system tends towards a selfsimilar state. The reference time is taken at t_{0} = 156 000τ and the distributions for t_{0}, t_{0}/ 2, t_{0}/ 4, t_{0}/ 8, and t_{0}/ 16 are plotted. The reference ideal MHD solution is indicated by the black dashed line. 

Open with DEXTER 
4.4.1. The effective thinning of the rarefaction wave
Figure 7 shows the evolution of the rarefaction wave as the system heads towards a final state. In the velocity reference frame used, the thickness of the rarefaction wave decreases with time, tending towards the reference ideal MHD solution. This effective thinning of the wave front is a result of the diffusive processes becoming less important with time when compared with the dynamic processes of the system as it evolves. This can be simply modelled by taking the thickness of the rarefaction wave (W_{R}) to be given by the sum of component that comes from the expansion of the system (At) and diffusion of the system (Bt^{1 / 2}), i.e. (25)where t> 0. Therefore, (26)which implies that a t gets larger, the influence of the diffusive component becomes smaller meaning that (27)i.e. the system is heading towards a selfsimilar state at the rate of t^{− 1 / 2}. As a note, the physical meaning of A is the velocity difference between the entry and exit points of the rarefaction wave for a fluid parcel that is given by the ideal MHD solution. For this particular simulation, by calculating the Dopplershifted fastmode wave speeds (in the simulation frame) upstream and downstream of the wave gives V_{fU} = 1.150V_{A} and V_{fD} = 0.926V_{A}, therefore A = 0.223V_{A}.
Fig. 8 Drift velocity through the shock front for the simulation with β = 0.3 and ξ_{i} = 0.01. The xaxis is shifted by x_{s}, where x_{s} is taken as the point in the shock where the drift velocity is at its peak, so that both distributions are aligned. It is clear that the distribution of the drift velocities is the same at both times indicating that the shock front has reached a steady state. The reference time is taken at t_{0} = 120 000τ. 

Open with DEXTER 
Another way of thinking about this relates to the drift velocities and relative frequencies of the momentum coupling as the rarefaction wave expands. If we think about a packet of the neutral fluid and its acceleration by the rarefaction wave towards the shock front entering the wave at time t_{0} and exiting at time t_{1}. The amount of momentum transferred to the neutral fluid by collisions with the ionised fluid (ρ_{n}v_{xn})_{C} must be constant in time to give the same velocity, i.e.(28)Therefore, as the time taken for the packet of neutral fluid (t = t_{1} − t_{0}) to cross the wave front gets larger, then the average velocity drift that the fluid packet experiences (⟨ v_{xn} − v_{xp} ⟩) tends to zero. As this happens, the nonlinear term involved with the dissipation of the kinetic energy associated with the drift velocity through frictional heating must also change as follows: (29)With this approach, we can see that two things are happening: 1) the system is tending towards a state of perfect coupling where the ratio of the dynamic timescale (τ_{D}) to the collisional coupling timescale τ_{D}ν → ∞; and 2) that the frictional heating term becomes 0, which means that the nonlinear terms cannot change the state of the upstream shock conditions.
Fig. 9 Spatial distributions of v_{xn}a) and v_{yn}b) for β = 0.1 and ξ_{i} = 0.1, 0.01, 0.001, 0.0001, and 0.00001 denoted by the green, purple, pink, red, and turquoise lines, respectively. 

Open with DEXTER 
4.4.2. Steadystate shock layer
Figure 8 shows the drift velocity distribution across the shock front for the simulation where β = 0.3 and ξ_{i} = 0.01. The central position of the shock has been normalised as the position of the peak drift velocity. The distribution is taken at two different times t = 60 000τ and t = 120 000τ, but the distribution of the velocity drift across the shock is the same at these two times. After a sufficient period of time the shock reaches a steady state.
It is worthwhile here to discuss the difference between the shock structure found for these calculations and those investigated for the ISM. Because this study deals with a slowmode shock and not a fastmode shock as well as plasma β values that are less than 1, it is very natural that there are differences in the two approaches. In the fastmode shock case, the extra magnetic pressure that comes from the compression of the magnetic field results in the expansion of the ionised fluid layer, which couples to the neutral fluid and expands that layer. In this case, as the slowmode shock has a decrease in magnetic pressure that is partly balanced by an increase in gas pressure, it is the neutral fluid that is overpressured by the shock and this region pushes forwards.
5. Parameter dependence of shock structure
In the previous section, we mainly focussed on the temporal evolution of one particular example. Now we focus on the structure of the shock region once the system has reached the quasiselfsimilar state as described in Sect. 4.4 and investigate the influence of some of the parameters of the system on the structure of the shock. Keeping all other parameters the same, we investigate the twodimensional parameter space of ion fraction ξ_{i} and plasma β.
5.1. Ionisation fraction dependence
Figure 9 gives the velocity distribution in the x (panel a) and y (panel b) directions for the case where β = 0.1 and ξ_{i} varies between 0.00001 and 0.1. The time that the snapshot is taken from each simulation is given by t_{snap} = 150τ/ξ_{i}. The axis of this figure is given as xξ_{i}. This is because the collision frequency of the neutrals to the ions is proportional to the ionisation fraction and so it could be simply expected that the thickness of the shock region (T_{shock}) would roughly scale as T_{shock} = V_{shock}/α_{c}ρ_{tot}ξ_{i} (e.g. Draine & McKee 1993). It is clear that this is not the case, but interestingly, for the v_{x} velocity the exit point of the shock position is a point of the system that evolves selfsimilarly (see Fig. 7); therefore, the time scaling used for the snapshot multiplied by ξ_{i} relates to the selfsimilar evolution of the system and that appears to be independent of ionisation fraction and hence aligns in this figure. This is not the case for the v_{y} velocity as both the entrance and exit points for the fluid from the shock move further outwards as the ionisation fraction decreases.
It is also clear in Fig. 9 that for the ionisation fraction of ξ_{i} = 0.1 there is an overshoot in the neutral velocity. Though such features often appear in shock simulations as a result of numerical effects, this particular feature is very well resolved (approximately 100 grid points) and as such can be taken as a genuine feature of this particular shock.
5.1.1. Structure of the drift velocity across the shock
Figure 10 gives the distribution of the drift velocity across the shock front for the range of ionisation fractions under study. For large values of ξ_{i} the drift velocity profile for the x velocities is approximately a Gaussian distribution, but as ξ_{i} gets smaller, the distribution changes and it develops a strong skew with the transition at the entry to the shock region becoming associated with a significantly sharper transition than the exit from the shock region. The ydirection velocities with the drift velocity felt by a fluid packet steadily increasing as it moves through the shock before suddenly returning to 0 for the case of ξ_{i} = 0.1. However, the skew of this distribution also becomes progressively more and more positive as the ionisation fraction is decreased.
Fig. 10 Spatial distributions of v_{xn} − v_{xp}a) and v_{yn} − v_{yp}; b) for ξ_{i} = 0.1, 0.01, 0.001, 0.0001, and 0.00001. 

Open with DEXTER 
5.1.2. Thickness of the shock region
The thickness of the shock region as a function of the ionisation fraction is given in Fig. 11. This thickness is calculated as the extent over which the drift velocity is creating momentum coupling in the x direction across the shock front (it is clear from Fig. 10 that the values calculated from using the velocity drift in the y direction would result in a different scaling). The scaling is given for β = 0.3, 0.1, 0.01, and 0.001.
The general trend shown in this figure is that the thickness of the layer decreases monotonically with increase in the ionisation fraction. However, as can be expected from the results in Figs. 9 and 10, the relation is not linear, but has a dependence that is approximately ~, as shown by the dashed line in Fig. 11
Fig. 11 Thickness of the shock region as a function of ξ_{i} for β = 0.3 (green), 0.1 (blue), 0.01 (pink), and 0.001 (red). The black dashed line shows a power law as a function of ξ_{i} with exponent − 1.2 

Open with DEXTER 
5.2. Plasma β dependence
The previous subsection dealt with the influence of the ionisation fraction; here we now look at the dependence on the shock front of the plasma β. Figure 12 is the same as Fig. 9 but the different curves represent different plasma β values instead of different ξ_{i} values. We are looking at the case of ξ_{i} = 0.01. The key point of this figure is that as plasma β gets smaller, the system undergoes a transition from producing cshocks to that of producing jshocks where there is the jump in physical variables associated with a shock (in the x velocity and also the pressure and density, but not the y velocity). The physics behind this are investigated in a following subsection.
It can be expected that the overshoot in the x velocity (something that is resolved by approximately 100 grid points) is a real feature of the system. In fact, we can estimate the downstream velocity of this shock by just performing a hydrodynamic analysis using the shock jump conditions on the shock, i.e. assuming that as this is a discontinuity in which the neutrals are not influenced by the ions and so the shock jump conditions should be just those of the hydrodynamic RankineHugoniot relation. For an ideal gas with adiabatic index γ = 5 / 3, the ratio of the upstream and downstream velocities in the shock frame is u_{U}/u_{D} → 4 as the shock gets progressively stronger. For the case of β = 0.01, we have a shock velocity of ~0.15V_{A} and an inflow velocity in the simulation frame of ~− 0.2V_{A}, which would give an upper estimate for the postshock velocity in the simulation frame of ~0.06V_{A}, i.e. the velocity undergoes a positive overshoot above the v_{x} = 0 level that is expected for the final state after the fluids have recoupled. The actual postshock velocities shown in Fig. 12 are ~0.02V_{A} so they are within the predicted limit. Once the shock has passed, the overshoot in the velocity decays exponentially, through collisional coupling over approximately one coupling length scale for neutrals to couple to the ions towards the v_{x} = 0 value.
It is also worth pointing out that there is an inflow dependence on β, which is driven by the rarefaction wave. The inflow depends on the physics of expansion of the system and the position of the shock front depends on the plasma β. It is likely that part of the change in position of shock front (i.e. propagation speed of the shock front) is a result of the Doppler shifting of the shock by the inflow. Compression of the shock also gives an explanation for what is happening.
Fig. 12 Spatial distributions of v_{xn}a) and v_{yn}b) for ξ = 0.01 and β = 0.1, 0.01 and 0.001 denoted by the green red and blue lines, respectively. 

Open with DEXTER 
Figure 13 focusses on the distributions for both the neutral and ionised fluid around the shock region for the case where β = 0.01 and ξ_{i} = 0.01. A sharp transition is clear in v_{x}, P, and ρ for both the neutral and ionised fluid. This is a distinct difference from the jshock solution as presented in Fig. 3 of Draine & McKee (1993), where the shock is only present in the neutral fluid because for this case when there is a jshock in the neutral fluid, it is also present in the ionised fluid. The pressure evolution through the shock front is rather interesting. Though the neutral fluid shocks, where the kinetic energy transfers to thermal energy, most of the energy conversion goes from magnetic energy to thermal energy and kinetic energy (for the velocity in the y direction). Therefore, even after the shock, the pressure steadily increases. The density overshoots the postshock value as the neutral fluid shocks, but because of the diverging neutral velocity field the density decays to the postshock value.
Fig. 13 Spatial distribution in the x direction around the shock front and the coupling region of a) B_{y}; b) gas pressure; c) density; d)v_{x}; e) v_{y}; and f) the temperature for the neutral (red) and ionised (blue) fluids at time t = 15 000τ where β = 0.01 and ξ_{i} = 0.01. The green line indicates the total (pressure and density) or the difference (x and y velocities and temperature) for the two fluids. The plasma pressure and density are increased by a factor of 50 to make their distributions clearly visible. The xaxis has been shifted to set the origin at the shock front. 

Open with DEXTER 
5.2.1. Structure of the drift velocity across the shock
Figure 14 shows the drift velocity distributions for both the x and y velocities for different plasma β values. Once the shock has transitioned from a cshock to a jshock, the structure of the drift velocity is characterised by the sharp jump of the shock and a monotonic drop towards the 0.
Fig. 14 Spatial distributions of v_{xn} − v_{xp}a) and v_{yn} − v_{yp}; b) for β = 0.1, 0.01 and 0.001 denoted by the green red and blue lines, respectively. 

Open with DEXTER 
5.2.2. Thickness of the shock region
Figure 15 gives the thickness of the shock region as a function of plasma β (right) and the inflow Mach number of the neutrals in the shock frame (left). Though there is a monotonic increase with plasma β, it does not become a clear power law. To understand what happens when in the shock front when β changes, it is worth looking at the upstream Mach number (calculated in the shock frame) of the neutral fluid. Once the inflow speed exceeds the sound speed, a shock is formed (see discussion in next paragraph and Fig. 16); this limits the amount that it can be further compressed, so the nature of the coupling (where the temperature jump increases the collision frequency) then becomes important for determining how thick the layer becomes.
Fig. 15 Thickness (in units of V_{A}/ν_{T}) of the shock for the simulations with ξ_{i} = 0.01 as a function of β (left) and , the hydrodynamic Mach number of the neutral fluid in the inflow region (right). The dependence of the thickness of the shock as a function of the plasma β can be explained by the Mach number of the inflow region to which this relates to the compression at the shock front. 

Open with DEXTER 
Fig. 16 Spatial distribution of  v_{xn}/Cs_{n}  and  v_{xp}/V_{A}  for a cshock (β = 0.3 and ξ_{i} = 0.01) in panel a), a weak jshock (β = 0.1 and ξ_{i} = 0.01) in panel b), and a jshock (β = 0.01 and ξ_{i} = 0.01) in panel c) around the shock front (calculated in the shock frame). The xaxis is shifted by x_{s}, where x_{s} is taken as the point in the shock where the drift velocity is at its peak, so that all distributions are aligned. 

Open with DEXTER 
Figure 16 shows the Mach number in the shock frame (calculated using v_{xn}) of the neutral fluid and the Alfvénic Mach number in the shock frame of the ionised fluid (calculated using v_{xp}) for β = 0.3, 0.1 and 0.01. This choice of plasma β gives a subsonic cshock, a transonic weak jshock, and a strong jshock. The xaxis is shifted by x_{s}, where x_{s} is taken as the point in the shock where the drift velocity is at its peak, so that all distributions are aligned. Here it is clear that for a jshock to occur there has to be a supersonic inflow velocity, otherwise a cshock forms. For the case where the inflow velocity is only weakly supersonic, the shock is weak and the coupling across the shock front is very similar to that of the cshock. The Mach number drops to less than 0.3 for all cases shown at postshock . At all times the ionised fluid is subAlfvénic.
5.2.3. Fine structure of the jshock
In the shock frame, the ionised fluid is travelling towards the shock front, in this case at a velocity of ~0.4V_{A}, so it can be expected that for a finite distance D the ionised fluid could completely decouple from the neutral fluid around the shock front. This distance can be estimated as D = V_{ion}/α_{c}ρ_{n} and for the simulation of β = 0.01 and ξ_{n} = 0.99 implies that D is of order 0.1. Figure 17 is a zoom in of the shock region of Fig. 13. Unlike Fig. 13, for this calculation, Δx = 0.01, and the units of the xaxis of the figure have not been rescaled by ξ_{i}.
For the neutral fluid the shock transition happens over a few grid points determined by the numerical dissipation of the scheme, and the ionised fluid has a much wider transition. The width of this transition is approximately the same length scale D as estimated in the previous paragraph.
The finite width of the ionised fluid shock may have great importance for understanding how other diffusive and dissipative effects play a role in heating in the shock front. If the thickness of the shock for the ionised fluid is determined by coupling processes then the Hall diffusion, magnetic diffusion, and viscosity would all be in regimes where their influence is small. As the ionised fluid maintains a finite thickness, it is likely that the neutral viscosity associated with the shock compression would become the other major dissipative process in the system.
Fig. 17 Spatial distribution in the x direction around the shock front of a) B_{y}; b) the gas pressure, c) the density; d)v_{x}; e) v_{y}; and f) the temperature for the neutral (red) and ionised (blue) fluids at time t = 15 000τ where β = 0.01 and ξ_{i} = 0.01. The green line indicates the total (pressure and density) or the difference (x and y velocities and temperature) for the two fluids. The plasma pressure and density are increased by a factor of 50 to make their distributions clearly visible. 

Open with DEXTER 
Fig. 18 Relation between the total energy that goes to heating the fluids through frictional heating as a package of fluid passes through the shock front for panel a)ξ_{i}; and panel b) β for ξ_{i} = 0.01. 

Open with DEXTER 
5.3. Frictional heating across the shock front
Figure 18 indicates the total heating through frictional heating across the shock front as a function of ξ_{i} (panel a) and as a function of β for ξ_{i} = 0.01 (panel b). To give a context of the values shown, we note that the normalisation is the reference magnetic energy of . Therefore, the equivalent of a maximum of approximately 2 per cent of the reference magnetic energy is converted to heat in the shock front through the frictional heating process.
There are trends present in both the dependence on ξ_{i} and plasma β. As the ionisation fraction decreases, the width of the shock front increases and, based on the arguments for the heating in the rarefaction wave presented in Sect. 4.4.1, we can expect that this leads to a decrease in the size of the nonlinear term resulting in a decrease in the overall heating. For the plasma β, the heating peaks at β = 0.1 (which is approximately where ). Because of the large temperature increase that results from the low β shocks, the collisional frequency increases and this results in the reduction of the nonlinear terms across the shock front.
6. Discussion
In this study we investigated the 1D partially ionised MHD Riemann problem initiated by a magnetic slingshot in which the fully ionised ideal MHD version forms two waves: a fastmode rarefaction wave and a slowmode shock wave. The partially ionised system evolves through four stages: the initiation, where most of the dynamics are in the ionised fluid and the coupling is mainly through the nonlinear energy terms; weak coupling, where the momentum coupling is beginning and the neutral dynamics are characterised by an explosive outflow; intermediate coupling, where the rarefaction wave and shock are forming as distinct entities in the system; and finally reaching a quasiselfsimilar state. In the quasiselfsimilar stage, the shock forms a steady state but the rarefaction wave undergoes a diffusive evolution as it tends towards its final state. A highspeed (approximately the Alfvén velocity) jet forms behind the shock. This is driven in the ionised fluid by the Lorentz force and in the neutral fluid by collisional coupling.
Two types of shock were found, a cshock solution for low Mach number inflow and a jshock solution when the inflow velocity becomes supersonic. These are characterised by smooth transitions in the x velocity, pressure, and density for both fluids in the cshock and sharp transitions in both fluids for the jshock. Though in the model used the inflow is driven by a fastmode rarefaction wave, which, as it is a compressional wave, gives a smaller inflow velocity when it can expand in threedimensions (3D) than in 1D, so the dependence of the transition between the c and jshocks on the plasma β may not be so applicable outside the model. The dependence on the Mach number, however, is a universal rule because it relates exactly to the conditions of the shock jump, which is an inherently 1D problem.
One important point to discuss is the similarities and differences found in this study from those of shocks in the ISM. Though there are some large differences in which physical processes are under focus for each study, there is a great deal of comparison that can be made. A key difference that should be highlighted is that because of the different geometries of the magnetic field under consideration, different shocks are formed. Those under consideration in the ISM are fastmode MHD shocks, where both the magnetic field strength and gas pressure increase downstream of the shock front, but we have been investigating a slowmode MHD shock which through magnetic tension has a very large component of the downstream velocity parallel to the shock front. The existence in both c and jshock solutions for both studies is one area of similarity, though the conditions for the appearance of the jshock solution differ between the two studies. For this study, jshocks appeared based on the hydrodynamic Mach number of the upstream flow: once this became supersonic, shocks formed. Charge exchange and ionisation/recombination were shown to strongly influence the shock structure in the ISM, and it would be no surprise if the shock studied here was the same. For the influence of charge exchange it is likely that this process would act as an effective increase in the collisional coupling (e.g. Terradas et al. 2015). It can be expected that, as the ionisation fraction is so important for determining the size of the coupling region around the shock, ionisation associated with the shock reduces the size of the shock region.
One application of this research is the investigation of shockdriven jets in the solar atmosphere. Reconnection driving shock formation low in the solar atmosphere has been proposed as a key part of cool jet formation (e.g. Shibata et al. 1982; Nakamura et al. 2012; Takasao et al. 2013). Our results suggest that owing to the changes in ionisation fraction and density with height in the solar atmosphere (Vernazza et al. 1981) that any slowmode shock that is propagating is constantly evolving as it travels through the atmosphere. Our results also suggest that though limited in their spatial extent, observationally shock fronts, offer the greatest chance of observing ionneutral drift in the solar atmosphere as shocks pass through the partially ionised regions of the atmosphere.
One interesting thought relates to the large gradients in v_{yn} with x that result from the magnetic field driving the plasma jet. Though this strong shear would only be experienced by an individual fluid packet for a finite amount of time as it passes through the shock front, there is the potential that it could have some interesting dynamic consequences. There is the possibility that it could drive the formation of velocity shear driven instabilities in the neutral fluid, which would be
associated with turbulence in the neutral fluid in the reconnection jet and potentially thickening the velocity transition region associated with the shock. It would be an interesting topic to investigate further.
Acknowledgments
A.H. is supported by his STFC Ernest Rutherford Fellowship grant number ST/L00397X/1. S.T. acknowledges support by the Research Fellowship of the Japan Society for the Promotion of Science (JSPS). The authors would like to thank Drs. Hiroaki Isobe, Kazunari Shibata, Alexander Russell, & David Chernoff for their insightful observations regarding this research. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk. This equipment was funded by a BIS National Einfrastructure capital grant ST/K00042X/1, STFC capital grant ST/K00087X/1, DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National EInfrastructure.
References
 Arber, T. D., Botha, G. J. J., & Brady, C. S. 2009, ApJ, 705, 1183 [NASA ADS] [CrossRef] [Google Scholar]
 Braginskii, S. I. 1965, Rev. Plasma Phys., 1, 205 [NASA ADS] [Google Scholar]
 Brandenburg, A., & Zweibel, E. G. 1994, ApJ, 427, L91 [NASA ADS] [CrossRef] [Google Scholar]
 Chernoff, D. F. 1987, ApJ, 312, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Chernoff, D. F., McKee, C. F., & Hollenbach, D. J. 1982, ApJ, 259, L97 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T. 1980, ApJ, 241, 1021 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T., & McKee, C. F. 1993, ARA&A, 31, 373 [NASA ADS] [CrossRef] [Google Scholar]
 Goedbloed, J. P. H., & Poedts, S. 2004, Principles of Magnetohydrodynamics, eds. J. P. H. Goedbloed, & S. Poedts (UK: Cambridge University Press), 5 [Google Scholar]
 Harten, A., Lax, P. D., & van Leer, B. 1983, SIAM Rev., 25, 35 [CrossRef] [MathSciNet] [Google Scholar]
 Hillier, A., Shibata, K., & Isobe, H. 2010, PASJ, 62, 1231 [NASA ADS] [Google Scholar]
 Huang, Y.M., & Bhattacharjee, A. 2012, Phys. Rev. Lett., 109, 265002 [NASA ADS] [CrossRef] [Google Scholar]
 Inoue, T., & Inutsuka, S.I. 2008, ApJ, 687, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Jefferies, J. T. 1968, Spectral Line Formation (Waltham, Mass. USA: Blaisdell) [Google Scholar]
 Leake, J. E., Lukin, V. S., Linton, M. G., & Meier, E. T. 2012, ApJ, 760, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Leake, J. E., Lukin, V. S., & Linton, M. G. 2013, Phys. Plasmas, 20, 061202 [NASA ADS] [CrossRef] [Google Scholar]
 Longcope, D. W., & Bradshaw, S. J. 2010, ApJ, 718, 1491 [NASA ADS] [CrossRef] [Google Scholar]
 Meier, E. T., & Shumlak, U. 2012, Phys. Plasmas, 19, 072508 [NASA ADS] [CrossRef] [Google Scholar]
 Mei, Z., Shen, C., Wu, N., et al. 2012, MNRAS, 425, 2824 [NASA ADS] [CrossRef] [Google Scholar]
 Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Nakamura, N., Shibata, K., & Isobe, H. 2012, ApJ, 761, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Petschek, H. E. 1964, NASA Sp. Publ., 50, 425 [NASA ADS] [Google Scholar]
 Petschek, H. E., & Thorne, R. M. 1967, ApJ, 147, 1157 [NASA ADS] [CrossRef] [Google Scholar]
 Priest, E. R., & Forbes, T. G. 1986, J. Geophys. Res., 91, 5579 [NASA ADS] [CrossRef] [Google Scholar]
 Sakai, J. I., & Smith, P. D. 2008, ApJ, 687, L127 [NASA ADS] [CrossRef] [Google Scholar]
 Sakai, J. I., & Smith, P. D. 2009, ApJ, 691, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (New York: Academic Press) [Google Scholar]
 Shibata, K., Nishikawa, T., Kitai, R., & Suematsu, Y. 1982, Sol. Phys., 77, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Shibayama, T., Kusano, K., Miyoshi, T., Nakabou, T., & Vekstein, G. 2015, Phys. Plasmas, 22, 100706 [NASA ADS] [CrossRef] [Google Scholar]
 Shibata, K., & Magara, T. 2011,Liv. Rev. Sol. Phys., 8 [Google Scholar]
 Shibata, K., & Tanuma, S. 2001, Earth, Planets, and Space, 53, 473 [Google Scholar]
 Singh, K. A. P., Hillier, A., Isobe, H., & Shibata, K. 2015, PASJ, 234 [Google Scholar]
 Sonnerup, B. U. Ö. 1970, J. Plasma Phys., 4, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Takasao, S., Isobe, H., & Shibata, K. 2013, PASJ, 65, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Takasao, S., Matsumoto, T., Nakamura, N., & Shibata, K. 2015, ApJ, 805, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Tanuma, S., Yokoyama, T., Kudoh, T., & Shibata, K. 2001, ApJ, 551, 312 [NASA ADS] [CrossRef] [Google Scholar]
 Terradas, J., Soler, R., Oliver, R., & Ballester, J. L. 2015, ApJ, 802, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Toro, E. F., Spruce, M., & Speares, W. 1994, Shock Waves, 4, 25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vekstein, G., & Kusano, K. 2013, MNRAS, 434, 1789 [NASA ADS] [CrossRef] [Google Scholar]
 Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635 [NASA ADS] [CrossRef] [Google Scholar]
 Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zweibel, E. G. 1989, ApJ, 340, 550 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: The (PIP) code
The (PIP) code is a 3D numerical code designed to investigate the dynamics of partially ionised plasma across a range of spatial and temporal scales that can occur in astrophysical systems. In this section we detail the full set of equations used in the code and a method to solve the collisional coupling terms to reduce the constraint they can have on the timestep of a simulation that is similar to the method presented in Inoue & Inutsuka (2008).
Appendix A.1: Equations
In this code, two sets of equations are solved separately, which are then coupled using collisions, ionisation, and recombination. The equations used are those of an MHD fluid and a hydrodynamic (HD) fluid that are joined through collisional coupling. The subscripts p and n are used to denote the variables for the MHD and HD equations, respectively. The formulation of these terms can be found in Leake et al. (2012), with Meier & Shumlak (2012) and Braginskii (1965) as supplementary material. For all the equations, the collisional, ionsation, and recombination terms are included on the RHS of the equation.
The full equations solved for the evolution of the neutral hydrogen fluid are written as
The full equations solved for the evolution of the charge neutral ionelectron plasma fluid are written as Equation (A.7) contains the Ambipolar diffusion term (η_{A}B^{2}J_{⊥}). However, this term can only be activated when single fluid simulations are performed. The thermal conduction terms are solved using the method detailed in Takasao et al. (2015). The ionisation and recombination rates are given by (Jefferies 1968) Two numerical schemes have been implemented: a fourthorder fourRungeKutta central difference scheme (Vögler et al. 2005) and an HLL scheme (Harten et al. 1983) with the HLLC scheme for the HD equations (Toro et al. 1994) and the HLLD scheme for the MHD equations (Miyoshi & Kusano 2005).
Appendix A.2: Time integration of collisional source terms
We have implemented two methods for calculating the update in the conserved variables. The first method is a simple numerical integration with time of the source terms (used in this paper because the dynamic timescales are less than the collisional timescales) and the second is a method that employs an exact solution for the integration of the equations that follow (this is appropriate when the collision timescale is much smaller than the other timescales of the system). A similar method can be found in Inoue & Inutsuka (2008).
First, we assume collisional coupling time is much smaller than dynamic and ionisation and recombination term, then Eqs. (A.1) to (A.7) become where R_{g} and μ_{β} are the gas constant and mean molecular weight. We can solve these equations analytically. By using this analytic solution, we can analytically integrate the equations over temporal spacing Δt and the physical variables at time t_{0} + Δt are given by
Assuming the partial derivative equations are in the form (A.29)where U,F_{mhd}, and S_{col} are state vector, flux vector of MHD equations, and source vector of collisional coupling term.
By using the operator splitting method, the time integration is performed in the following two steps: In the fourstep RungeKutta time integration, to get updated (t = t_{0} + Δt) state vector U_{1} from state vector U_{0}(t = t_{0}), we use the scheme below as follows: If the collisional temporal scale is much smaller than the MHD temporal scale, the analytic integration of collisional term should interrupt RungeKutta time integration. Therefore, fourstep RungeKutta time integration with analytic integration of collisional term becomes
Figure A.1 shows the result of numerical test of our semianalytic method. This is a 0dimensional simulation to test the thermal coupling of the plasma and neutrals. At the start of the simulation, neutral temperature T_{n} = 2 and plasma temperature T_{p} = 0.5. This difference of temperature between neutrals and plasma decreases with time as a result of the thermal coupling terms in Eqs. (3) and (7). The upper panel of the figure shows the results of explicit methods using a temporal spacing that is comparable to the collisional timescale. The closer the timestep to the collsional timescale, the more the results deviate from the analytic solution. On the other hand, the results of our semianalytic method, shown in the lower panel of the figure, correspond to the analytic solution even if temporal spacing is comparable to the collisional timescale. Therefore, our semianalytic method is useful for reducing the computational time by increasing the temporal spacing Δt in a partially ionised plasma simulation. Our test suggest that time steps of up to 10τ_{c} can be accurately solved with this method. This method also has the great benefit of removing the numerical integration of the collisional terms, which are stiff and as such need a safety factor in the numerical integration of approximately 0.2 to maintain accuracy.
Fig. A.1 Comparison between explicit method (upper panel) and analytic method (bottom panel). Threedotted lines and dashed lines show analytic solution of neutral temperature and plasma temperature. Cross, triangle, and square symbols show numerical solutions with temporal spacing Δt of 0.25, 0.5, and 1.0 of the collisional time τ_{c}. 

Open with DEXTER 
All Figures
Fig. 1 Spatial distribution in the x direction of B_{y}a); gas pressure b); density c); v_{x}d); v_{y}e); and temperature f) for the ideal MHD simulation. There is a fastmode rarefaction wave at approximately x/t = 1 that drives a flow of material towards the shock front at x/t = 0.175. 

Open with DEXTER  
In the text 
Fig. 2 Spatial distribution in the x direction of B_{y}a); gas pressure b); density c); v_{x}d); v_{y}e); and temperature f) for the neutral (red) and ionised (blue) fluids at time t = 1τ. The green line indicates the total (pressure and density) or the difference (x and y velocities and temperature) for the two fluids. The dashed black line shows the reference MHD solution. The ionised plasma has taken on characteristics that are similar to the ideal MHD solution (but with faster wave speeds); the neutrals, however, are undergoing a violent expansion similar to that of a 1D point explosion. 

Open with DEXTER  
In the text 
Fig. 3 Panel a): temporal evolution of the temperature; panel b): magnitude of the heating terms relating to collisional coupling at x = 0 for a timescale that follows the transition from the initiation to the weak coupling regimes (see Sect. 4.2). The normalisation of the power is . Because the temperatures cross at t ~ 3.5τ the thermal damping term become 0 at this point then reverse sign (i.e. the ionised fluid is losing heat to the neutral fluid before this time, but gaining heat from the neutral fluid after this point). 

Open with DEXTER  
In the text 
Fig. 4 Spatial distribution in the x direction of B_{y}a); gas pressure b); density c); v_{x}d); v_{y}e); and temperature f) for the neutral (red) and ionised (blue) fluids at time t = 10τ. The green line indicates the total (pressure and density) or the difference (x and y velocities, and temperature) for the two fluids. The dashed black line indicates the reference MHD solution. 

Open with DEXTER  
In the text 
Fig. 5 Spatial distribution in the x direction of B_{y}a); gas pressure b); density c); v_{x}d); v_{y}e); and temperature f) for the neutral (red) and ionised (blue) fluids at time t = 100τ. The green line indicates the total (pressure and density) or the difference (x and y velocities, and temperature) for the two fluids. The dashed black line indicates the reference MHD solution. 

Open with DEXTER  
In the text 
Fig. 6 Spatial distribution in the x direction of B_{y}a); gas pressure b); density c); v_{x}d); v_{y}e); and temperature f) for the neutral (red) and ionised (blue) fluids at time t = 2000τ. The green line indicates the total (pressure and density) or the difference (x and y velocities, and temperature) for the two fluids. The dashed black line indicates the reference MHD solution. The solution for the partially ionised case has reverted to that of the ideal MHD case outside of the wave fronts. 

Open with DEXTER  
In the text 
Fig. 7 Quasiselfsimilar evolution of the shockrarefaction wave system. To highlight how the system is tending towards selfsimilarity, the xaxis is normalised by the time t of the snapshot of the simulation, therefore shows the velocity at which the waves are propagating in normalised units. As the time is increased, the system tends towards a selfsimilar state. The reference time is taken at t_{0} = 156 000τ and the distributions for t_{0}, t_{0}/ 2, t_{0}/ 4, t_{0}/ 8, and t_{0}/ 16 are plotted. The reference ideal MHD solution is indicated by the black dashed line. 

Open with DEXTER  
In the text 
Fig. 8 Drift velocity through the shock front for the simulation with β = 0.3 and ξ_{i} = 0.01. The xaxis is shifted by x_{s}, where x_{s} is taken as the point in the shock where the drift velocity is at its peak, so that both distributions are aligned. It is clear that the distribution of the drift velocities is the same at both times indicating that the shock front has reached a steady state. The reference time is taken at t_{0} = 120 000τ. 

Open with DEXTER  
In the text 
Fig. 9 Spatial distributions of v_{xn}a) and v_{yn}b) for β = 0.1 and ξ_{i} = 0.1, 0.01, 0.001, 0.0001, and 0.00001 denoted by the green, purple, pink, red, and turquoise lines, respectively. 

Open with DEXTER  
In the text 
Fig. 10 Spatial distributions of v_{xn} − v_{xp}a) and v_{yn} − v_{yp}; b) for ξ_{i} = 0.1, 0.01, 0.001, 0.0001, and 0.00001. 

Open with DEXTER  
In the text 
Fig. 11 Thickness of the shock region as a function of ξ_{i} for β = 0.3 (green), 0.1 (blue), 0.01 (pink), and 0.001 (red). The black dashed line shows a power law as a function of ξ_{i} with exponent − 1.2 

Open with DEXTER  
In the text 
Fig. 12 Spatial distributions of v_{xn}a) and v_{yn}b) for ξ = 0.01 and β = 0.1, 0.01 and 0.001 denoted by the green red and blue lines, respectively. 

Open with DEXTER  
In the text 
Fig. 13 Spatial distribution in the x direction around the shock front and the coupling region of a) B_{y}; b) gas pressure; c) density; d)v_{x}; e) v_{y}; and f) the temperature for the neutral (red) and ionised (blue) fluids at time t = 15 000τ where β = 0.01 and ξ_{i} = 0.01. The green line indicates the total (pressure and density) or the difference (x and y velocities and temperature) for the two fluids. The plasma pressure and density are increased by a factor of 50 to make their distributions clearly visible. The xaxis has been shifted to set the origin at the shock front. 

Open with DEXTER  
In the text 
Fig. 14 Spatial distributions of v_{xn} − v_{xp}a) and v_{yn} − v_{yp}; b) for β = 0.1, 0.01 and 0.001 denoted by the green red and blue lines, respectively. 

Open with DEXTER  
In the text 
Fig. 15 Thickness (in units of V_{A}/ν_{T}) of the shock for the simulations with ξ_{i} = 0.01 as a function of β (left) and , the hydrodynamic Mach number of the neutral fluid in the inflow region (right). The dependence of the thickness of the shock as a function of the plasma β can be explained by the Mach number of the inflow region to which this relates to the compression at the shock front. 

Open with DEXTER  
In the text 
Fig. 16 Spatial distribution of  v_{xn}/Cs_{n}  and  v_{xp}/V_{A}  for a cshock (β = 0.3 and ξ_{i} = 0.01) in panel a), a weak jshock (β = 0.1 and ξ_{i} = 0.01) in panel b), and a jshock (β = 0.01 and ξ_{i} = 0.01) in panel c) around the shock front (calculated in the shock frame). The xaxis is shifted by x_{s}, where x_{s} is taken as the point in the shock where the drift velocity is at its peak, so that all distributions are aligned. 

Open with DEXTER  
In the text 
Fig. 17 Spatial distribution in the x direction around the shock front of a) B_{y}; b) the gas pressure, c) the density; d)v_{x}; e) v_{y}; and f) the temperature for the neutral (red) and ionised (blue) fluids at time t = 15 000τ where β = 0.01 and ξ_{i} = 0.01. The green line indicates the total (pressure and density) or the difference (x and y velocities and temperature) for the two fluids. The plasma pressure and density are increased by a factor of 50 to make their distributions clearly visible. 

Open with DEXTER  
In the text 
Fig. 18 Relation between the total energy that goes to heating the fluids through frictional heating as a package of fluid passes through the shock front for panel a)ξ_{i}; and panel b) β for ξ_{i} = 0.01. 

Open with DEXTER  
In the text 
Fig. A.1 Comparison between explicit method (upper panel) and analytic method (bottom panel). Threedotted lines and dashed lines show analytic solution of neutral temperature and plasma temperature. Cross, triangle, and square symbols show numerical solutions with temporal spacing Δt of 0.25, 0.5, and 1.0 of the collisional time τ_{c}. 

Open with DEXTER  
In the text 