High order accurate shock capturing schemes for two-component Richtmyer-Meshkov instabilities in compressible magnetohydrodynamics

We design a conservative and entropy satisfying numerical scheme to perform numerical simulations of two component Richtmyer-Meshkov (RM) instabilities in compressible magnetohydrodynamics (MHD). We first formulate a conservative model of a two-component compressible MHD fluid ruled under two ideal gases with different adiabatic exponents. The formulation includes a level set function that allows to evolve the two components of the plasma in a conservative and consistent way. We present a set of examples including two-component Riemann problems and high Mach shock wave interactions with entropy contact waves that validate the high order accurate numerical scheme. We observe that turbulent regimes are completely developed in different examples where shocks, contacts and rarefactions waves propagate with correct speed.


I. Introduction
Richtmyer-Meshkov (RM) instabilities 9,10 arise when a shock wave encounters a fluid interface.Small perturbations at the interface grow into nonlinear structures in the form of bubbles and spikes.The impulsive acceleration generated by the shock wave induces the perturbation to increase linearly in time.Interesting scenarios where the physics of these complex wave phenomena appear include astrophysical flows where the interaction of a bow shock wave with an interstellar medium ruled by an ideal gas under different adiabatic exponents occurs.In particular Wu and Roberts 16 considered the RM instability in MHD and discussed its possible role in the dynamics of the magnetosphere.They modeled the magnetopause as a tangential discontinuity, i.e., a contact discontinuity such that the normal component of the magnetic field to the interface vanishes.They considered the single mode RM instability based on an ideal gas with adiabatic exponent γ = 5/3.
Numerical simulations of emerging instabilities in the presence of shock wave interactions with unstable contact discontinuities in plasmas represent a computational challenge.In this research work we formulate a conservative model of a two-component compressible MHD fluid ruled under two ideal gases with different adiabatic exponents.The formulation includes a level set function that allows to evolve the two components of the plasma in a conservative way in a similar fashion as proposed by Marquina and Mulet 8 for compressible hydrodynamic flows.In order to explore the complex dynamics of the two-component MHD model through numerical simulations we design a conservative and entropy satisfying numerical scheme to approximate the solution of the proposed two component MHD model and perform numerical simulations of high Mach shock wave interactions in one and two dimensions.
This paper is organized as follows.In section II we propose a two-component compressible MHD model for the dynamics of the mixture of two ideal gases.In section III we present a spectral decomposition of the two-component MHD model for the numerical design.Section IV is devoted to the design of a high order shock capturing characteristic-based numerical scheme for the approximation of the two-component MHD model proposed.In section V we present a set of numerical experiments including two-component Riemann problems and high Mach shock wave interactions in one and two dimensions.We draw our conclusions in section VI.

II. A two-component compressible MHD model
In this section we propose a model for the dynamics of the mixture of two ideal gases in compressible MHD.Let a level curve φ represent the mass fraction of the first component and 1 − φ the mass fraction of the second one.We assume that both components are in thermal equilibrium and are calorically perfect gases with specific heats at constant volume Cv 1 and Cv 2 , specific heats at constant pressure Cp 1 and Cp 2 and ratios of specific heats γ 1 , γ 2 .The ratio of specific heats of the mixture of gases can be expressed as Let us consider a two-dimensional hyperbolic system of conservation laws of the form ) are the fluxes.The two-dimensional ideal MHD equations is a system of the form (2). We model the dynamics of the mixture of the compressible MHD equations adding an equation to the MHD system formed by the equations representing conservation of mass, momentum, magnetic field and energy.The new equation expresses mass conservation of the first component which coupled to the conservation of the total mass implies conservation of the second component.The vector of conserved variables becomes u = (ρ, ρv, B, E, ρφ) T being v = (u, v, w) the velocity field, B = (B x , B y , B z ) the magnetic field and E = 1 2 ρq 2 + 1 2 B 2 + P γ(φ)−1 the total energy where q 2 = u 2 + v 2 + w 2 and B 2 = B 2 x + B 2 y + B 2 z .The hydrodynamic pressure P is defined through the ideal gas equation of state (EOS) in terms of the mass fraction φ as P = (γ(φ) − 1)ρǫ where ǫ is the specific internal energy.The total pressure is represented as P * = P + 1 2 B 2 .The fluxes for the two-component model in each direction are then defined as In addition to this system of equations the magnetic field satisfies the divergence-free constraint In order to explore the complex dynamics of the presented two-component MHD model we propose a characteristic based numerical scheme that is designed taking into account full information of the nonlinearity of the wave structure of the system through the spectral decomposition of the Jacobians of the fluxes.In the next section we propose a complete system of eigenvectors and the corresponding eigenvalues of the Jacobians for the two-component MHD fluxes in terms of the thermodynamic magnitudes and the mass fraction.

III. A spectral decomposition of the two-component MHD fluxes
A system of equations is hyperbolic if the Jacobians of the fluxes are diagonalizable matrices with real eigenvalues and a complete set of eigenvectors in each neighborhood of the solution (Lax 7 ).The diagonalization of the Jacobians decouples the original hyperbolic system in m scalar conservation laws defining the so-called characteristic fields and their corresponding characteristic fluxes, (Lax 6 ).
The eigenvalues of the Jacobian represent the characteristic speeds of the characteristic fluxes analogous to the first derivative of the flux of a scalar conservation law.
Let us consider a hyperbolic system of type (2).Let λ f 1 (u), • • • , λ f m (u) be the eigenvalues of the Jacobian f ′ (u) counting each one as many times as its multiplicity and u)} a complete system of right and left eigenvectors respectively diagonalizing f ′ (u) such that Following this formalism let us study the spectral structure of the two-component MHD equations.There are nine characteristic fields which characteristic velocities are: z and the acoustic sound speed a = γ(φ)P ρ .In the calculation of the eigenvectors of the Jacobians we follow a similar procedure carried out by Serna 13 where the eigenvectors are obtained as smooth functions of the conserved variables.We show a complete system of eigenvectors of the Jacobian of f (x direction).The corresponding decomposition for the Jacobian of g (y direction) is similar.
We define sgn(t) = 1 for t ≥ 0 and sgn(t) = −1 otherwise and set β y and β z values from the expressions otherwise The eigenvectors associated to λ 2 , λ 4 , λ 5 , λ 7 and λ 9 , are The eigenvectors associated to λ 1 , λ 3 , λ 6 and λ 8 , can be expressed in an unified way for k = 1, 3, 6, 8 as where c, c, α and ᾱ are defined such that the eigenvectors are continuous functions with respect to the conserved variables as x α f ; otherwise α f and α s are defined from the following expressions,

; otherwise
The system of eigenvectors for the Jacobian of g can be obtained by interchanging B x by B y , u by v, the second component by the third and the fifth by the sixth.
Next we design a numerical scheme based on the local characteristic decomposition described above.

IV. High order characteristic-based numerical scheme
In this section we extend the general purpose shock capturing characteristic-based entropy-fix upwind scheme proposed for single mode compressible MHD equations in 13 to the two-component MHD conservative formulation introduced in section II.The numerical scheme in 13 handles an entropy correction through a prescription of a local viscosity ensuring convergence to the entropy solution of the system of equations.We use the proposed spectral decomposition of the Jacobians of the fluxes for the two-component case and design a first order shock capturing numerical scheme following the interface splitting strategy used in the Marquina flux formula. 4 numerical scheme in conservation form for a system of conservation laws in two dimensions can be written as where u n jk ≈ u(x j , y k , t n ) is a numerical approximation of the solution in the computational cell (x j , y k ) where x j = j∆x, y k = k∆y and t n = n∆t where ∆x, ∆y and ∆t are the spatial and time step sizes respectively and f and g represent the numerical fluxes in each direction consistent with the fluxes of the system.We perform the local characteristic decomposition of the system decoupling the equations in linearly independent characteristic fields.
Numerical schemes for hyperbolic conservation laws based on characteristic field decomposition are required to satisfy Rankine-Hugoniot relations to approximate the numerical fluxes of the scheme 6,7 .Numerical schemes using one linearization at interfaces need of explicit formulas to compute an average state at interfaces. 11There exist exact formulas for ideal MHD system only for the case where the adiabatic exponent γ is equal to 2 (as used in the paper by Brio and Wu 2 ).For other values of γ there are not such formulas to satisfy Rankine-Hugoniot relations.Most numerical schemes use arithmetic averages instead to define what can be non-physical intermediate states. 2,5 ince we are focusing on the approximation of interfaces between fluids with different adiabatic exponent, we will use Marquina's flux splitting that allows to satisfy approximately Rankine-Hugoniot relations at interfaces.We approximate the numerical fluxes at cell interfaces by means of two linearizations at each side of the interface following the interface splitting strategy used in the Marquina flux formula. 4The strategy avoids averaging at interfaces and therefore thermodynamic inconsistencies at intermediate states.It has been proved to be very convenient to avoid numerical instabilities in the simulation of complex fluids as hydrodynamic and relativistic flows. 3,4,8,12,14 To mpute f in terms of two linearizations at each interface we use the first order flux splitting formula where ψ p + and ψ p − represent the lateral numerical characteristic fluxes.These are computed at the interface following the entropy-fix upwind procedure proposed in 13 from the local characteristic fluxes and variables The calculation of g at interface (x j , y k+ 1 2 ) is similar.In order to achieve high order accuracy in space we perform an extension of the basic first order numerical scheme applying a third order accurate reconstruction function on local characteristic fluxes and variables following the so-called "flux formulation".This is performed evaluating at the interface the reconstruction function that is determined via primitive function and satisfies the conservation property as proposed by Shu and Osher 15 .We use the third order accurate piecewise hyperbolic reconstruction procedure presented in 12 and used in 13 for single mode MHD.We apply a high order total variation diminishing (TVD) Runge-Kutta time stepping procedure for the integration in time 15 .
To satisfy the divergence-free constraint we apply the correction on the magnetic field at the end of every time step following the prescription proposed by Brackbill et al. 1 Next we present preliminary numerical results validating the above shock-capturing scheme for the approximation of the solution of the proposed two-component compressible MHD.

V. Numerical results
We present a set of numerical experiments for the system of two-component MHD equations in one and two spatial dimensions using the third order accurate numerical scheme proposed in the previous section.
The proposed numerical method is stable under a CFL restriction determined by for one-dimensional experiments and ∆t = C max(|u|+c x f ) ∆x + max(|v|+c y f ) ∆y (9)   for two-dimensional ones where C is 0.8 in our calculations and the maximums are taken over all computational cells.
We first test our scheme with the standard Brio-Wu shock tube problem with two plasmas.In second place we propose an example simulating the interaction of a shock wave with a contact discontinuity between plasmas that obey ideal gas laws with different adiabatic exponents.In the last case we simulate a Richtmyer-Meshkov instability in one and two dimensions.
For each one-dimensional problem we display the profiles of the density, mass fraction, x-component of the velocity, total pressure, y-component of the magnetic field and the fast acoustic impedance.
We observe the latter to better understand the wave propagation of the specific scenarios under study.The acoustic impedance is considered a measure of the resistance to the propagation of hydrodynamic shock waves in a specific medium and in hydrodynamics is defined as Z = ρc where c represents the sound speed.
We use the fast acoustic impedance in compressible plasmas as a measure of this resistance defined as Z = ρc f .The lower the acoustic impedance the faster shock waves propagate through the region of interest.In the two component compressible plasma model under study ruled by ideal gases with different adiabatic constants the fast acoustic impedance of a region allows to predict how fast travels a shock wave when crossing an interface from one gas to the other.

A. Two-component Brio-Wu shock tube Riemann problem
We approximate the solution of the one dimensional Brio-Wu 2 shock tube Riemann problem for two component MHD.This example was proposed for single mode ideal MHD gases to show the formation of a compound wave.The two initial constant states, u L and u R are (ρ, u, v, w, B x , B y , B z , P, φ) = (1, 0, 0, 0, 0.75, 1, 0, 1, 0); x ≤ 0 (0.125, 0, 0, 0, 0.75, −1, 0, 0.1, 1); x > 0 The adiabatic exponents are γ L = 5/3 and γ R = 2 for the left and right states respectively.We solve the one-dimensional MHD system for x ∈ [−1, 1] with N = 800 equally spaced grid points.We evolve the numerical scheme using a CFL number C = 0.5.
In Figure 1 we display the solution at time t = 0.2.We observe that the numerical scheme behaves robust and accurate showing similar dynamics than the ones for the single mode ideal MHD case. 2 We obtain two fast rarefaction waves and a slow compound wave interacting with the tail of the left rarefaction.

B. Two-component Mach 2 shock-contact problem
We propose a one-dimensional shock contact problem that mimic the interaction of a shock wave with the magnetopause represented by an interface that separates two plasmas modeled by the MHD equations ruled by two ideal gas laws with adiabatic exponents γ = 5/3 and γ = 2 respectively.Across the interface the total pressure is constant.We consider the domain interval [−5, 5] and locate the interface at x = 0.At time t = 0 a shock wave is initiated at x = −2.5 with a fast shock Mach number M f = 2 that travels through the γ = 5/3 plasma to the right.Table 1 displays initial data with nonzero quantities.Variables v = w = 0 and B x = B z = 0.
These initial data are obtained from the Rankine-Hugoniot equations derived from the MHD equations presented by Zhuang et al. 18 using as starting values the ones by Wu 17 for γ = 5/3.
Figure 2 displays the solution at time t = 1.1543 computed with 800 grid points and CFL coefficient 0.05.The evolution shows the shock wave located initially at x = −2.5 traveling with constant speed equal to 2.62 towards the interface located at x = 0.This is hit at time 0.9543.From this time the shock wave generates two waves: one transmitted shock wave that travels through the γ = 2 plasma and a reflected rarefaction wave traveling to the left through the γ = 5/3 plasma.Results in figure 2 are the solution at time 0.2 after the shock wave collides the interface.We observe that the fast acoustic impedance in the γ = 2 region (the right hand side of the interface) is much lower than the one in the γ = 5/3 region.This implies that the incident shock wave that travels towards the interface with speed 2.62 generates a transmitted shock wave (after the interaction) that travels at speed 8.37 through the γ = 2 region.
The Mach 2 fast shock wave initially located at x=0.3 hits the bubble at time t = 0.01908.The shock wave interaction with the interface generates a reflected rarefaction wave and a transmitted shock wave traveling through the bubble with speed much faster than the original one.At time t = 0.05492 the transmitted shock wave hits the right hand side interface of the bubble.A reflected rarefaction wave is generated and a new transmitted shock wave travels towards the right outside of the bubble.At this point the wave dynamics inside the bubble becomes complex involving several shock waves and rarefactions interacting.The regime becomes unstable.Figure 4 shows the fully generated complex wave dynamics for the one dimensional case.
In figure 4 we display the results at time t = 0.16 for the density and the absolute value of the magnetic field for the two-dimensional bubble interaction where the complex wave dynamics generated inside and around the bubble is observed.

VI. Conclusion
We have presented a two-component compressible MHD model consisting of including a level set function in the formulation that allows to evolve two components of the plasma in a conservative way.We design a conservative and entropy satisfying characteristic-based numerical scheme to approximate the solution of the proposed two component MHD model and perform numerical simulations.To explore the complex dynamics of the plasma system we introduce a two-component Mach 2 shock-contact problem mimicking the interaction of a shock wave with the magnetopause represented by an interface that separates two plasmas and simulate a Richtmyer-Meshkov instability in one and two dimensions.
where the Alfven velocity c a = |b x | and the fast and slow velocities are given by

Figure 1 .
Figure 1.Two-component Brio-Wu shock tube Riemann problem at time t = 0.2

Table 1 .
Two-component Mach 2 shock-contact problem initial data