EDP Sciences
Free Access
Issue
A&A
Volume 591, July 2016
Article Number A112
Number of page(s) 16
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201628215
Published online 22 June 2016

© 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 post-reconnection. 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 slow-mode shocks created by magnetic field relaxation. Priest & Forbes (1986) developed a general class of steady-state two-dimensional (2D) reconnection models that included both the diffusion-dominated 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 slow-mode magnetohydrodynamic (MHD) shock known as the switch-off 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 slow-mode 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 slow-mode shocks forming (e.g. Longcope & Bradshaw 2010; Takasao et al. 2015).

Another example that highlights the importance of reconnection-driven, slow-mode 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 slow-mode 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). Fast-mode 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 c-shock; Draine 1980); this allows radiative cooling to more efficiently cool the gas as it shocks, significantly reducing the temperature of the post-shock 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 j-shock). 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 steady-state shock solutions to model molecular line intensities and predict the magnetic field strength of BN-KL. 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 Harris-type current sheet under the influence of partial ionisation, through the one-fluid, ambipolar diffusion approximation, was found to form a power law with the magnetic field scaling as Bx1 / 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 (lnew) that can be estimated as lnew = loldBg/Bex, where Bg and Bex are the original guide and external magnetic fields, respectively (Singh et al. 2015). A similar problem related to the current sheet thickness of the Kippenhahn-Schlüter prominence model under ion-neutral drift was investigated by Hillier et al. (2010), finding that the thickness of the current sheet under ambipolar diffusion would become lnew = loldBhor/Bver, where Bhor and Bver 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 multi-fluid simulations, the coupling and decoupling of ion and neutral fluids in a reconnecting current sheet, finding that the reconnection inflow became decoupled, but the high-velocity 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 Hall-mediated 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 slow-mode reconnection shocks in a partially ionised plasma are yet to be understood. In this paper we look at a simple one-dimensional (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 charge-neutral, 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 Pn = ρnRgTn and Pp = 2ρpRgTp, 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 ion-neutral 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 slow-shocks 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(T0)ρtot = 1 = τ-1 and the bulk Alfvén velocity the normalisation of the magnetic field is used such that . The value of αc(T0) is normalised to one, but αc(Tn,Tp) 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 Lnorm = VA/ν. 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 slow-mode, switch-off 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 one-half 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 = xmax 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 pressure-scale 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 self-similarity, 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 Sedov-Taylor point explosion Sedov 1959).

thumbnail Fig. 1

Spatial distribution in the x direction of Bya); gas pressure b); density c); vxd); vye); and temperature f) for the ideal MHD simulation. There is a fast-mode 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 By magnetic field (the Bx magnetic field is always constant to satisfy ∇·B = 0 and Bx = 0.3 in this case), gas pressure, density, temperature, vy and vx. 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 fast-mode rarefaction wave and a slow-mode shock wave. The fast-mode 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 post-wave values for By, pressure, density, vx and vy.

Post shock there is a high-velocity jet in the y direction; this can be seen as analogous to a reconnection jet. The post-shock velocity becomes vx = 0 in the simulation frame, as does the post-shock y-component of the magnetic field (By). 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 vxU/vxD, 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 post-shock 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 pU in this model for low plasma β.

thumbnail Fig. 2

Spatial distribution in the x direction of Bya); gas pressure b); density c); vxd); vye); 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

thumbnail 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

thumbnail Fig. 4

Spatial distribution in the x direction of Bya); gas pressure b); density c); vxd); vye); 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 fast-mode rarefaction wave and a slow-mode 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 fast-mode 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 slow-mode 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 vD in the jet region.

The large vD in the jet region results in the nonlinear terms of the coupling (those of order ) becoming much larger than those of order vD. Therefore, the coupling of the fluids comes about via the nonlinear coupling terms in the energy equations. The left-hand 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(Tn,Tp)ρnρp(vnvp)2 (solid line) and thermal damping term that drives thermal equilibrium given by 3αc(Tn,Tp)ρnρp(pn/ρnpp/ 2ρp) (dashed line) and the two terms for the work performed on the fluid by the drift velocity given by | αc(Tn,Tp)ρnρpvn·(vnvp) | (blue line) and | αc(Tn,Tp)ρnρpvp·(vnvp) | (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 non-reversible 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, high-temperature region, which, other than the external plasma not being cold, results in a Sedov-Taylor like explosive expansion of the inner high-temperature 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 post-shock 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 pre-shock region and post-shock 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 neutral-shock 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

thumbnail Fig. 5

Spatial distribution in the x direction of Bya); gas pressure b); density c); vxd); vye); 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 c-shock, that the high levels of drift velocity are found (~0.15VA). The two waves that we expect to form in this system, the slow-mode shock and fast-mode rarefaction wave, are visible, but are yet to form as completely distinct elements.

It can be seen that in the region of the c-shock, the neutral fluid is undergoing the velocity transition before the ionised fluid. It is worth noting that this is a clear difference from the c-shocks studied that are relevant to the ISM. This leads to the pre-heating in the shock front relative to the reference MHD solution.

4.4. Quasi-self-similar state

Once the system has sufficiently evolved then it reaches a state that is approaching that of self-similarity. 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 fast-mode rarefaction wave and the slow-mode shock. The pre- and post-shock 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 post-shock 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.

thumbnail Fig. 6

Spatial distribution in the x direction of Bya); gas pressure b); density c); vxd); vye); 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 steady-state 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 steady-state shock. In the following, we provide an analysis of why this quasi-self-similar 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 self-similarity. The x-axis 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 post-shock region evolve in a self-similar 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.

thumbnail Fig. 7

Quasi-self-similar evolution of the shock-rarefaction wave system. To highlight how the system is tending towards self-similarity, the x-axis 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 self-similar state. The reference time is taken at t0 = 156 000τ and the distributions for t0, t0/ 2, t0/ 4, t0/ 8, and t0/ 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 (WR) to be given by the sum of component that comes from the expansion of the system (At) and diffusion of the system (Bt1 / 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 self-similar 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 fast-mode wave speeds (in the simulation frame) upstream and downstream of the wave gives VfU = 1.150VA and VfD = 0.926VA, therefore A = 0.223VA.

thumbnail Fig. 8

Drift velocity through the shock front for the simulation with β = 0.3 and ξi = 0.01. The x-axis is shifted by xs, where xs 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 t0 = 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 t0 and exiting at time t1. The amount of momentum transferred to the neutral fluid by collisions with the ionised fluid (ρnvxn)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 = t1t0) to cross the wave front gets larger, then the average velocity drift that the fluid packet experiences (vxnvxp) 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.

thumbnail Fig. 9

Spatial distributions of vxna) and vynb) 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. Steady-state 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 slow-mode shock and not a fast-mode 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 fast-mode 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 slow-mode shock has a decrease in magnetic pressure that is partly balanced by an increase in gas pressure, it is the neutral fluid that is over-pressured 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 quasi-self-similar 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 two-dimensional 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 tsnap = 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 (Tshock) would roughly scale as Tshock = Vshock/αcρtotξi (e.g. Draine & McKee 1993). It is clear that this is not the case, but interestingly, for the vx velocity the exit point of the shock position is a point of the system that evolves self-similarly (see Fig. 7); therefore, the time scaling used for the snapshot multiplied by ξi relates to the self-similar 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 vy 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 y-direction 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.

thumbnail Fig. 10

Spatial distributions of vxnvxpa) and vynvyp; 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

thumbnail 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 c-shocks to that of producing j-shocks 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 Rankine-Hugoniot relation. For an ideal gas with adiabatic index γ = 5 / 3, the ratio of the upstream and downstream velocities in the shock frame is uU/uD → 4 as the shock gets progressively stronger. For the case of β = 0.01, we have a shock velocity of ~0.15VA and an inflow velocity in the simulation frame of ~− 0.2VA, which would give an upper estimate for the post-shock velocity in the simulation frame of ~0.06VA, i.e. the velocity undergoes a positive overshoot above the vx = 0 level that is expected for the final state after the fluids have recoupled. The actual post-shock velocities shown in Fig. 12 are ~0.02VA 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 vx = 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.

thumbnail Fig. 12

Spatial distributions of vxna) and vynb) 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 vx, P, and ρ for both the neutral and ionised fluid. This is a distinct difference from the j-shock 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 j-shock 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 post-shock value as the neutral fluid shocks, but because of the diverging neutral velocity field the density decays to the post-shock value.

thumbnail Fig. 13

Spatial distribution in the x direction around the shock front and the coupling region of a) By; b) gas pressure; c) density; d)vx; e) vy; 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 x-axis 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 c-shock to a j-shock, the structure of the drift velocity is characterised by the sharp jump of the shock and a monotonic drop towards the 0.

thumbnail Fig. 14

Spatial distributions of vxnvxpa) and vynvyp; 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.

thumbnail Fig. 15

Thickness (in units of VA/ν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

thumbnail Fig. 16

Spatial distribution of | vxn/Csn | and | vxp/VA | for a c-shock (β = 0.3 and ξi = 0.01) in panel a), a weak j-shock (β = 0.1 and ξi = 0.01) in panel b), and a j-shock (β = 0.01 and ξi = 0.01) in panel c) around the shock front (calculated in the shock frame). The x-axis is shifted by xs, where xs 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 vxn) of the neutral fluid and the Alfvénic Mach number in the shock frame of the ionised fluid (calculated using vxp) for β = 0.3, 0.1 and 0.01. This choice of plasma β gives a subsonic c-shock, a transonic weak j-shock, and a strong j-shock. The x-axis is shifted by xs, where xs 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 j-shock to occur there has to be a supersonic inflow velocity, otherwise a c-shock 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 c-shock. The Mach number drops to less than 0.3 for all cases shown at post-shock . At all times the ionised fluid is sub-Alfvénic.

5.2.3. Fine structure of the j-shock

In the shock frame, the ionised fluid is travelling towards the shock front, in this case at a velocity of ~0.4VA, 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 = Vion/α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 x-axis 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.

thumbnail Fig. 17

Spatial distribution in the x direction around the shock front of a) By; b) the gas pressure, c) the density; d)vx; e) vy; 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

thumbnail 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 fast-mode rarefaction wave and a slow-mode 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 quasi-self-similar state. In the quasi-self-similar stage, the shock forms a steady state but the rarefaction wave undergoes a diffusive evolution as it tends towards its final state. A high-speed (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 c-shock solution for low Mach number inflow and a j-shock 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 c-shock and sharp transitions in both fluids for the j-shock. Though in the model used the inflow is driven by a fast-mode rarefaction wave, which, as it is a compressional wave, gives a smaller inflow velocity when it can expand in three-dimensions (3D) than in 1D, so the dependence of the transition between the c- and j-shocks 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 fast-mode MHD shocks, where both the magnetic field strength and gas pressure increase downstream of the shock front, but we have been investigating a slow-mode 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 j-shock solutions for both studies is one area of similarity, though the conditions for the appearance of the j-shock solution differ between the two studies. For this study, j-shocks 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 shock-driven 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 slow-mode 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 ion-neutral 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 vyn 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 E-infrastructure 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 E-Infrastructure.

References

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 ion-electron plasma fluid are written as Equation (A.7) contains the Ambipolar diffusion term (ηAB2J). 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 fourth-order four-Runge-Kutta 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 Rg 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 t0 + Δt are given by

Assuming the partial derivative equations are in the form (A.29)where U,Fmhd, and Scol 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 four-step Runge-Kutta time integration, to get updated (t = t0 + Δt) state vector U1 from state vector U0(t = t0), 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 Runge-Kutta time integration. Therefore, four-step Runge-Kutta time integration with analytic integration of collisional term becomes

Figure A.1 shows the result of numerical test of our semi-analytic method. This is a 0-dimensional simulation to test the thermal coupling of the plasma and neutrals. At the start of the simulation, neutral temperature Tn = 2 and plasma temperature Tp = 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 semi-analytic 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 semi-analytic 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.

thumbnail Fig. A.1

Comparison between explicit method (upper panel) and analytic method (bottom panel). Three-dotted 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

thumbnail Fig. 1

Spatial distribution in the x direction of Bya); gas pressure b); density c); vxd); vye); and temperature f) for the ideal MHD simulation. There is a fast-mode 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
thumbnail Fig. 2

Spatial distribution in the x direction of Bya); gas pressure b); density c); vxd); vye); 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
thumbnail 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
thumbnail Fig. 4

Spatial distribution in the x direction of Bya); gas pressure b); density c); vxd); vye); 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
thumbnail Fig. 5

Spatial distribution in the x direction of Bya); gas pressure b); density c); vxd); vye); 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
thumbnail Fig. 6

Spatial distribution in the x direction of Bya); gas pressure b); density c); vxd); vye); 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
thumbnail Fig. 7

Quasi-self-similar evolution of the shock-rarefaction wave system. To highlight how the system is tending towards self-similarity, the x-axis 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 self-similar state. The reference time is taken at t0 = 156 000τ and the distributions for t0, t0/ 2, t0/ 4, t0/ 8, and t0/ 16 are plotted. The reference ideal MHD solution is indicated by the black dashed line.

Open with DEXTER
In the text
thumbnail Fig. 8

Drift velocity through the shock front for the simulation with β = 0.3 and ξi = 0.01. The x-axis is shifted by xs, where xs 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 t0 = 120 000τ.

Open with DEXTER
In the text
thumbnail Fig. 9

Spatial distributions of vxna) and vynb) 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
thumbnail Fig. 10

Spatial distributions of vxnvxpa) and vynvyp; b) for ξi = 0.1, 0.01, 0.001, 0.0001, and 0.00001.

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 12

Spatial distributions of vxna) and vynb) 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
thumbnail Fig. 13

Spatial distribution in the x direction around the shock front and the coupling region of a) By; b) gas pressure; c) density; d)vx; e) vy; 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 x-axis has been shifted to set the origin at the shock front.

Open with DEXTER
In the text
thumbnail Fig. 14

Spatial distributions of vxnvxpa) and vynvyp; b) for β = 0.1, 0.01 and 0.001 denoted by the green red and blue lines, respectively.

Open with DEXTER
In the text
thumbnail Fig. 15

Thickness (in units of VA/ν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
thumbnail Fig. 16

Spatial distribution of | vxn/Csn | and | vxp/VA | for a c-shock (β = 0.3 and ξi = 0.01) in panel a), a weak j-shock (β = 0.1 and ξi = 0.01) in panel b), and a j-shock (β = 0.01 and ξi = 0.01) in panel c) around the shock front (calculated in the shock frame). The x-axis is shifted by xs, where xs 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
thumbnail Fig. 17

Spatial distribution in the x direction around the shock front of a) By; b) the gas pressure, c) the density; d)vx; e) vy; 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
thumbnail 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
thumbnail Fig. A.1

Comparison between explicit method (upper panel) and analytic method (bottom panel). Three-dotted 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.