Transverse bifurcation of viscous slow MHD shocks

TL;DR

This paper combines analytical and numerical Evans function methods to analyze the stability and transverse bifurcation of multi-D viscous and inviscid slow MHD shocks, revealing a preference for nonplanar traveling wave solutions.

Contribution

It presents the first multi-D numerical Evans function study for viscous MHD shocks and demonstrates the agreement between viscous and inviscid stability transitions.

Findings

Nonplanar traveling waves are typical in slow Lax MHD shocks.

Viscous and inviscid stability transitions align.

First multi-D numerical Evans function analysis for viscous MHD.

Abstract

We study by a combination of analytical and numerical Evans function techniques multi-D viscous and inviscid stability and associated transverse bifurcation of planar slow Lax MHD shocks in a channel with periodic boundary conditions. Notably, this includes the first multi-D numerical Evans function study for viscous MHD. Our results suggest that, rather than a planar shock, a nonplanar traveling wave with the same normal velocity is the typical mode of propagation in the slow Lax mode. Moreover, viscous and inviscid stability transitions appear to agree, answering (for this particular model and setting) an open question of Zumbrun and Serre.

Figures

Tables

Table 1. Table 1 . Shock types in 2-D MHD.

Parameter range	Num of positive e-vals (Left,Right)	Shock type
$h_1>H^{\ast }:=u_1^{-}\sqrt{{\rho }^{-}}$	(4,3)	Lax 2-shock (Slow)
	(5,4)	Lax 1-shock (Fast)
	(5,3)	Doubly overcompressive

Table 2. Table 2. Comparison of exact coefficient with numerical coefficient. Parameters u 1 − = 1 superscript subscript 𝑢 1 1 u_{1}^{-}=1 , ρ − = 1 superscript 𝜌 1 \rho^{-}=1 are fixed; recall that u 1 + = u 1 − / R superscript subscript 𝑢 1 superscript subscript 𝑢 1 𝑅 u_{1}^{+}=u_{1}^{-}/R . For each value of γ 𝛾 \gamma and R 𝑅 R we record the exact answer, the numerical estimate (est), and the relative error (err) between the two. Numerical estimates of the coefficient were determined by computing the roots λ ​ ( h 1 ) 𝜆 subscript ℎ 1 \lambda(h_{1}) of the Lopatinsky determinant for several values of h 1 ⊂ [ 2 , 16 ] subscript ℎ 1 2 16 h_{1}\subset[2,16] and then using curve fitting of log ⁡ ( λ ) 𝜆 \log(\lambda) and log ⁡ ( 1 / h 1 ) 1 subscript ℎ 1 \log(1/h_{1}) . The maximum relative error is 4.70e-1 and the average relative error is 4.83e-2.

$\gamma$ $\backslash R$	1.5	2	2.5	3	3.5	4	4.5	5	5.5	6
$\gamma =5/4$, exact	7.40e-2	1.01e-1	1.13e-1	1.18e-1	1.20e-1	1.20e-1	1.19e-1	1.18e-1	1.16e-1	14e-1
$\gamma =5/4$, est	6.65e-2	9.76e-2	1.10e-1	1.18e-1	1.16e-1	1.16e-1	1.16e-1	1.18e-1	1.13e-1	1.10e-1
$\gamma =5/4$, err	1.01e-1	3.64e-2	2.98e-2	2.54e-4	2.64e-2	2.72e-2	2.16e-2	2.26e-3	2.29e-2	3.52e-2
$\gamma =7/5$, exact	7.33e-2	0.100	1.11e-1	1.15e-1	1.17e-1	1.17e-1	1.16e-1	1.14e-1	1.13e-1	1.11e-1
$\gamma =7/5$, est	6.65e-2	9.97e-2	1.10e-1	1.15e-1	1.16e-1	1.16e-1	1.14e-1	1.15e-1	1.10e-1	1.10e-1
$\gamma =7/5$, err	9.44e-2	2.32e-3	8.06e-3	6.05e-3	4.70e-1	3.36e-3	1.51e-2	4.11e-3	2.02e-2	6.98e-3
$\gamma =5/3$ exact	7.24e-2	9.77e-2	1.08e-1	1.11	1.12e-1	1.12e-1	1.11e-1	1.09e-1	1.07e-1	1.05e-1
$\gamma =5/3$ est	6.92e-2	9.84e-2	1.09e-1	1.12	1.16e-1	1.11e-1	1.12e-1	112e-1	1.08e-1	1.05e-1
$\gamma =5/3$ err	4.37e-2	7.23e-3	7.70e-3	1.73e-3	3.43e-2	3.02e-3	1.01e-2	2.54e-2	4.65e-3	1.88e-3
$\gamma =3$ exact	6.78e-2	8.77e-2	9.40e-2	9.53e-2	9.47e-2	9.32e-2	9.13e-2	8.93e-2	8.73e-2	8.54e-2
$\gamma =3$ est	7.11e-2	9.41e-2	9.98e-2	1.02e-1	1.02e-1	1.05e-1	1.01e-1	9.60e-2	9.84e-2	9.29e-2
$\gamma$ = 3 err	4.89e-2	7.30e-2	6.17e-2	7.26e-2	7.98e-2	1.23e-1	1.10e-1	7.50e-2	1.26e-1	8.86e-2

Table 3. Table 3 . Table providing computational details of the viscous stability study. The fifth through ninth columns respectively show the winding number of the computation, the location of the root (if applicable) computed with absolute tolerance of 5 × 10 − 7 5 superscript 10 7 5\times 10^{-7} , the outer radius of the contour on which the Evans function was computed, the number of points on the contour, and the time in seconds the computation took to run. The Evans function was computed on a semi-annulus with inner radius 10 − 5 superscript 10 5 10^{-5} and outer radius as stated. A * indicates that the radius was taken large enough that curve fitting the Evans function with its asymptotic behavior yields a relative error no greater than 0.2.

$\gamma$	$u_{1+}$	$h_1$	$\xi$	WND	Root	Radius	Num Pnts	Run time
5/3	0.0001	1.1	0.001	0	NA	128	539	197
5/3	0.0001	1.1	0.1	1	1.2280e-04	128	491	162
5/3	0.0001	1.5	0.05	1	7.0364e-05	128	517	171
5/3	0.0001	8	0.2	0	NA	128	543	381
7/5	0.0001	1.5	0.005	0	NA	128	577	182
7/5	0.0001	4	0.05	1	1.1829e-05	128*	599	282
7/5	0.0001	16	0.1	0	NA	128	527	747
5/3	0.01	1.1	0.001	1	2.8902e-05	128*	341	161
5/3	0.01	1.1	0.1	1	0.0024	128	297	128
5/3	0.01	1.5	0.05	1	7.5206e-04	64*	281	128
5/3	0.01	8	0.2	0	NA	128	361	295
7/5	0.01	1.1	0.8	1	0.0050	128*	287	139
7/5	0.01	2	1.6	0	NA	64*	237	132
7/5	0.01	8	1.6	0	NA	8*	133	163
5/3	0.2	1.1	0.005	1	4.7035e-04	16*	137	93.2
5/3	0.2	1.1	0.2	1	0.0137	4*	141	84.2
5/3	0.2	1.5	0.1	1	0.0044	32*	159	112
5/3	0.2	8	0.8	0	NA	8*	171	203
7/5	0.2	1.5	0.05	1	0.0024	32*	165	113
7/5	0.2	4	0.1	1	3.7694e-04	4*	221	138
7/5	0.2	16	0.1	0	NA	8*	251	521

Equations

$${\rho }_t+\mathrm{div}(\rho u)=$$

$$(\rho u{)}_t+\mathrm{div}(\rho u\otimes u-h\otimes h)+\nabla q=$$

$$h_t-\nabla \times (u\times h)=$$

$$\mathrm{div}(h)\Big{|}_{t=0}=0,$$

$$H_{\ast }=u_1^{+}\sqrt{{\rho }^{+}},H^{\ast }=u_1^{-}\sqrt{{\rho }^{-}}$$

$${\rho }_t+\mathrm{div}(\rho u)$$

$$(\rho u{)}_t+\mathrm{div}(\rho u\otimes u-h\otimes h)+\nabla q$$

$$h_t-\nabla \times (u\times h)+\beta \mathrm{div}(h)e_1$$

$$\nabla \times (h\times u)$$

$${\rho }_t+\mathrm{div}(\rho u)=0,(\rho u{)}_t+\mathrm{div}(\rho u\otimes u)+h\times (\nabla \times h)+\nabla p=0,h_t-\nabla \times (u\times h)=0.$$

$$f_0(\mathcal{W}{)}_t+\sum _{j=1}^d{f}_j(\mathcal{W}{)}_{x_j}=0,\mathcal{W}\in {\mathbb{R}}^n,$$

$$f_0$$

$$f_1$$

$$f_2$$

$$A_0{\mathcal{W}}_t+\sum _{j=1}^d{A}_j{\mathcal{W}}_{x_j}=0,A_j:=D{f}_j(\overline{\mathcal{W}}).$$

$${\mathcal{V}}_t+\mathcal{F}(\mathcal{V})=0,$$

$$\mathcal{F}(\mathcal{V}):=\sum _{j=1}^d{\mathcal{F}}_j(\mathcal{V}{)}_{x_j}{\mathcal{F}}_j(\mathcal{V}):=f_j(\mathcal{V}(\mathcal{W}))=f_j(f_0^{-1}(\mathcal{W})),$$

$${\mathcal{V}}_t+\sum _{j=1}^d{\mathcal{A}}_j{\mathcal{V}}_{x_j}=0,{\mathcal{A}}_j:=D{\mathcal{F}}_j(\overline{\mathcal{V}})=A_j{A}_0^{-1}(\overline{\mathcal{W}})).$$

$$\Gamma \mathcal{V}:=\sum _{j=1}^d{\Gamma }_j{\partial }_{x_j}\mathcal{V}\equiv 0,{\Gamma }_j\in {\mathbb{R}}^{m\times n}\equiv \text{constant},$$

$$\Gamma \mathcal{F}(\mathcal{V})=-\mathcal{M}\Gamma \mathcal{V},$$

$$\mathcal{M}:=\sum _{j=1}^d{\mathcal{M}}_j{\partial }_{x_j},{\mathcal{M}}_j\in {\mathbb{R}}^{m\times m}\equiv \text{constant}.$$

$$\sigma \left(\sum _{j=1}^d{k}_j{\mathcal{M}}_j\right)real\; and\; semisimple,\forall k\in {\mathbb{R}}^d.$$

$$\Gamma \mathcal{L}=-\mathcal{M}\Gamma ,$$

$$\mathcal{A}(k)=\sum _{j=1}^d{k}_j{\mathcal{A}}_j^{\pm },\Gamma (k)=\sum _{j=1}^d{k}_j{\Gamma }_j,\mathcal{M}(k)=\sum _{j=1}^d{k}_j{\mathcal{M}}_j,$$

$$\Gamma (k)A(k)=\mathcal{M}(k)\Gamma (k),$$

$$\sigma (A(k))=\sigma (A(k)|\ker \Gamma (k))\cup \sigma (\mathcal{M}(k){|}_{\mathrm{Range}\Gamma (k)}).$$

$$\mathcal{A}(k)=(k_1{A}_1+k_2{A}_2)A_0^{-1}=\left(\begin{array}{ccccc}0& k_1& k_2& 0& 0\\ p_{\rho }{k}_1& 0& 0& -{\bar{h}}_1{k}_1-{\bar{h}}_2{k}_2& {\bar{h}}_2{k}_1-{\bar{h}}_1{k}_2\\ p_{\rho }{k}_2& 0& 0& {\bar{h}}_2{k}_2-{\bar{h}}_2{k}_1& -{\bar{h}}_1{k}_1-{\bar{h}}_2{k}_2\\ 0& -\frac{{\bar{h}}_2}{\bar{\rho }}{k}_2& \frac{{\bar{h}}_1}{\bar{\rho }}{k}_2& \beta k_1& \beta k_2\\ 0& \frac{{\bar{h}}_2}{\bar{\rho }}{k}_1& -\frac{{\bar{h}}_1}{\bar{\rho }}{k}_1& 0& 0\end{array}\right),$$

$$\mathcal{A}(k)=\left(\begin{array}{ccccc}0& 0& 1& 0& 0\\ 0& 0& 0& 0& -{\bar{h}}_1\\ p_{\rho }& 0& 0& 0& 0\\ 0& 0& \frac{{\bar{h}}_1}{\bar{\rho }}& 0& \beta \\ 0& 0& 0& 0& 0\end{array}\right),\mathcal{A}(k){|}_{\ker \Gamma (k)}=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& -{\bar{h}}_1\\ p_{\rho }& 0& 0& {\bar{h}}_1\\ 0& 0& \frac{{\bar{h}}_1}{\bar{\rho }}& 0\end{array}\right).$$

$${\mathcal{F}}_j(\mathcal{V})={\mathcal{N}}_j\mathcal{V}+{\mathcal{Q}}_j(\mathcal{V}),$$

$${\Gamma }_j{\mathcal{N}}_j$$

$${\Gamma }_j{\mathcal{N}}_k+{\Gamma }_k{\mathcal{N}}_j$$

$$\sum _{j,k=1}^d⟨w_j,d^2({\Gamma }_j{\mathcal{Q}}_k+{\Gamma }_k{\mathcal{Q}}_j)w_k⟩=0$$

$$\sum _{j,k=1}^d⟨{\partial }_{x_k}{\partial }_{x_j}\varphi ,{\Gamma }_j{\mathcal{F}}_k(\mathcal{V})⟩=\sum _{j,k=1}^d⟨{\partial }_{x_k}{\partial }_{x_j}\varphi ,{\mathcal{M}}_j{\Gamma }_k\mathcal{V}⟩;a.e.$$

$$\mathrm{div}(h{)}_t+\mathrm{div}(u\mathrm{div}(h))=0,$$

Full text

Research of B.B. was partially supported under NSF grants no. DMS-1700279, DMS-1400555, and NSF Postdoctoral Fellowship No. DMS-1400872.

Research of R.M. was partially supported under NSF grants no. DMS-0300487 and DMS-0801745, and by an IU COAS Dissertation Year Fellowship (2014-2015).

Research of K.Z. was partially supported under NSF grants no. DMS-1700279 and DMS-1400555.

1. 1 Introduction

2. 1.1 Problem and background

3. 1.2 Main results and outline of the paper

4. 1.3 Discussion and open problems

5. 1.4 Notation

6. 2 The $\beta$-model

7. 2.1 The inviscid $\beta$-model and its properties

8. 2.1.1 Hyperbolicity

9. 2.1.2 Weak involution and persistence of constraints

10. 2.1.3 Compatibility of constraints I: shock waves

11. 2.1.4 Compatibility with constraints II: normal modes

12. 2.1.5 Weak Lopatinsky stability

13. 2.1.6 Uniform Lopatinsky stability

14. 2.2 The viscous case

15. 2.2.1 Dissipativity

16. 2.2.2 Persistence of constraints

17. 2.2.3 Compatibility of shock profiles

18. 2.2.4 Compatibility of normal modes

19. 2.2.5 Weak and strong Evans-Lopatinski stability

20. 2.2.6 Uniform Evans-Lopatinsky stability and nonlinear stability

21. 3 Inviscid stability analysis: the Lopatinsky determinant

22. 3.1 Lopatinsky determinant: construction and asymptotic analysis

23. 3.1.1 Large-$h_{1}$ asymptotics

24. 3.1.2 Asymptotic instability

25. 3.1.3 Comparison with previous results: parallel case

26. 3.2 Full inviscid stability diagram

27. 3.2.1 The critical destabilization parameter.

28. 3.3 Remarks on the nonparallel case

29. 4 Viscous stability analysis: Evans function

30. 4.1 Computing the profile

31. 4.2 Computation of decaying manifolds and eigenfunctions

32. 4.3 Viscous stability diagram

33. 4.4 On the symmetry of eigenfunctions and equivariance of the Evans function

34. 4.5 Finding the critical destabilization parameter $h_{1}$

35. 4.6 Verifying concavity

36. 4.6.1 Finer points: analyticity and glancing

37. 5 Additional description of numerics

38. A On winding number computations

39. B Guide to [26] and [15] # 1: persistence of constraint condition, a second proof

40. C Guide to [26] and [15] # 2: a different approach towards the derivation of Rankine-Hugoniot conditions

41. C.1 The dynamic jump condition of [26]

1. Introduction

In this paper, continuing and extending investigations of [2, 6, 26, 46, 51, 52], we study by a combination of analytical and numerical Evans function techniques multi-D viscous and inviscid stability and associated transverse bifurcation of planar viscous slow Lax magnetohydrodynamic (MHD) shocks in a channel with periodic boundary conditions. Notably, this includes the first multi-D numerical Evans function study for viscous MHD, a computationally intensive problem representing the current state of the art, and, together with the treatment of gas-dynamical shocks in [33], the first such study for viscous shock profiles of any physical system in multi-D.

We obtain also new detail on the inviscid stability problem, while at the same time unifying and somewhat simplifying previously obtained results. In particular, we give a general framework for the treatment of constraints, or involutions [18], such as arise in multi-D MHD or elasticity, recovering and extending to the viscous case the fundamental results obtained by Blokhin et al (see [15, 26] and references therein) for inviscid MHD in a way apparently special to that case. The latter result answers in the affirmative the fundamental open problem posed in [45] whether Evans function stability, in the sense defined there, is necessary as well as sufficient for nonlinear multi-dimensional stability of viscous shocks in the presence of a constraint.

Our main physical conclusions are two. First, we make a mathematical connection between the spectral instability observed for slow inviscid MHD shocks [15, 26] and “corrugation instabilities” observed in the astrophysics community [38, 50], via a viscous bifurcation analysis as in [46, 51, 52]. Namely, we demonstrate transitions from stability to instability satisfying the bifurcation hypotheses proposed in [46] (parallel magnetic field case) and [52] (nonparallel case), implying bifurcation in a mode transverse to the direction of shock propagation, i.e., lying in the direction parallel to the front, the first examples for which these scenarios have been shown to occur. Our results suggest that, rather than a planar shock, a nonplanar “wrinkled” or “corrugated” traveling wave with nearby normal velocity is the typical mode of propagation in the slow Lax mode.

Second, continuing 1-D investigations of [2], we show numerically that the transition to instability for viscous multi-D slow Lax MHD shocks coincides with the transition to instability observed in the inviscid case. As shown in [60, 56, 57], for a rather general class of physical systems generalizing the “Kawashima class” of [41], viscous stability implies inviscid stability: that is, viscous effects may destabilize, but never stabilize a planar shock wave. The question posed in [60, 56, 57] whether and under what circumstances these two conditions coincide is a fundamental open problem in the theory of shock waves. On the one hand, it is much simpler to determine inviscid as opposed to viscous stability, so that pre-knowledge that they coincide would be a great help in applications; on the other, destabilization due to viscous effects would be extremely interesting physically. Our results here give the first information in this direction for multi-D, and (along with [33]) a first set of data for multi-D viscous systems. They are obtained by the introduction of an algorithm for numerical determination of the “refined stability condition” of [60, 57, 58, 12], detecting concavity at transition of the associated “critical” spectral curve through the origin, a new tool of general use. Our results, as in the 1-D study [2], suggest that viscous and inviscid stability transitions coincide.

For simplicity, our investigations here are confined to the 2-D case. However, the methods used apply in general dimensions. In the parallel case that is our main object of study, by rotational symmetry about the longitudinal axis, or direction of propagation of the front, multi-D stability- more generally, multi-D spectrum- reduces to the 2-D case, so there is no loss of generality in restricting to 2-D (cf. discussion of [26]). For the nonparallel case, discussed briefly here, 2-D stability is necessary but not sufficient for multi-D stability in dimensions $d\geq 3$. Likewise, there could conceivably be an earlier bifurcation involving spectral modes associated with the neglected third dimension. This, and the nonparallel case in general, would be interesting for further study.

1.1. Problem and background

The Navier–Stokes, or viscous, equations for isentropic 2-D MHD are given, in vectorial notation, by

[TABLE]

where $u=(u_{1},u_{2})$ is the velocity field, $h=(h_{1},h_{2})$ is the magnetic field, and $q=p+\frac{|h|^{2}}{2}$, where $p=p(\rho)$ is gas-dynamical pressure [9, 39, 17, 41, 19]. Here, $(x_{1},x_{2})$ is spatial location and $t$ is time. We take $p(\rho)=a\rho^{\gamma}$ corresponding to a $\gamma$-law, or polytropic equation of state. The corresponding Euler, or inviscid, equations are given by Eqs. (1.1) with $\mu=\eta=\nu=0$. In either (viscous or inviscid) case, the magnetic field must satisfy in addition the constraint

[TABLE]

which if satisfied at initial time $t=0$, may be seen to persist for all $t>0$.

Our aim is to study the spectral stability of both viscous and inviscid planar (without loss of generality standing) shock waves $u(x,t)\equiv\overline{u}(x_{1})$ in dimension $2$, either as solutions on the whole space, or - which amounts to restricting Fourier modes to a (discrete) lattice - as solutions on a two-dimensional channel, $x_{1}\in\mathbb{R},x_{2}\in[0,1]$, with periodic boundary conditions in $x_{2}$. In the parallel case $u=(u_{1},0)$, $h=(h_{1},0)$, Eqs. (1.1) decouple into the equations of nonmagnetic isentropic gas dynamics in $(\rho,u)$ and a heat equation for $h$, from which we may readily deduce that the set of parallel planar MHD shocks consists precisely of the set of nonmagnetic gas-dynamical shocks in $(\rho,u)$, adjoined with $h_{1}\equiv\text{\rm constant}$ (for details, see for instance [6, 45, 26]). As existence/transversality of traveling wave profiles for viscous polytropic gas dynamics is well known both in the isentropic [7] and nonisentropic [29, 55] case, one obtains thereby immediately existence/transversality of parallel MHD profiles, and, by perturbation, of near-parallel profiles as well; see [6] for details.

Shocks in MHD have different types defined by the number of characteristics at plus and minus infinity moving inward toward the shock.111A standard detail suppressed here is that the equations must first be recast in a form that is hyperbolic and noncharacteristic with respect to shock speed; see [45] or Sections 2 below. For parallel shocks, this is determined by the strength $|h_{1}|$ of the normal magnetic field $(h_{1},0)^{T}$, being “fast Lax” type for $0\leq|h_{1}|<H_{*}$, “intermediate” type for $H_{*}<|h_{1}|<H^{*}$, and “slow Lax” type for $H^{*}<|h_{1}|$ [26, 45], where

[TABLE]

(for $0<u_{1}^{+}<u_{1}^{-}$; see Lemma 3.1(iii)). Fast shocks are somewhat analogous to gas-dynamical shocks, and indeed reduce to this case in the zero-magnetic field limit $|h|\equiv 0$. Intermediate shocks are of nonclassical “overcompressive” type not appearing in gas dynamics [24, 56, 57, 6]. Slow shocks are of classical Lax type, but separated in parameter space from the fast type and exhibiting somewhat different properties. The patterns of incoming characteristics at plus and minus infinity for each type are displayed in Table 1 for the case of parallel shocks.

Inviscid numerical studies [21, 53] indicate that fast parallel shocks are typically stable, while slow parallel shocks are typically unstable. Intermediate shocks, since overcompressive, are always inviscid unstable, and will not be discussed here. (Nonetheless, they appear to play an important role in viscous behavior [24, 60, 58] where they have been seen numerically to be at least 1-D stable [6]). Indeed, it has been shown analytically [14, 27] that fast parallel MHD shocks are stable under the gas-dynamical stability condition of Majda [44], notably for a polytropic equation of state; likewise, it has been shown analytically [13, 15, 26] that slow MHD shocks are unstable in the infinite-magnetic field limit $|h|\to\infty$. In particular, in the brief but suggestive paper [26] Freistühler and Trakhinin, among other things, showed analytically the inviscid instability of slow Lax shocks for parallel MHD for sufficiently large magnetic field, extending to the parallel case (degenerate in this context [26]) the fundamental results of Blokhin et al [13, 15].

The result of Freistühler and Trakhinin [26] corroborates and puts on more solid mathematical ground an earlier study on instability of slow planar shocks in MHD performed by Stone and Edelman in the setting of astrophysics [50], where it is thought to play a role for example in the dynamics of accretion disks of binary dwarf stars. According to [50], the loss of stability through oscillations in the slow magnetosonic shock front is known as corrugation instability [50]. Though most of Stone and Edelman’s results rely on formal linear stability analysis, they also study the phenomena numerically through a time evolution code. The latter experiments suggest that the observed oscillatory linear instabilities result at nonlinear level in “fingers” that end up destroying the planar structure of the shock front. That is, the numerical results [50, Section 3.1 and 3.3] of Stone and Edelman indicate, further, that onset of instability is associated with loss of planar structure of the viscous profile, i.e., appearance of the above-mentioned corrugations.

The possibility of this latter phenomenon has been verified rigorously in the form of a steady transverse bifurcation in a general $\mathcal{O}$(2)-symmetric strictly parabolic system of conservation laws [46] relevant to the parallel MHD case, under appropriate spectral bifurcation conditions, namely, that transition to instability occurs through a pair of real eigenvalues corresponding to nonzero transverse Fourier modes passing through the origin. In the non-$\mathcal{O}$(2) symmetric case, corresponding to nonparallel magnetic field, a similar Hopf bifurcation result was shown in [52], under the assumption that loss of stability occurs through the passage of a complex conjugate pair of eigenvalues associated with nonzero transverse modes. Our ultimate goal is to verify these spectral scenarios by a detailed numerical study of the eigenmodes of the linearized operator about the shock.

1.2. Main results and outline of the paper

The first logical step in this work consists of combining the analytical conclusions of [26] of inviscid instability in the infinite-magnetic field limit with numerical Evans function results showing that slow Lax shocks can be stable for smaller magnetic fields. Putting these observations together, one may conclude the existence of a stability transition, associated with which one might hope to observe bifurcation in wave structure. This is far from obvious at the inviscid level, where such transitions are associated with infinitely many Fourier modes simultaneously entering the right half of the complex plane (see [10]); nor is it clear a priori that there is a corresponding stability transition at the viscous level, since viscous and inviscid spectra can differ greatly at mid and high frequencies. Nevertheless, by the result of Zumbrun and Serre [60, Proposition 5.3] connecting viscous and inviscid spectra in the low frequency regime, one may conjecture the associated appearance of more standard bifurcations in the better-behaved viscous case, involving finitely many low-frequency modes; see Section 1.3 or [59] for further discussion. In the simplest situation that the single (necessarily real) double eigenvalue pair (double by $\mathcal{O}$(2) symmetry) for large magnetic field moves into the stable half plane as magnetic field is decreased, without meeting any other eigenvalues along the way, this would necessarily be a “steady” spectral bifurcation, passing through $\lambda=0$. This simple scenario is likewise not a priori guaranteed, but our numerical investigations confirm that it is indeed what occurs.

To carry out these numerical investigations requires some interesting extensions of the standard Evans function framework [28, 56, 57, 58, 60] to handle the presence of constraints such as (1.2), similar to what was done for inviscid MHD by Blokhin et al [15, 26], and (partially) for viscous MHD by Métivier et al [45]. In the process, we unify and simplify these previous analyses, at the same time obtaining a new formulation of the MHD equations- the “$\beta$-model-” that is particularly convenient for numerics, combining in one model the desirable properties of noncharacteristicity, hyperbolicity, and conservation form; this is developed in Section 2.

To describe the main issues, a starting point is the observation that, in the presence of constraint (1.2), the equations of MHD are not prescribed uniquely, but only up to the addition of multiples of the constraint. Indeed, in the study of spectral stability, one could adjoin the constraint as an additional equation if desired. Thus, one must take care to choose a form of the equations possessing properties under which the standard Evans function and Lopatinsky determinants used to study viscous and inviscid stability are well defined, namely noncharacteristicity, hyperbolicity, and conservative form. There is a standard reformulation of the equations in which they become symmetric hyperbolic-parabolic and noncharacteristic [15, 17, 45] and another, different, formulation in which they become conservative; the standard approach has been to use ad hoc combinations of these in the analysis, depending on the need at hand. Here, we introduce for our single formulation a different analytical framework encompassing both viscous and inviscid cases. This gives necessity and sufficiency of the Evans condition for viscous MHD stability in the presence of constraint (1.2), answering an open problem posed in [45], where sufficiency but not necessity was established. In passing, we rederive and further illuminate the inviscid results of [15, 26]. These issues are discussed in Section 2, where the $\beta$-model and basic analytic framework are introduced.

In Section 3, we provide a description of the Rankine-Hugoniot conditions and the Lopatinsky determinant condition, giving the foundations for a careful study of inviscid instabilities. In particular, we (i) recapitulate in the more convenient $\beta$-model framework the large-magnetic field asymptotics of [26] showing instability, at the same time correcting certain computation errors in [26]; and (ii) carry out a numerical Lopatinsky study both verifying our asymptotics and extending the analysis to the small-magnetic field regime (see Figs. 3 and 2(b), respectively).

Two notable conclusions are that:

1. (i)

the large-magnetic field asymptotics are quite accurate, extending even down to rather small magnetic field strengths (see Fig. 1); 2. (ii)

there do exist inviscid stable slow parallel MHD shocks for sufficiently small magnetic field and sufficiently small amplitude of the shock.

The latter conclusion has the important implication that a stability transition, with potential for bifurcation, occurs. Finally, in Section 3.1.3, we compare our analytical and numerical multi-D results with those in [26] and [50]. Notably, we find that our large-magnetic field asymptotics improve by as much as $20\%$ on the accuracy of previous analyses. Specifically, as shown in Fig. 1, when compared with results in [26, Equation (61)], the calculation (3.10) of $\lambda_{2}$ shows a better agreement with numerical results; see also Section 3. This discrepancy is partly explained by the computations of Appendix C, where, among other things, we show that the dynamic Rankine-Hugoniot conditions of [26, Equation (44)] are incorrectly stated. This error, acknowledged by the authors in [54], leads to quantitative but not qualitative changes in the asymptotic results obtained in [26]. More important, the corrected version yields a result that is valid for all $\gamma\geq 1$ for isentropic-law, extending the result of [26, §3.4] constrained to $\gamma\in[1,2]$.

In Section 4, we connect the inviscid analysis with the viscous theory of stability of planar shocks based on the work of [60], through a study of the low-frequency limit (Table 3 and Fig. 5); this gives an additional check on our Lopatinsky computations through their asymptotic agreement with our numerical computations of the Evans function, an object defined in a completely different way.

We go on to carry out a complete numerical Evans function study of viscous stability over all parameters and frequencies (see Fig. 2(a)), verifying that there occurs the same stability transition, at the same parameter values, that were seen in the inviscid case. We present, further, numerical results verifying the bifurcation conditions assumed in the abstract results in [52] and [46] (see also the related [47]), consisting of the absence of other neutrally stable eigenvalues (see Fig. 6) and nonzero speed of crossing of the imaginary axis as the magnetic field is varied (see Fig. 8). We compute also approximate zero-eigenfunctions at the bifurcation point, yielding the approximate shape of the bifurcating nonplanar wave (see Fig. 4).

We finalize the paper with an appendix: in A we briefly explain how numerical winding number computations were carried out; in B, following [15, 26], we present an alternative proof by direct computation of Proposition 2.7 in the special case of inviscid MHD, showing that the constraint \mathrm{div}(h)\Big{|}_{t=0}=0 persists throughout the dynamics, namely, $\mathrm{div}\left(h(t)\right)\equiv 0$, for all $t\geq 0$; last, in C we give another perspective on the dynamic Rankine-Hugoniot condition of [15, 26], at the same time correcting details of some related calculations in [26].

1.3. Discussion and open problems

Our general results on Evans functions for systems with constraints pave the way for a unified treatment of multi-dimensional viscous shock stability in MHD, elasticity, and related equations arising in continuum mechanics. Our introduction of the $\beta$-model, though more special to MHD, by putting the equations into a standard symmetrizable conservative form, has the tremendous advantage that it allows computations using existing “off-the-shelf” code in the numerical stability package STABLAB [5]. As discussed in [33, 3, 4], success or failure of multi-dimensional computations is highly dependent on the specific algorithm used, with a number of catastrophic possible pitfalls that must be avoided. Thus, the ability to use existing, already-tested algorithms is of significant practical advantage.

We point out some further background and implications from a more general perspective contrasting viscous and inviscid stability. As shown in [60], viscous stability is closely related at low frequencies to inviscid stability, hence uniform viscous stability implies uniform inviscid stability. This means that, as shock or magnetic field strength is increased from a stable regime, a transition to inviscid instability implies a corresponding transition to viscous instability, occurring in low-frequency modes. The reverse is not true, as it is possible that a transition to viscous instability could occur in advance of the transition to inviscid instability due to destabilization of an intermediate- or high-frequency mode unrelated to the inviscid problem. Indeed, let $\xi\in\mathbb{R}$ be the Fourier frequency parameter in the direction $x_{2}$ transversal to the shock front. Then, as depicted in Fig. 7, there are essentially 3 different scenarios for viscous vs. inviscid stability transitions in a finite-cross section channel, depending mainly on concavity vs. convexity of the neutral spectral curve $\lambda=\lambda(\xi)$ for the viscous case, bifurcating from $\lambda(0)=0$, given by the sign of $\lambda^{\prime\prime}(\cdot)$ near $\xi=0$.

Recall [60], that this curve is tangent at $(\xi,\lambda)=(0,0)$ to the corresponding inviscid stability curve, given by homogeneity by a ray through the origin. The viscous spectral curve $\lambda(\cdot)$ is depicted in various cases in Fig. 7 along with its tangent inviscid ray. In the first case, $\lambda(\cdot)$ is concave ($\lambda^{\prime\prime}(\cdot)<0$) and we see that the transition to instability occurs simultaneously in the whole space (for which $\xi\in\mathbb{R}$) for viscous and inviscid problems, and slightly later for a duct of finite cross-section (for which $\xi\in{\mathbb{Z}}$, hence viscous spectra lags behind inviscid by a fixed finite amount). In the second, $\lambda(\cdot)$ is convex ($\lambda^{\prime\prime}(\cdot)>0$) and viscous instability occurs slightly before inviscid instability as the bifurcation parameter is varied, for either the whole space or finite cross-section. In the third, $\lambda(\cdot)$ is concave, but high-frequency instabilities cause the viscous problem to destabilize first. This is consistent with the results of [60] in the whole space, where it is shown that the viscous transition occurs not later than the inviscid one. For the whole space problem ($\xi\in\mathbb{R}$), our discussion above shows that simultaneously precisely in case 1, and strictly sooner in cases 2 and 3.

Our numerical results (see Section 4.5) for typical parameters $\gamma=5/3$, $u_{1}^{+}=0.86$ indicate that the viscous transition occurs at $H^{*}\approx 2$, while the inviscid transition occurs at $H^{*}\approx 1.995$. Likewise, one can see from Fig. 5 that the second derivative of the spectral parameter $\lambda=\lambda(\xi)$ with respect to $\xi$ is negative. Indeed, we find this to be the case for all parameter values; see Section 4.6 and particularly Fig. 11(a). That is, we appear to be in the first case depicted in Fig. 7. This has the important consequence that, considered in the whole space, viscous and inviscid stability transitions coincide for the 2-D transverse instabilities considered here, similarly as was seen in [2] for the 1-D longitudinal instabilities considered there. Thus, though viscosity can in principle according to the results of [60] hasten the onset of instability, for the two studies carried out so far (in [2] and here) for gas dynamics and MHD, this possibility was not in practice observed. It remains a very interesting open problem whether such “viscosity-enhanced instability” can occur for physically relevant models of gas dynamics or MHD, both philosophically, and- since inviscid transitions may often be computed explicitly, whereas viscous transitions require substantial numerical computation- from a practical point of view. Indeed, we note that it has been shown in [23] that the neutral inviscid stability curve, where the Lopatinsky condition precisely vanishes for $\lambda=0$, $\xi\neq 0$, may be determined explicitly, making this a practical condition indeed.

Our numerical study confirms not only stability transition but also that the spectral bifurcation assumptions assumed in [46, 51, 52] indeed hold: the loss of stability happens through a double eigenvalue crossing the imaginary axis. For fully parabolic “artificial viscosity” versions of MHD, the actual nonlinear bifurcation would follow from the results presented in [46, 51, 52] restricted to the space of divergence free magnetic field functions, provided that the corresponding spectral scenario were verified. Extending this result to the physical, “real” viscosity case considered here is an important open problem. Presumably, one could expect a similar spectral bifurcation for the artificial viscosity case (giving the full nonlinear result), but we do not investigate this here. The extension of our investigations to a complete “all-parameters” study of MHD shocks analogous to that done for gas dynamics in [33] is another important direction for further study: likewise, spectral stability of small amplitude nonextreme shocks as considered here, both inviscid and viscous.222 For extreme shocks, i.e., $1$- or $5$-shocks in the artificial viscosity case, see [25, 22].

It is worthwhile to emphasize the loss of planar structure observed here (see Section 4), in the numerical experiments of Edelman and Stone (cf. [50, Section 3.1 and 3.3]) and in the context of steady bifurcations in a $\mathcal{O}$(2)-symmetric strictly parabolic system of conservation laws; see [46]. This indicates that the resolution of the 2-D-Riemann problem for slow shocks in MHD is not realized through planar shocks, but generically involve nonplanar “corrugated” fronts as component slow waves; see also the recent numerical results of [20]. Spectral and nonlinear stability of these nonplanar waves is another very interesting open problem for further investigation.

Finally, we mention an interesting related analysis of Freistühler, Kleber, and Schropp [23] for the inviscid isothermal ($\gamma=1$) case, in which they find by explicit computation the inviscid stability boundary $\Delta(0,\pm 1)=0$ (in the notation of Section 3.1) for slow parallel shocks; a similar computation should be possible for general $\gamma$, sharpening our description of the inviscid boundary in Fig. 2-(b).333Namely, the dashed line in the figure, here computed numerically as described in Section 3.2.1. Freistühler et al investigate numerically also the fast shock case, showing that parallel isothermal fast shocks are uniformly stable, but nonparallel ones experience a transition to instability across a particular parameter surface. Similarly, for the general isentropic case ($\gamma>1$), Trakhinin [53] has shown that fast nonparallel shocks can be unstable in some regimes. A very interesting further study would be to carry out the corresponding analysis of viscous stability transition in these cases for fast shocks as we have done here for slow shocks.

1.4. Notation

In this paper we denote the real part (resp., imaginary part) of a number $z\in\mathbb{C}$ by $\mathrm{Re}(z)$ (resp., $\mathrm{Im}(z)$). We let $u\in\mathbb{R}^{n}$ be a vector of states assuming values $u^{\pm}$ accross a shock. Given any function $u\mapsto f(u)$ we write $[f(u)]=f(u^{+})-f(u^{-})$. The partial derivative of a differentiable mapping $(x_{1},x_{2},\ldots,x_{n})\mapsto g(x_{1},x_{2},\ldots,x_{n})$ with respect to the variable $x_{i}$ is written as $g_{x_{i}}.$ Given a matrix $A\in\mathbb{C}^{n\times n}$, we write $A^{-1}$ to denote its inverse, and $\sigma\left(A\right)$ to represent its eigenvalues.

Acknowledgments

R.M. would like to thank A. Mailybaev, D. Marchesin and C. Rohde for interesting conversations. K.Z. thanks David Lannes, Jinghua Yao, and Alin Pogan for interesting conversations regarding constraints. Numerical computations in this paper were made with the use of STABLAB, a MATLAB-based stability platform developed by B.B and K.Z. together with Jeffrey Humpherys and Joshua Lytle; SYMPY; and MATHEMATICA.

2. The $\beta$-model

Focusing our attention without loss of generality (by Galillean invariance of (1.1)) on shocks propagating in the $x_{1}$ direction, we now introduce the $\beta$-model, obtained by adding a constant multiple of $\mathrm{div}(h)$ to the equation (1.1c) in (1.1), yielding in place of equations (1.1) the system

[TABLE]

where $\beta$ is a fixed real-valued parameter and $e_{1}=(1,0)^{t}$. Recalling the constraint $\mathrm{div}(h)=0$ we see that (1.1) and (2.1) are equivalent for smooth solutions. However, (2.1) has certain practical advantages, as we now explain.

An important reason to modify (1.1) is that, considered without the constraint $\mathrm{div}(h)=0$, the inviscid version $\mu=\eta=\nu=0$, of (1.1) is not hyperbolic. A standard resolution of this problem is to substitute in (1.1) relations $\nabla(|h|^{2}/2)-\mathrm{div}(h\otimes h)=h\times(\nabla\times h)-h\mathrm{div}(h)=h\times(\nabla\times h)$ and

[TABLE]

which amounts to adding the nonconstant multiples $h\mathrm{div}(h)$ and $u\,\mathrm{div}(h)$ of the constraint $\mathrm{div}(h)$ to (1.1b) and (1.1c), to obtain the symmetrizable hyperbolic system [26, 15, 45]

[TABLE]

This has the advantage of both recovering hyperbolicity, and providing a symmetrizable system of equations. On the other hand, this handy device results in loss of conservative form of the equations, as a result of which jump conditions across shocks are not computable for this version of the equations. Meanwhile, the jump conditions for the conservative version (1.1) are degenerate (not full rank), and must be supplemented with that of the constraint in order to obtain the correct number of boundary conditions for the shock problem. See [45, §7] for further discussion. We refer elsewhere to this model, used in [14, 45, 26], as the “hybrid” or “standard” model.

It is clear that the inviscid $\beta$-model preserves the conservative structure of the equations; however, it is not symmetrizable. Nonetheless, as we will see shortly, it maintains (weak) hyperbolicity, among other useful properties, in particular those needed for linearized stability analysis by the study of Majda’s Lopatinsky determinant [44, 60, 45]. Moreover, it may be used “as is” for both interior equations and jump conditions, without additional modifications, allowing the use of standard numerical schemes for investigation of Lopatinsky stability. Likewise, the viscous version of the $\beta$-model is conservative and, though not of classical symmetrizable “Kawashima” form [41], exhibits the same favorable type of dispersion relation enjoyed by models of that type, among other properties needed for linearized stability analysis by the study of the Evans function [1, 28, 56, 57, 58, 60, 45].

2.1. The inviscid $\beta$-model and its properties

To explore the properties of the $\beta$-model we find it convenient to work in an abstract setting, following the discussion in [19, §5.4], afterwards specializing to our model. We start with the more complicated inviscid case, then finish by indicating briefly the treatment of the viscous case following [47, pp. 2 and 62-64].

System (2.1), in the inviscid case $\mu=\eta=\nu$, may be written in the general form

[TABLE]

with $d=2$, $n=5$, $\mathcal{W}=(\rho,u_{1},u_{2},h_{1},h_{2})^{T}$, and

[TABLE]

The linearization of (2.3) about a constant solution $\mathcal{W}\equiv\overline{\mathcal{\bar{W}}}$ is given by

[TABLE]

Taking $\mathcal{V}=f_{0}\left(\mathcal{W}\right)$, where $f_{0}(\cdot)$ is assumed to be a (local) diffeomorphism from $\mathbb{R}^{n}$ to itself, we can rewrite (2.3) in the standard form

[TABLE]

[TABLE]

with linearization about a constant solution $\mathcal{V}\equiv\overline{\mathcal{V}}=f_{0}(\overline{\mathcal{W}})$ given by

[TABLE]

Likewise, the divergence-free constraint may be written in the general form

[TABLE]

satisfying Dafermos’ compatibility condition (involution) (cf. [18, Chap. V, Section 5.4]):

[TABLE]

where $\mathcal{M}$ is first-order constant-coefficient:

[TABLE]

and hyperbolic:

[TABLE]

Note that (2.10) implies also the linearized version:

[TABLE]

where $\mathscr{L}:=-\sum_{j=1}^{d}\mathcal{A}_{j}^{\pm}\partial_{x_{j}}$, $x_{1}\gtrless 0$, is the linearized operator about a background constant solution.

Remark 2.1*.*

For (2.1) we have $d=2$, $m=1$, $n=5$, $\mathcal{V}=\left(\rho,\rho u_{1},\rho u_{2},h_{1},h_{2}\right)$, $\Gamma_{1}=(0,0,0,1,0)$, $\Gamma_{2}=(0,0,0,0,1)$, $\mathcal{M}_{1}=\beta$, and $M_{2}=0$. Thus, $\Gamma\mathcal{V}=\mathrm{div}(h)$ and $\mathcal{M}=\beta\partial_{x_{1}}$.

2.1.1. Hyperbolicity

For a given reference state $\overline{\mathcal{V}}$ and direction vector $k\in\mathbb{R}^{d}$, define

[TABLE]

From (2.13), we have

[TABLE]

whence $\ker\Gamma(k)$ is a right invariant subspace and ${\rm Range}\,\Gamma(k)^{T}$ a left invariant subspace of $A(k)$. Likewise, $\hbox{\rm Coker}\Gamma(k)$ is a left invariant subspace and ${\rm Range}\,\Gamma(k)$ a right invariant subspace of $\mathcal{M}(k)$.

Hyperbolicity of (2.6) with constraint (2.9) is defined as the property that the eigenvalues $\sigma(A(k)|_{\ker\Gamma(k)})$ of $A(k)|_{\ker\Gamma(k)}$ be real and semisimple for all $k\in\mathbb{R}^{d}$, and “weak hyperbolicity” of (2.6) with (2.9) as the property that $\sigma(A(k)|_{\ker\Gamma(k)})$ be real for all $k\in\mathbb{R}^{d}$. Hyperbolicity of (2.6) without constraint (2.9) is defined as the property that $\sigma(A(k))$ be real and semisimple for all $k\in\mathbb{R}^{d}$, and “weak hyperbolicity” without (2.9) as the property that $\sigma(A(k))$ be real for all $k\in\mathbb{R}^{d}$.

Proposition 2.2**.**

Under the above assumptions, it holds that

1. (i)

System (2.6) is weakly hyperbolic with constraint (2.9) if and only if it is weakly hyperbolic without the constraint. In particular, $\mu$, $r$ are an eigenvalue, eigenvector pair of $\mathscr{L}(k)$ if and only if either (a) $\Gamma(k)r=0$, or (b) $\gamma:=\Gamma(k)r$ is an eigenvector of $\mathcal{M}(k)$ with eigenvalue $\mu$, $0\neq k\in\mathbb{R}^{d}$, whence

[TABLE] 2. (ii)

The equation (2.6) is hyperbolic with the constraint if it is hyperbolic without the constraint. If (2.6) is hyperbolic with the constraint, then for every $k$ for which $\sigma(\mathcal{M}(k)|_{{\rm Range}\,\Gamma(k)})\cap\sigma(A(k)|_{\ker\Gamma(k)})=\emptyset$, $\sigma(A(k))$ is imaginary and semisimple. 3. (iii)

If (2.6) is hyperbolic with the constraint, and $\sigma(A(k))$ is semisimple for all $k\in\mathbb{R}^{d}$, then (2.6) is hyperbolic without the constraint. However, in general hyperbolicity with constraint (2.9) does not imply hyperbolicity without the constraint.

Proof.

Applying (2.14) to $(A(k)-\mu)r=0$, we obtain $(\mathcal{M}(k)-\mu)\gamma=0$, with $\gamma=\Gamma(k)r$, and thus the second assertion in (i), dichotomy (a)-(b). In case (a), $r\in\ker\Gamma(k)$ and so $\mu\in\sigma(A(k)|_{\ker\Gamma(k)})$; in case (b), $\mu\in\sigma(\mathcal{M}(k)|_{{\rm Range}\,\Gamma(k)})$. Combining, we obtain (2.15). The first assertion in (i) then follows by hyperbolicity of $\mathcal{M}$, which implies $\sigma(\mathcal{M}(k))$ real, hence $\sigma(\mathcal{M}(k)_{S})$ real on any invariant subspace of $\mathcal{M}(k)$.

The first assertion in (ii) follows by the fact that the eigenvalues and eigenvectors of $A(k)$ restricted to an invariant subset are a subset of the eigenvalues and eigenvectors of $A(k)$. The second assertion in (ii) holds because $\sigma(\mathcal{M}(k)|_{{\rm Range}\,\Gamma(k)})\cap\sigma(A(k)|_{\ker\Gamma(k)})=\emptyset$ implies that the number of eigenvectors and eigenvalues of $\mathcal{M}(k)|_{{\rm Range}\,\Gamma(k)}$ and $A(k)|_{\ker\Gamma(k)}$ is then $n$, together accounting for all eigenvalues of $A(k)$, along with the fact that $(M-\mu)\Gamma S=0$ implies $\Gamma(A-\mu)S=0$ and so $(A-\mu)S\in\ker\Gamma$- this, along with separation of $\sigma(M(k)|_{{\rm Range}\,\Gamma})$ and $\sigma(A(k)|_{\ker\Gamma})$ gives $(A-\mu)(S+R)=0$ for some $R\in\ker\Gamma$, yielding semisimplicity of eigenvalues lying in $\sigma(M(k)|_{{\rm Range}\,\Gamma})$.

The first assertion in (iii) follows by the observation that weak hyperbolicity plus semisimplicity of $\sigma(A(k))$ for all $k\in\mathbb{R}^{d}$ implies hyperbolicity. Finally, the example $\mathcal{A_{1}}=\begin{pmatrix}0&1\\ 0&0\end{pmatrix}$, $\Gamma_{1}=(0,1)$, $\ker\Gamma_{1}=(1,0)^{T}$, $\mathcal{M}_{1}=0$ for $d=1$, $n=2$, and $m=1$ satisfies hyperbolicity with constraint (since $A|_{\ker\Gamma_{1}}=0$) and hyperbolicity of $\mathcal{M}$, but is not hyperbolic without constraint, demonstrating the second assertion in (iii). ∎

Remark 2.3*.*

Assuming hyperbolicity of $\mathcal{A}$ with constraint $\Gamma$ and hyperbolicity of $\mathcal{M}$, consider the generalized eigenvector equation $(\mathcal{A}(k)-\mu)^{i}r=0$. Applying $\Gamma(k)$ on the left, we obtain $(\mathcal{M}(k)-\mu)^{i}\Gamma(k)r=0$, giving $(\mathcal{M}(k)-\mu)\Gamma(k)r=0$ by hyperbolicity of $\mathcal{M}$. But, this in turn gives $0=\Gamma(\mathcal{A}(k)-\mu)r$, or $(\mathcal{A}(k)-\mu)r\in\ker\Gamma(k)$, from which we may determine that $r$ is a generalized eigenvector of $\mathcal{A}(k)$ of height $2$ over a genuine eigenvector in $\ker\Gamma(k)$. Thus, the counterexample given in the proof is indicative of the general case, the only obstruction to hyperbolicity being possible appearance of Jordan blocks of height $2$.

Proposition 2.4**.**

The $\beta$-model is weakly hyperbolic without constraint (2.9) and hyperbolic with the constraint. In the parallel case $u=(\bar{u}_{1},0)$, $h=(\bar{h}_{1},0)$, it is not hyperbolic without the constraint if $\beta\neq\bar{u}_{1}$.

Proof.

Hyperbolicity of the MHD equations with constraint (2.9) is well known; see, e.g., [45]. Likewise, $\mathcal{M}$ in this case, since scalar ($m=1$), is automatically hyperbolic. Thus, we obtain weak hyperbolicity without constraint by Proposition 2.2. On the other hand, taking without loss of generality $\bar{u}=0$, by Galillean invariance, direct computation using (2.4), (2.8)(ii) gives $A_{0}={\rm diag}\{1,\rho,\rho,1,1\}$, and

[TABLE]

hence, in the parallel case, for $k=(0,1)$, $\Gamma(k)=(0,0,0,0,1)$, $\mathcal{M}(k)=0$, and

[TABLE]

Noting that the left kernel ${\rm Span}\{(-\bar{h}_{1},0,0,\bar{\rho})\}$ of $\mathcal{A}(k)|_{\ker\Gamma(k)}$ is not orthogonal to the upper righthand block $(0,0,0,\beta)^{T}$ of $\mathcal{A}(k)$ unless $\beta=0$, we find that (2.1) is not hyperbolic if $\beta\neq 0$, or, in original coordinates, $\beta\neq\bar{u}_{1}$. ∎

The appearance of Jordan blocks with maximum height 2 suggests that one might establish local existence for (2.6) with loss of one derivative, as for example in the model (linear) case $u_{t}=v_{x}$, $v_{t}=0$. However, we mention this, and the details of Proposition 2.4 for general interest only; for our purposes (namely, analysis of linear shock stability), weak hyperbolicity is what is needed.

2.1.2. Weak involution and persistence of constraints

For $C^{2}$ solutions $\mathcal{V}$ of (2.6), it is easy to see that constraint (2.9) is preserved under the evolution of (2.6). For, applying $\Gamma$ to (2.6), we obtain $(\partial_{t}-\mathcal{M})(\Gamma\mathcal{V})=0$, whence, by well-posedness (hyperbolicity) of the constant-coefficient equation $(\partial_{t}-\mathcal{M})\mathcal{V}=0$, if $\Gamma\mathcal{V}\equiv 0$ at time $t=0$, then $\Gamma\mathcal{V}\equiv 0$ for all $t\geq 0$ as well. We show now that the same holds true for weak solutions of (2.6); indeed, (2.10) holds (weakly) also for discontinuous $\mathcal{V}$.

The key to these results is to observe that condition (2.10), relating a nonlinear expression on the left to a linear constant-coefficient one on the right, implies substantial structure. Decomposing

[TABLE]

where $\mathcal{Q}_{j}(\mathcal{V})=O(|\mathcal{V}|^{2})$ denotes the nonlinear part of $\mathcal{F}_{j}$, we have the following algebraic relations.

Lemma 2.5**.**

With $\mathcal{N}_{j}$, $\mathcal{Q}_{j}$ as above, involution property (2.10) is equivalent to

[TABLE]

Proof.

Relations (2.18)(i) may be deduced from (2.13) evaluated at $\mathcal{V}=0$, comparing coefficients of like derivatives of $\mathcal{V}$. Subtracting (2.18)(i) from (2.10) gives $\Gamma\sum_{j}\partial_{x_{j}}\mathcal{Q}(\mathcal{V})=0$, whence (2.18)(ii) follows by taking $\mathcal{V}(x)=\sum w_{j}x_{j}$ for arbitrary $w\in\mathbb{R}^{n}$ to find that

[TABLE]

for all $w_{j}$, $w_{k}$, whence each coefficient $d^{2}(\Gamma_{j}\mathcal{Q}_{k}+\Gamma_{k}\mathcal{Q}_{j})$ vanishes separately and so $d(\Gamma_{j}\mathcal{Q}_{k}+\Gamma_{k}\mathcal{Q}_{j})=\text{\rm constant}$, and $(\Gamma_{j}\mathcal{Q}_{k}+\Gamma_{k}\mathcal{Q}_{j})$ must be linear in $\mathcal{V}$. By $\mathcal{Q}_{j}(\mathcal{V})=O(|\mathcal{V}|^{2})$, this implies $(\Gamma_{j}\mathcal{Q}_{k}+\Gamma_{k}\mathcal{Q}_{j})=0$, or (2.18)(ii). ∎

Corollary 2.6**.**

For $\mathcal{V}\in L^{\infty}$, (2.10) holds in weak sense; indeed, for $\phi\in C^{2}$,

[TABLE]

Proof.

Relation (2.19) follows immediately from (2.18), comparing like $\phi$-derivatives term by term. Weak satisfaction of (2.10) then follows by integration in $x$ of (2.19) for $\phi\in C^{\infty}_{0}$. ∎

Corollary 2.7**.**

An $L^{\infty}$ weak solution $\mathcal{V}$ of (2.6) satisfying (2.9) weakly at time $t=0$ satisfies $(\partial_{t}-\mathcal{M})\Gamma\mathcal{V}=0$ and (2.9) weakly for all $t\geq 0$.

Proof.

If $\mathcal{V}$ is a weak solution of (2.6), then it is a weak solution of $\Gamma(\partial_{t}\mathcal{V}-\mathcal{F}(\mathcal{V})=0$, hence, by Corollary 2.6, a weak solution of $(\partial_{t}-\mathcal{M})\Gamma\mathcal{V}=0$. Next, we recall that a function $U\in L^{p}$ satisfies a constant-coefficient linear differential differential equation $LU=0$ in $z=(z_{1},z_{2})\in\mathbb{R}^{m+n}$ weakly if and only if $LU^{\varepsilon}=0$ for all $\varepsilon>0$, where $U^{\varepsilon}:=U*\eta^{\varepsilon}$ is the mollification of $U$ in coordinate $z_{1}$ by a symmetric smoothing kernel $\eta^{\varepsilon}(z_{1}):=\varepsilon^{-1/m}\eta(|z_{1}|/\varepsilon)$. For, denoting $\langle g,h\rangle_{*}=\int_{\mathbb{R}^{m}}g^{*}f(z)dz$, the relation $\langle L^{*}\phi,U^{\varepsilon}\rangle_{*}=\langle L^{*}\phi^{\varepsilon},U\rangle_{*}$ together with $U^{\varepsilon}\to U$ in $L^{1}$ as $\varepsilon\to 0^{+}$ gives the result. Letting $\mathcal{V}^{\varepsilon}(\cdot,t)\in C^{\infty}$ denote mollification in $x$, we thus have $\phi^{\varepsilon}:=\Gamma V^{\varepsilon}\equiv 0$ (strongly) at $t=0$ while $(\partial_{t}-\mathcal{M})\phi^{\varepsilon}=0$ weakly for $t\geq 0$. It follows by hyperbolicity/uniqueness of weak solutions that $\phi^{\varepsilon}\equiv 0$ and thus $\Gamma\mathcal{V}^{\varepsilon}=0$ for all $t\geq 0$. Taking $\varepsilon\to 0^{+}$, we find that $\Gamma\mathcal{V}=0$ weakly for all $t\geq 0$. ∎

Corollary 2.8**.**

For any $\beta\in\mathbb{R}$, an $L^{\infty}$ weak solution $\mathcal{V}$ of the $\beta$-model (2.1) satisfying ${\rm div}(h)=0$ weakly at $t=0$ satisfies (2.9) weakly for all $t\geq 0$.

Remark 2.9*.*

Corollary 2.8 includes the case $\beta=0$ for which the $\beta$-model coincides with the original, conservative form (1.1) of the MHD equations. For weak solutions with ${\rm div}(h)$ initially weakly vanishing, this justifies the manipulations used in deriving both the general $\beta$-model, and the hybrid symmetrizable hyperbolic/augmented jump condition system of [26, 15, 45], based on weak satisfaction of the constraint ${\rm div}(h)$. In particular, it shows that weak solutions of all three models coincide, where defined, for weakly divergence-free data. We note for the hybrid model that weak solutions are by its nature defined only for piecewise smooth solutions with entropic shock discontinuities. Within this class, persistence of the constraint is shown, e.g., in [45, Proposition 7.1], via weak satisfaction of

[TABLE]

the associated jump condition $[(s-u_{1}){\rm div}(h)]=0$ for the constraint at a shock traveling with speed $s$ being obtained from the jump condition for shocks of the augmented equation (2.20). For a linearized version of this result, see [26, Remark 3.2]. We offer a second proof of persistence for the linearized $\beta$-model in Appendix B, through an argument like that in [26, 45].

Remark 2.10* (Smooth vs. weak persistence).*

It is important to notice that our framework here is rather special, even in situations when constraints are preserved for smooth solutions: for example, when the operator $\mathcal{M}$ in (2.10) is quasilinear hyperbolic rather than linear, constant-coefficient. In this case, satisfaction of the constraint may be “broken” by appearance of shock discontinuities; that is, persistence of constraints could hold for smooth, but not general weak solutions. A very interesting example pointed out to us by D. Lannes is 2-D shallow water flow, for which vorticity can be created in initially irrotational flow through wave breaking/formation of shock waves [42].

2.1.3. Compatibility of constraints I: shock waves

Corollary (2.7) suggests that shock waves of (2.6) associated with characteristics of the constrained model (2.6), (2.9) should be compatible with the constraint in the sense that they satisfy $\Gamma\mathcal{V}$ weakly, since such shocks may be expected to form from data that is initially smooth and satisfying the constraint. That is, for a shock wave

[TABLE]

propagating in direction $k\in S^{d-1}$ with speed $s$, and satisfying Rankine-Hugoniot conditions

[TABLE]

where $[h]:=h(\mathcal{V}_{+})-h(\mathcal{V}_{-})$ denotes jump in $h$ across the shock, whose speed is not associated with the fixed, complementary set of characteristics $\sigma(\mathcal{M}(k)|_{{\rm Range}\,\Gamma})$ governing propagation of constraints (recall decomposition (2.15)), we expect there to hold also the jump conditions

[TABLE]

associated with (2.9). We now verify this intuition by direct computation using Lemma 2.5.

Let us first remark that small-amplitude shocks do certainly exist in simple characteristic fields $\alpha_{j}(k)\in\sigma(\mathcal{A}(k)|_{\ker\Gamma(k)})$ of the constrained model (2.6), (2.9), by the standard bifurcation analysis of Lax [43, 48], with jump $[\mathcal{V}]$ in $\mathcal{V}$ across the shock front lying approximately in the associated eigendirection $r(k)\in\ker\Gamma$, hence satisfying (2.23) to first order at least. Likewise, in the complementary case of simple characteristic fields $\alpha_{j}(k)\in\sigma(\mathcal{M}(k)|_{{\rm Range}\,\Gamma(k)})$ associated with propagation of $\Gamma\mathcal{V}$. However, in the latter case, it is easy to see that these are contact discontinuities with $[\mathcal{V}]$ lying approximately in the associated eigendirection $r(k)$, which, moreover, is specifically not annihilated by $\Gamma(k)$. Hence, for such discontinuities, (2.9) definitely does not hold weakly. For large-amplitude shocks, or multiple characteristics, shock existence must be studied on a case-by-case basis. For MHD, global existence is well known, as recorded for the parallel case in Lemma 3.1 below.

Recall finally [60] that shocks are classified as different types– e.g., “Lax,” “overcompressive,” or “undercompressive” – according to the difference between the numbers of incoming and outgoing characteristics at their endstates $\mathcal{V}_{\pm}$.

Proposition 2.11**.**

Assuming involution property (2.10), a shock solution (2.21) of (2.6) satisfying (2.22) with speed $s\not\in\sigma(\mathcal{M}(k)|_{{\rm Range}\,\Gamma})$ automatically satisfies (2.9) in the weak sense (2.23). Moreover, the type of the shock is the same considered with respect to characteristics of the system with or without constraint.

Proof.

Expanding $\mathcal{F}(k)=\mathcal{N}(k)+\mathcal{Q}(k)$, with ${N}(k):=\sum_{j=1}^{d}\mathcal{N}_{j}k_{j}$ and $\mathcal{Q}(k):=\sum_{j=1}^{d}\mathcal{Q}_{j}k_{j}$, and using (2.18), we find that $0=\Gamma(k)[\mathcal{F}(k)-s\mathcal{V}]=\Gamma(k)(\mathcal{N}(k)-s)[\mathcal{V}]=(\mathcal{M}(k)-s)[\Gamma(k)\mathcal{V}]$, whence (2.23) follows by the assumed invertibility of $\mathcal{M}(k)-s$ on ${\rm Range}\,\Gamma(k)$. Moreover, by assumption, the additional characteristics $\sigma(\mathcal{M}(k)|_{{\rm Range}\,\Gamma})$ present for the unconstrained system besides the characteristics $\sigma(\mathcal{A}|_{\ker\Gamma})$ of the constrained one are distinct from $s$, hence noncharacteristic. Moreover, since $\mathcal{M}(k)$ is independent of $\mathcal{V}_{\pm}$, the number of resulting new characteristics incoming to the shock is equal to the number of new characteristics outgoing from the shock, each incoming characteristic on one side being matched by an outgoing characteristic of the same sign relative to $s$ on the other side. As shock type is defined by the difference between the numbers of incoming and outgoing characteristics [60], it is therefore unaffected by presence or absence of the constraint. ∎

Corollary 2.12**.**

For stationary MHD shocks in the $x_{1}$ direction $\mathcal{W}(x,t)=\overline{\mathcal{W}}(x_{1})=\mathcal{W}_{\pm}$ for $x_{1}\gtrless 0$ of the $\beta$-model (2.1) and any $\beta\neq 0$, ${\rm div}h=0$ holds weakly, i.e., $[h_{1}]=0$ across the shock. Moreover, the type of the shock is the same considered with or without the constraint.

Proof.

For the $\beta$-model, we have (see Remark 2.1) $\mathcal{M}_{1}=\beta$, $\mathcal{M}_{2}=0$, hence $\sigma(\mathcal{M}(k))=\beta$ for shock speed $s=0$ and direction $k=(1,0)$, and so $s\not\in\sigma(\mathcal{M}(k))$ precisely when $\beta\neq 0$. ∎

Remark 2.13*.*

In the case $\beta=0$, corresponding to (1.1), the jump condition $[h_{1}]=0$ for ${\rm div}h=0$ is consistent with but not implied by the Rankine-Hugoniot conditions, which are in this case degenerate (not full rank) [15, 45, 26]. We note that any inviscid shock obtainable as the limit of smooth viscous shock profiles of (1.1) with $\mu,\eta,\nu\neq 0$ necessarily satisfies ${\rm div}h=0$, so this assumption is not only internally consistent (as just shown), but also consistent with the vanishing viscosity point of view.

2.1.4. Compatibility with constraints II: normal modes

For definiteness, restrict now (without loss of generality) to the case of a zero-speed planar shock in direction $x_{1}$, of standard Lax type, in particular satisfying

[TABLE]

Note that (2.24), by (2.15), automatically implies $0\not\in\sigma(\mathcal{M}_{1}|_{{\rm Range}\,\Gamma_{1}})$. Recall from, e.g., [60], the normal modes equations after shifting to a frame $(z,x_{2},\dots,x_{d})$ with discontinuity at $x_{1}=0$:

[TABLE]

(interior equation) and

[TABLE]

(linearized jump conditions), where $[\cdot]$ denotes jump across $x_{1}=0$. Here, $\widehat{\mathcal{V}}\in\mathbb{C}^{n}$ denotes the Laplace–Fourier transform of the coordinated-shifted solution $\mathcal{V}$ and $Y\in\mathbb{C}$ the Laplace–Fourier transform of the shift, or front location in the original spatial coordinates, with $\xi=(\xi_{2},\dots,\xi_{d})\in\mathbb{R}^{d-1}$ the Fourier frequency in directions $(x_{2},\dots,x_{d})$ and $\lambda\in\mathbb{C}$ the Laplace frequency in time $t$. Associated with each solution $\mathcal{V}\in L^{2}(\mathbb{R})$ of (2.25)-(2.26) is a normal mode

[TABLE]

of the linearized coordinate-shifted equations, with $\mathrm{Re}(\lambda)>0$ corresponding to linear instability.

Hereafter, we replace the shifted coordinate $z$ by its original designation $x_{1}$. Define symbols

[TABLE]

where $\widetilde{\mathcal{A}}(\xi)=\sum_{j=2}^{d}i\xi_{j}\mathcal{A}_{j}$, $\widetilde{\Gamma}(\xi)=\sum_{j=2}^{d}i\xi_{j}\Gamma_{j}$, and $\widetilde{\mathcal{M}}(\xi)=\sum_{j=2}^{d}i\xi_{j}\mathcal{M}_{j}$. Likewise, denote

[TABLE]

Lemma 2.14**.**

Assuming (2.10), (2.24), for $\widehat{\mathcal{V}}$ piecewise smooth and satisfying (2.25)-(2.26):

1. (i)

$(\lambda+\widehat{\mathcal{M}})\widehat{\Gamma}\widehat{\mathcal{V}}=0$* for $x_{1}\gtrless 0$.* 2. (ii)

$[\mathcal{M}_{1}\widehat{\Gamma}\widehat{\mathcal{V}}]=0$* at $x_{1}=0$, hence $(\lambda+\widehat{\mathcal{M}})\widehat{\Gamma}\widehat{\mathcal{V}}=0$ weakly on $\mathbb{R}$.*

Proof.

From (2.13), (2.18), we obtain the corresponding relations

[TABLE]

[TABLE]

Recall the jump condition (2.22) for the shock, specialized to the $x_{1}$-directional case $k=(1,0,\dots,0)$:

[TABLE]

and the jump condition (2.23) for the constraint:

[TABLE]

obtained by applying $\Gamma_{1}$ to (2.30) and using (2.18) and $0\not\in\sigma(\mathcal{M}_{1})$. Similarly, applying $\tilde{\Gamma}$ to (2.30) and using (2.29), we obtain $0=\tilde{\Gamma}[\mathcal{N}_{1}\overline{\mathcal{V}}+\mathcal{Q}_{1}(\overline{\mathcal{V}})]=-\Gamma_{1}[\widetilde{\mathcal{N}}\overline{\mathcal{V}}+\mathcal{Q}_{1}(\overline{\mathcal{V}})]+[\widetilde{\mathcal{M}}\Gamma_{1}+\mathcal{M}_{1}\tilde{\Gamma}],$ or

[TABLE]

Applying now $\Gamma_{1}$ to (2.26), and using (2.28)-(2.29) and (2.31), gives

[TABLE]

yielding, by invertibility of $\mathcal{M}_{1}$ (a consequence of (2.24), as noted above),

[TABLE]

Writing (2.25) as $(\lambda+\mathcal{A}_{1}\partial_{x_{1}}+\widetilde{\mathcal{A}})\widehat{\mathcal{V}}=0$, applying $\hat{\Gamma}$ on the left, and using (2.28), we obtain immediately $(\lambda+\widehat{\mathcal{M}})\hat{\Gamma}\widehat{\mathcal{V}}$, verifying (i). Using (2.25) to express $\partial_{x_{1}}\widehat{\mathcal{V}}=\mathcal{A}_{1}^{-1}(\lambda+\widetilde{\mathcal{A}})\widehat{\mathcal{V}},$ we obtain

[TABLE]

or, using $\mathcal{M}_{1}\Gamma_{1}=\Gamma_{1}\mathcal{A}_{1}$ to express $\mathcal{M}_{1}\Gamma_{1}\mathcal{A}_{1}^{-1}=\Gamma_{1}$, \mathcal{M}_{1}[\hat{\Gamma}\widehat{\mathcal{V}}]=\big{(}-\Gamma_{1}(\lambda+\widetilde{\mathcal{A}})+\mathcal{M}_{1}\widetilde{\Gamma}\big{)}[\widehat{\mathcal{V}}]. Using (2.28) to express $\mathcal{M}_{1}\tilde{\Gamma}-\Gamma_{1}\widetilde{\mathcal{A}}=\tilde{\Gamma}\mathcal{A}_{1}-\widetilde{\mathcal{M}}$, we obtain finally

[TABLE]

or, substituting for $[\Gamma_{1}\widehat{\mathcal{V}}]$ using (2.31), and then applying (2.29):

[TABLE]

which vanishes by (2.26). By invertibility of $\mathcal{M}_{1}$, we thus obtain $[\widehat{\Gamma}\widehat{\mathcal{V}}]=0$, verifying (ii). ∎

Corollary 2.15**.**

Assuming (2.10), (2.24), for $\widehat{\mathcal{V}}\in L^{2}$ piecewise smooth and satisfying (2.25)-(2.26), we have $\widehat{\Gamma}\widehat{\mathcal{V}}\equiv 0$. That is, decaying normal modes automatically satisfy the Laplace-Fourier transformed version of constraint (2.9).

Proof.

By Lemma 2.14, $\widehat{\varphi}:=\hat{\Gamma}\widehat{\mathcal{V}}$ is a weak solution of $(\lambda+\widehat{\mathcal{M}}(\xi))\widehat{\varphi}=0$. Mollifying $\widehat{\varphi}$ by convolution with a standard smoothing kernel $\eta^{\varepsilon}$, we obtain a family of $C^{\infty}\cap L^{2}$ solutions $\widehat{\varphi}^{\varepsilon}:=\widehat{\varphi}*\eta^{\varepsilon}$ of $(\lambda-\widehat{\mathcal{M}}(\xi))f=0$, converging in $L^{2}$ to $\varphi$ as $\varepsilon\to 0^{+}$. These solutions, if nontrivial, would represent eigenfunctions of $-\widehat{\mathcal{M}}(\xi)$ with eigenvalue $\lambda$. However, $\mathcal{M}$, being constant coefficient, has only continuous spectrum, and so each $\widehat{\varphi}^{\varepsilon}$ must vanish identically, as therefore does $\widehat{\varphi}$ in the limit as $\varepsilon\to 0$. ∎

Corollary 2.16**.**

For stationary MHD shocks in the $x_{1}$ direction of the $\beta$-model (2.1), $\beta\neq 0$, and $\widehat{\mathcal{V}}\in L^{2}$ piecewise smooth and satisfying (2.25)-(2.26), we have $\tilde{\Gamma}\widehat{\mathcal{V}}=\beta(\partial_{x_{1}}\widehat{\mathcal{V}}_{4}+i\xi\widehat{\mathcal{V}}_{5})\equiv 0$. That is, decaying normal modes automatically satisfy the (transformed) divergence-free constraint.

Remark 2.17*.*

Relation $(\lambda+\widetilde{\mathcal{M}}\tilde{\Gamma})\widetilde{\mathcal{V}}=0$ may be recognized as the Laplace-Fourier transformed version of $(\partial_{t}+\mathcal{M})\Gamma\mathcal{V}=0$ of Corollary 2.7. For the hybrid model of [14, 45, 26], weak satisfaction of the Fourier-Laplace transform analog

[TABLE]

of 2.20, Remark 2.9, pointed out in [26, (47)-(48), Remark 3.2], likewise gives the result that decaying normal modes satisfy the divergence-free constraint, a fundamental result first noted in [15, 14].

Remark 2.18*.*

For both the $\beta$-model and the standard (hybrid) MHD model of [14, 45, 26], the fact that the mode of propagation of the constraint is scalar convection/advection shows that not all decaying modes of $(\mathcal{A}_{1}\partial_{x_{1}}+(\lambda+i\xi\mathcal{A}_{2}))\widetilde{\mathcal{V}}=0$ satisfy the constraint $\tilde{\Gamma}\widetilde{\mathcal{V}}=0$. For, inverting the dispersion relation $\lambda(k_{1},\xi)=ik_{1}\bar{u}_{1}+i\xi\bar{u}_{2}$ to solve for the associated eigenvalue $\mu:=ik=(\lambda-i\xi\bar{u}_{2})/\bar{u}_{1}$ of $\mathcal{A}_{1}^{-1}(\lambda+i\xi\mathcal{A}_{2})$, we see for $\xi=0$ that $\mathrm{Re}(\mu)=\mathrm{Re}(\lambda)$, and so the mode carrying the constraint is decaying as $x\to-\infty$, but does not satisfy the constraint. Thus, the situation of Corollary 2.16 is more subtle than just checking that the constraint is satisfied for all relevant modes.

2.1.5. Weak Lopatinsky stability

Gathering information, we are now ready to make conclusions regarding linearized and nonlinear stability with and without constraints. We first recall the weak stability, or “weak Lopatinsky” condition of Majda [44], that (i) the dimensions of the subspaces $\mathcal{E}^{\pm}$ of decaying modes of (2.25) on $x\gtrless 0$ sum to $n-1$ (consistent splitting [1, 60]) and (ii) for $\mathrm{Re}(\lambda)>0$, there exist no decaying normal modes $\widehat{\mathcal{V}}\in L^{2}$ piecewise smooth and satisfying (2.25)-(2.26). Condition (i) asserts effectively that essential spectrum of the linearized evolution operator lies in $\mathrm{Re}(\lambda)\leq 0$, condition (ii) that point spectrum lies in $\mathrm{Re}(\lambda)\leq 0$; see [31] for further discussion. Failure of weak Lopatinsky stability, also known as strong instability, implies exponential linear instability.

Corollary 2.19**.**

Assuming involution property (2.10), a Lax type shock solution (2.21) of (2.6) in direction $x_{1}$ with speed $s\not\in\sigma(\mathcal{M}_{1}|_{{\rm Range}\,\Gamma_{1}})$ satisfies the weak Lopatinsky condition with constraint (2.9) if and only if it satisfies the condition without the constraint.

Proof.

In the one-dimensional case $\xi=0$, condition (i) reduces to the requirement that the shock be Lax type. Since, by Proposition2.11, the type of the shock is the same with or without constraint, this verifies the result in dimension one. But, by Hersh’s Lemma (cf. [32], or [11, Chap. 4, Lemma 4.1]), satisfaction of this condition in multi-dimensions is equivalent for weakly hyperbolic systems to satisfaction in dimension one, whence by Proposition 2.2 the result follows also in multi-d.

Namely, one observes that existence of a pure imaginary eigenvalue $\mu=ik_{1}$ of $-\mathcal{A}_{1}^{-1}(\lambda+i\sum_{j\neq 1}\xi_{j}\mathcal{A}_{j})$ with associated eigenvector $\mathcal{R}$, implies $-\mathcal{A}_{1}^{-1}(\lambda+i\sum_{j\neq 1}\xi_{j}\mathcal{A}_{j})\mathcal{R}=ik_{1}\mathcal{R},$ and thus $(ik_{1}A_{1}+\sum_{j\neq 1}i\xi_{j}\mathcal{A}_{j})\mathcal{R}=-\lambda\mathcal{R}$, yielding $\lambda$ real by weak hyperbolicity, hence contradicting $\mathrm{Re}(\lambda)>0$. But this implies that no eigenvalues cross the imaginary axis as $\xi$ and $\lambda$ are varied within $\xi\in\mathbb{R}^{d-1}$, $\mathrm{Re}(\lambda)>0$, and so the number of positive/negative real part eigenvalues remains constant, in particular equal to that in the dimension one case $\xi=0$. Combining our conclusions in one- and multi-d, we find that condition (i) holds automatically, with or without constraint. But, by Corollary 2.15, condition (i) holds without constraint if and only if it holds with constraint. ∎

Corollary 2.20**.**

For stationary MHD shocks in the $x_{1}$ direction of either the $\beta$-model (2.1), $\beta\neq 0$, or of the standard (hybrid) model of [14, 45, 26], weak Lopatinsky stability holds with the divergence free constraint if and only if it holds without the constraint.

Proof.

The result for the $\beta$-model follows by Corollary 2.19. The result for the hybrid model follows by the same argument, noting that all relevant properties hold for the hybrid model as well. ∎

A convenient way to test for weak Lopatinsky stability is via the Lopatinsky determinant

[TABLE]

defined for all $\xi\in\mathbb{R}^{d-2}$, $\mathrm{Re}(\lambda)>0$, where $\mathscr{E}^{\pm}(\lambda,\xi)$ denote manifolds of (spatially) decaying solutions of (2.25) on $x_{1}\gtrless 0$; see [45, 60], or Section 3.1 and (3.9) for further discussion. Evidently, vanishing of the determinant $\Delta$ is equivalent to existence of a decaying normal mode, or solution of (2.25)–(2.26); hence, nonvanishing of the Lopatinsky determinant on $\mathrm{Re}(\lambda)>0$ is equivalent to weak Lopatinsky stability, a condition that is useful both numerically and analytically. In the present context, by Corollary 2.20, one may test weak Lopatinsky stability using either the hybrid model of [14, 45, 26] or, as we do here, the $\beta$-model, both without imposition of the divergence-free constraint. The advantage of the $\beta$ model is that it is of a form determined by a single set of equations, so treatable by existing, well-tested code: for example, the numerical stability package STABLAB [5], applicable to general models (2.3).

2.1.6. Uniform Lopatinsky stability

We next recall the uniform Lopatinsky condition of Majda [44], which requires, further, that for $\mathrm{Re}(\lambda)>0$ and $|\lambda,\xi|=1$, the Lopatinsky determinant is uniformly bounded, $|\Delta(\lambda,\xi)|\geq c>0$. Under various additional structural conditions on (2.6) (cf. [45]), the subspaces $\mathcal{E}^{\pm}$ of decaying modes of (2.25) on $x\gtrless 0$, hence also the Lopatinsky determinant $\Delta$, may be extended continuously to the closure $\mathrm{Re}(\lambda)\geq 0$, in which case uniform Lopatinsky stability is equivalent to nonvanishing of $\Delta$ on $\mathrm{Re}(\lambda)\geq 0$ except at the origin $(\lambda,\xi)=(0,0)$. Under additional structural conditions on (2.6), one may conclude further, by the symmetrizer method of Kreiss, that uniform Lopatinsky stability implies uniform resolvent bounds on nearby perturbations of the stationary shock, yielding also nonlinear stability, defined as well-posedness of the associated time-evolution problem [44, 45].

For MHD, such conditions were verified for the standard (hybrid) model in [45], yielding the following sufficient condition for nonlinear stability.

Proposition 2.21** ([45]).**

For shock waves of MHD, uniform Lopatinsky stability for the linearization of the hybrid model of [14, 45, 26], ignoring the divergence-free constraint, yields nonlinear stability of the the perturbation equations about a Lax-type shock, with or without the constraint.

Proposition 2.21 does not assert that uniform Lopatinsky stability with the constraint is equivalent to stability without the constraint. And, indeed, Corollary 2.15, though still true for $\mathrm{Re}(\lambda)=0$, is not much use, since the continuous extension to $\mathrm{Re}(\lambda)=0$ of decaying subspaces $\mathcal{E}^{\pm}$ are no longer necessarily decaying, but may contain neutral, oscillatory modes. Thus, vanishing of the Lopatinsky determinant implies only a solution of (2.26) belonging to special subspaces at plus and minus spatial infinity, and not existence of a decaying normal mode.

However, we can remedy this with a bit of further investigation. Following [45], define $(i\tau,\xi_{0})$, $\xi_{0}),\tau\in\mathbb{R}$, to be a nonglancing point for (2.25) if the stable eigenvalues $\mu^{\pm}(\lambda,\xi)$ of (2.25) and the unstable eigenvalues $\nu^{\pm}(\lambda,\xi)$ have distinct limits as $(\lambda,\xi)\to(i\tau,\xi_{0})$ from $\{\mathrm{Re}(\lambda)>0\}$; otherwise, call it a glancing point for (2.25). (This simplified definition may be seen to be equivalent to the notions of glancing defined in [45] for the cases considered there, including MHD.) Call $(i\tau,\xi_{0})$, $\xi_{0}),\tau\in\mathbb{R}$, a point of continuity for (2.25) if the associated stable subspaces $\mathcal{E}^{\pm}$ have limits as $(\lambda,\xi)\to(i\tau,\xi_{0})$ from $\{\mathrm{Re}(\lambda)>0\}$. At a point of continuity, call $\widehat{\mathcal{V}}$, piecewise smooth, an extended normal mode if it satisfies (2.25)-(2.26), and approaches $\lim_{(\lambda,\xi)\to(i\tau+0^{+},\xi)}\mathcal{E}^{\pm}$ as $x_{1}\to\pm\infty$.

Proposition 2.22**.**

Let $(i\tau,\xi_{0})$, $\tau,\xi_{0}\in\mathbb{R}$, be a point of continuity for (2.25), and a nonglancing point for

[TABLE]

Then, assuming (2.10), (2.24), for an extended normal mode $\widehat{\mathcal{V}}$, we have $\widehat{\Gamma}\widehat{\mathcal{V}}\equiv 0$. That is, extended normal modes automatically satisfy the Laplace-Fourier transformed version of constraint (2.9).

Proof.

Let $e_{\pm}$ and $f_{\pm}$ denote the subspaces of decaying and growing solutions of (2.34) as $x_{1}\to\pm\infty$ for $\Re\lambda>0$, and their continuous extensions to $\Re\lambda\geq 0$ (well-defined, by the nonglancing assumption for (2.34)).

For $\mathrm{Re}(\lambda)>0$, subspaces $\mathcal{E}^{\pm}$ correspond to decaying modes $\widehat{\mathcal{V}}^{\pm}$ at $x_{1}\to\pm\infty$ of

[TABLE]

whence, by Lemma 2.14, $\hat{\phi}^{\pm}:=\hat{\Gamma}\widehat{\mathcal{V}}^{\pm}=(\Gamma_{1}G^{\pm}+\sum_{j\neq 2}i\xi_{j}\Gamma_{j})\widehat{\mathcal{V}}^{\pm}$ is a decaying weak solution of (2.34), hence lies in the decaying subspaces $e^{\pm}$ of $-\mathcal{M}_{1}(\lambda j+\sum_{j\neq 2}i\xi_{j}\mathcal{M}_{j})$ as $x_{1}\to\pm\infty$, with jump condition $[\mathcal{M_{1}}\hat{\phi}]=\mathcal{M}_{1}\hat{\phi}=0$, or, by invertibility of $\mathcal{M}_{1}$ (recall: a consequence of (2.24), $[\hat{\phi}]=0$.

By the assumed continuity of $\mathcal{E}^{\pm}$, we have at $(i\tau,\xi_{0})$ that normal modes $\widehat{\mathcal{V}}$ may be realized as limits of decaying normal modes on $\mathrm{Re}(\lambda)>0$ in the limit as $(\lambda,\xi)\to(i\tau,\xi_{0})$, hence $\hat{\phi}=\hat{\Gamma}\widehat{\mathcal{V}}$ (piecewise continuous) is the limit of decaying solutions of (2.34), i.e., vectors lying in $e^{\pm}$ for $x_{1}\gtrless 0$. By the assumed nonglancing of $(i\tau,\xi_{0})$ with respect to (2.34)- specifically, the implied continuity of $e^{\pm}$ as $(\lambda,\xi)\to(i\tau,\xi_{0})$ from $\mathrm{Re}(\lambda)>0$- $\hat{\phi}$ must thus lie in $e^{\pm}$ for all $x_{1}$, with $[\hat{\phi}]=0$ forcing $\hat{\phi}(0^{+})=\hat{\phi}(0^{-})$ at $x_{1}$=0. Noting, since $\mathcal{M}_{j}$ are constant, that $e^{\pm}=f^{\mp}$, and applying the nonglancing assumption as second time- specifically, the fact that $e^{\pm}$ and $f^{\pm}$ converge to transverse limiting subspaces, or $e^{+}\cap f^{+}=\emptyset$ and $e^{-}\cap f^{-}=\emptyset$, we find that $\hat{\phi}(0)$ must lie in $(e\cap f)^{\pm}=\{0\}$, and thus $\hat{\phi}=\hat{\Gamma}\widehat{\mathcal{V}}\equiv 0$ as claimed. ∎

Before discussing continuity of subspaces, we establish a preliminary result of interest in its own right, comprising the Laplace–Fourier transform analog of (2.13).

Lemma 2.23**.**

Assuming (2.10), (2.24), let $G_{\pm}(\lambda,\xi)$ as in (2.35) be defined as $-\mathcal{A}_{1}^{-1}(\lambda+\sum_{j\neq 2}i\xi_{j}\mathcal{A}_{j})_{\pm}$, and $H_{\pm}(\lambda,\xi)$ and $L_{\pm}$ as $-\mathcal{M}_{1}^{-1}(\lambda+\sum_{j\neq 2}i\xi_{j}\mathcal{M}_{j})_{\pm}$ and $(\Gamma_{1}G+\sum_{j\neq 2}i\xi_{j}\Gamma_{j})_{\pm}$. Then,

[TABLE]

Proof.

Using (2.28) to express

[TABLE]

we may represent $L$ alternatively (dropping $\pm$ for ease of writing) as

[TABLE]

or, expanding and recombining using (2.37),

[TABLE]

Thus,

[TABLE]

which is equal to $H(\Gamma_{1}G+\sum_{j=2}^{n}i\xi_{j}\Gamma_{j})$ if and only if $M_{1}^{-1}\sum_{j=2}^{n}\Gamma_{j}A_{1}G=H\sum_{j=2}^{n}i\xi_{j}\Gamma_{j}),$ or

[TABLE]

as follows by inspection using (2.37) and comparing coefficients of $\lambda$, $\xi$, and $\xi^{2}$ terms. ∎

Corollary 2.24**.**

Assuming (2.10), (2.24), $\ker L_{\pm}$ is a right invariant subpace and ${\rm Range}L_{\pm}^{*}$ a left invariant subspace of $G_{\pm}$. Assuming $L_{\pm}(\lambda,\xi)$ is full rank $\ell$ for all $\lambda,\xi$, then we have the decomposition $\mathbb{C}=\ker L_{\pm}\oplus{\rm Range}L_{\pm}^{*}$, with $\ker L$ an $n-\ell$ dimensional analytically varying invariant subspace of $G_{\pm}$, consisting precisely of those normal modes $\hat{V}$ satisfying the constraint $\hat{\Gamma}\hat{V}=0$, while ${\rm Range}L^{T}$ gives a complementary antiholomorphic space.

For both the hybrid MHD system and the $\beta$-model, it is clear by inspection that $L_{\pm}$ is full rank one, hence the Corollary applies. Corollary 2.24 has the interesting interpretation as defining a “reduced” interior system in the Laplace-Fourier transform setting analogous to that of the Fourier settings, pertaining precisely to those modes satisfying the constraint $\Gamma$. Moreover, coordinatizing via the splitting $\ker L\oplus{\rm Range}L^{*}$, we may convert $G$ by a $C^{\infty}$ transformation ($G$, $L$ holomorphic in $\lambda$, $\xi$, and $L^{*}$ anti-holomorphic ) to the canonical upper block triangular form

[TABLE]

where the lower right-hand block is obtained from (2.36), via the computation $(LL^{*})^{-1}LGL^{*}=(LL^{*})^{-1}HLL^{*}=H.$

Defining nonglancing and points of continuity for the reduced systems $G|_{\ker L}$ and $H$ in the obvious way, we may at a point of continuity for $G_{\ker L}$, by a further continuous (not necessarily $C^{\infty}$) transformation, reduce to form

[TABLE]

where $g_{+}$ is the reduction of $G_{\ker L}$ to its stable subspace, and, at a nonglancing point of $H$, still further to

[TABLE]

where $h_{\pm}$ are the reductions of $H$ to its stable and unstable subspaces.

We connect continuity of the full and reduced sytems in two important cases.

Proposition 2.25**.**

Assume (2.10), (2.24), and $L_{\pm}$ is full rank for all $\lambda,\xi$.

1. (i)

If, local to $(i\tau,\xi_{0})$, $\tau,\xi_{0}\in\mathbb{R}$, there is a continuous right invariant space ${\rm Span}\{R_{1},\dots,R_{\ell}\}$ of $G$ such that $LR$ is invertible for $R$ defined as $R(\lambda,\xi)=(R_{1},\dots,R_{\ell})$, then $(i\tau,\xi_{0})$ is a point of continuity for (2.25) if and only if it is a point of continuity of both the associated reduced system and the complementary system $\partial_{x}z=Hz$. 2. (ii)

A point $(i\tau,\xi_{0})$, $\tau,\xi_{0}\in\mathbb{R}$ that is a nonglancing point for (2.34) is a point of continuity for (2.25) if it is a point of continuity of the associated reduced system $G|_{\ker L}$.

Proof.

(i) In this case, coordinatizing via the splitting $\ker L\oplus{\rm Range}(R)$ reduces $G$ to block-diagonal form $G={\rm blockdiag}\{G|_{\ker L},H\}$, whence the result follows trivially. (ii) Using the canonical coordinatization (2.41), we find, for $\mathrm{Re}(\lambda)>0$ and local to $(i\tau,\xi_{0})$, that $E_{+}=\begin{pmatrix}{\rm Id}&0&0&0\end{pmatrix}\oplus r_{+}$, where $r$ is of form $r_{+}(\lambda,\xi)=\begin{pmatrix}*&0&{\rm Id}&0\end{pmatrix}$; equivalently, that there exists $r_{+}$ of this form with range a right-invariant subspace of $G$ modulo $e_{+}$, where $e_{+}$ denotes the continuous extension of the stable subspace of $H$. (Here, we are using in an important way the fact that $\sigma(g_{-})$ and $\sigma(h_{+})$ are separated for $\mathrm{Re}(\lambda)>0$, in eliminating the second row of $r_{\pm}$.) But this evidently is equivalent to

[TABLE]

yielding continuity as a trivial consequence. ∎

Corollary 2.26**.**

For stationary MHD shocks in the $x_{1}$ direction of either the $\beta$-model (2.1), $\beta\neq 0$, or of the standard (hybrid) model of [14, 45, 26], all $(i\tau,\xi_{0})$ are points of continuity for the associated interior equation (2.25), for both reduced and full models, so long as

[TABLE]

Proof.

(hybrid model) It was shown in [45] that, assuming (2.42), all points $(i\tau,\xi_{0})$ are points of continuity for the hybrid model without constraint, a consequence of [45, Lemma 7.2] specialized to the 2D case (that is, ignoring middle, Alfven modes $\lambda_{2}^{\pm}$), together with [45, Theorem 5.2]. Moreover, the hyperbolic characteristic $\lambda_{0}(k)=\bar{u}_{1}k_{1}+\bar{u}_{2}k_{2}$ carrying nonzero divergence of $h$, with associated eigenvector $R_{0}(k)=(0,0,0,k_{1},k_{2})$, may be directly inverted to obtain an analytic eigenvalue $\mu_{0}(\lambda,\xi)=\bar{u}_{1}^{-1}(\lambda+\bar{u}_{2}i\xi)$ of $G=-\mathcal{A}_{1}^{-1}(\lambda+i\xi_{0}\mathcal{A}_{2})$ with associated eigenvector $S_{0}=(0,0,0,\mu/i,\xi)$; see [45, (A.4)]. For $\mu_{0}(i\tau,\xi_{0})=ik_{1}$ pure imaginary, and $(\lambda,\xi)\to(i\tau,\xi_{0})$, this is not annihilated by the constraint, since

[TABLE]

as $\mu_{0}\to ik_{1}$. But, this implies also $LS_{0}\neq 0$, since $LS_{0}=(\Gamma_{1}\mu+i\xi\Gamma_{2})S_{0}$. Such points are therefore points of continuity for the reduced system by Proposition 2.25(i). If $\mu_{0}$ is not pure imaginary, on the other hand, then for $(\lambda,\xi)$ near $(i\tau,\xi_{0})$, $\lambda_{0}$ is strictly stable, and so, along with other strictly stable modes, continues analytically and transversally to the space $\ker(L)$, so that continuity for the full system reduces also in this case to continuity of the stable subspace for the reduced problem. Thus, we conclude from continuity for the full system continuity for the reduced system as well.

($\beta$-model) For both models, the complementary system associated with $H$, being scalar, is necessarily nonglancing. Since we have already shown continuity for the reduced system via our treatment of the hybrid model, we thus obtain continuity for the full $\beta$-model by Proposition 2.25(ii). ∎

We define reduced weak Lopatinsky stability as nonexistence of solutions of (2.25)-(2.26) lying in $\mathcal{E}^{\pm}$ for $x_{1}\gtrless\infty$. By Corollary 2.19, this is equivalent to weak Lopatinski stability for the full system, without constraint.

When the decaying subspaces $\mathcal{E}^{\pm}_{0}$ of the reduced system extend continuously for all $\lambda=i\tau$, we may define a notion of reduced uniform Lopatinsky stability as nonexistence of solutions of (2.25)-(2.26) lying in the limiting spaces $\mathcal{E}^{\pm}_{0}$ for $x_{1}\gtrless 0$. When the subspaces $\mathcal{E}^{\pm}$ of the full system extend continuously as well, we have by Proposition 2.22 that full and reduced uniform Lopatinski stability are equivalent.

Combining results, we have as an immediate consequence of Corollary 2.26 the following definitive conclusion, elucidating and significantly extending the fundamental observation of Blokhin-Trakhinin [14, 26] that weak stability for constrained MHD is equivalent to weak stability without constraint. In particular, this shows for the first time that the nonlinear stability condition of Proposition 2.21 is sharp.

Proposition 2.27**.**

For stationary MHD shocks in the $x_{1}$ direction of either the $\beta$-model (2.1), $\beta\neq 0$, or of the standard (hybrid) model of [14, 45, 26], uniform Lopatinsky stability holds with the divergence free constraint if and only if it holds without the constraint, if and only if it holds for the reduced model defined by $G|_{\ker L}$.

2.2. The viscous case

We conclude this section by a brief treatment of the analogous but simpler viscous case, generalizing and expanding on the approach used in [47, pp. 2 and 62-64] for the equations of viscoelasticity (there introduced somewhat implicitly and without particular emphasis in the course of other computations).

System (2.1), in the viscous case $\mu,\eta,\nu>0$, may be expressed in the general form

[TABLE]

augmenting the inviscid system (2.6), with linearization about a constant solution $\mathcal{V}\equiv\mathcal{V}_{0}$ given by

[TABLE]

where $\mathcal{A}_{j}=D\mathcal{F}_{j}(\mathcal{V}_{0})$ and $\mathcal{B}_{j,k}=\mathcal{B}_{j,k}(\mathcal{V}_{0})$. Likewise, Dafermos’ involution condition (2.10) is replaced by a viscous counterpart

[TABLE]

where $\mathcal{M}$ is now second-order constant-coefficient:

[TABLE]

and (strictly) parabolic:

[TABLE]

As in the inviscid case, (2.44) induces the linearized version

[TABLE]

2.2.1. Dissipativity

At a constant state $\mathcal{V}_{0}$, the viscous condition analogous to weak hyperbolicity in the inviscid case is weak dissipativity:

[TABLE]

a condition on the dispersion relation for $L$ related to the genuine coupling condition of Kawashima [41, 56, 57, 30]. From (2.46)–(2.47), we find, applying $\Gamma(k)$ to the eigenvalue equation $(L(k)-\lambda)\mathcal{V}=0$ that $(-\mathcal{M}(k)-\lambda)\Gamma(k)\mathcal{V}=0$ for $\mathcal{M}(k):=\sum_{j=1}^{d}\mathcal{M}^{1}_{j}k_{j}+\sum_{i,j=1}^{d}\mathcal{M}^{2}_{i,j}k_{i}k_{j}$, whence $\Gamma(k)\mathcal{V}=0$ or else $\lambda\in\sigma(-\mathcal{M}(k))$ verifying the weak dissipativity condition by (2.46). Thus, similarly as in Proposition 2.2 of the inviscid case, weak dissipativity holds with constraint $\Gamma(k)\mathcal{V}=0$ if and only if it holds without the constraint.

It was verified in [45] that (2.48) holds for the hybrid MHD model with and without constraint, whence it holds for the $\beta$-model as well.

2.2.2. Persistence of constraints

From (2.45)–(2.46), we obtain immediately persistence of the constraint $\Gamma\mathcal{V}=0$ for smooth solutions of (2.43), similarly as in Corollary 2.7 for the inviscid case. For, applying $\Gamma$ to (2.43) gives

[TABLE]

with initial data $\mathcal{V}\equiv 0$, whence we obtain $\Gamma\mathcal{V}=0$ for all $t\geq 0$ by well-posedness of (2.49). Recall that the standard viscous existence/stability theory concerns smooth ($H^{s}$) solutions [56].

Remark 2.28*.*

In the case of the MHD $\beta$-model, $Mz=-\beta\partial_{x_{1}}z+\nu\Delta z$, and $M\Gamma\mathcal{V}=-\beta\partial_{x_{1}}(\rm{div}h)+\nu\Delta({\rm div}h)$, so that solutions of (2.43) satisfy

[TABLE]

For the hybrid model of [14, 45, 26], solutions satisfy $0=\partial_{t}({\rm div}h)+\nabla\cdot(u(\rm{div}h))+\nu\Delta({\rm div}h)$, giving persistence by a similar argument to the above [45].

2.2.3. Compatibility of shock profiles

Likewise, similarly as in Proposition 2.11, we find for smooth viscous shock profiles

[TABLE]

$|k|=1$, that shock profiles of the unconstrained system automatically satisfy the constraint. For, applying $\Gamma$ to the traveling-wave ODE $-sk\cdot\nabla\mathcal{V}+\mathcal{F}(\mathcal{V})=0$, rearranging, and integrating from $-\infty$ to $z$ gives

[TABLE]

where $\gamma=\Gamma\mathcal{V}=\Gamma(k)\overline{\mathcal{V}}^{\prime}$ and $\mathcal{N}(k)=(\mathcal{M}^{2}(k))^{-1}(\mathcal{M}^{1}(k)-s)$ is constant-coefficient. (Here, $\mathcal{M}^{2}(k)(k)$ is invertible by parabolicity.) Noting that $\gamma(z)\to 0$ as $z\to\pm\infty$ by convergence of the shock profile to endstates $\mathcal{V}_{\pm}$, we find that $\gamma\equiv 0$, since there are no nontrivial solutions of a constant-coefficient ODE decaying at both $\pm\infty$.

2.2.4. Compatibility of normal modes

Restricting without loss of generality to the case of a zero-speed planar shock in direction $x_{1}$, $\mathcal{V}=\overline{\mathcal{V}}(x_{1})$, linearizing, and taking the Fourier transform in transverse directions $x_{2},\dots,x_{d}$, we obtain [60, 56, 57] the generalized eigenvalue, or normal modes equation

[TABLE]

with $(\xi_{2},\dots,\xi_{d})\in\mathbb{R}^{d-1}$ corresponding to Fourier frequencies, and existence of decaying solutions $\widehat{\mathcal{V}}$ for $\mathrm{Re}(\lambda)>0$ corresponding to exponential linear instability.

Denoting by $\widehat{\mathcal{M}}(\xi)$ the Fourier transform in directions $x_{2},\dots,x_{d}$ of $\mathcal{M}$, we have from (2.44), evidently, $\hat{\Gamma}\widehat{\mathcal{L}}=-\widehat{\mathcal{M}}\hat{\Gamma}$, whence normal modes, similarly as in the inviscid case, satisfy also the constant-coefficient ODE

[TABLE]

By the absence of nontrivial solutions of constant-coefficient ODE that decay at both $\pm\infty$, we find immediately as in Corollary 2.15 for the inviscid case, that decaying normal modes automatically satisfy the constraint, i.e., $\hat{\Gamma}\widehat{\mathcal{V}}=0$.

2.2.5. Weak and strong Evans-Lopatinski stability

In general, (2.48) is assumed to hold at the endstates $\mathcal{V}_{\pm}$ of a shock, whence, by a standard Hersch-type lemma [60, 56, 57], one obtains “consistent splitting” of the eigenvalue equations (2.51) on $\{(\xi,\lambda):\,\mathrm{Re}(\lambda)\geq 0\}\setminus\{0,0\}$, that is, the property that the dimensions of the subspaces of decaying solutions at $+\infty$ and $-\infty$ of (2.51) sums to the total dimension $N$ of the solution space, as do the dimensions of the subspaces of exponentially growing solutions at $+\infty$ and $-\infty$. It follows that one may define an Evans function $\mathcal{D}$ consisting of the Wronskian of $N$ solutions comprising bases of the decaying solutions at $\pm\infty$, with vanishing of $\mathcal{D}(\xi,\lambda)$ corresponding to existence of a decaying normal mode for $\xi,\lambda$. Following [60, 56, 57, 30], we define strong Evans-Lopatinsky stability [60, 56, 57, 30] as the absence of decaying normal modes- or, equivalently, zeros of the Evans function $\mathcal{D}$- on $\{(\xi,\lambda):\,\mathrm{Re}(\lambda)\geq 0\}\setminus\{0\}$, and weak Evans-Lopatinsky stability as the absence of zeros of the Evans function on $\{(\xi,\lambda):\,\mathrm{Re}(\lambda)>0\}$. Failure of weak Evans-Lopatinski stability is also called strong Evans-Lopatinski instability, implying exponential instability of the background shock.

From the definition, and the observations in Section 2.2.4, we obtain immediately the viscous analogs of Corollaries 2.19 and 2.20, showing that both weak and strong Evans-Lopatinsky stability hold for the $\beta$- or hybrid model of MHD with constraint if and only if they hold without constraint.

2.2.6. Uniform Evans-Lopatinsky stability and nonlinear stability

The uniform Evans-Lopatinsky condition, analogous to the uniform Lopatinsky condition of the inviscid case, is a uniform lower bound $|\mathcal{D}(\xi,\lambda)|\geq\theta/|(\xi,\lambda)|$, $\theta>0$, on $\{(\xi,\lambda):\,\mathrm{Re}(\lambda)\geq 0\}\setminus\{0\}$. In nice cases, in particular, for fast Lax shocks in MHD, this is equivalent to strong Evans-Lopatinsky stability plus transversality of the shock profile as a solution of the traveling-wave ODE, plus satisfaction of the uniform Lopatinsky condition for the associated inviscid shock; moreover, uniform Evans-Lopatinsky stability is equivalent to nonlinear stability of the underlying planar shock both with respect to large time and small viscosity asymptotics [56, 57, 30].

In these favorable cases, uniform Evans-Lopatinsky stability with constraint is evidently (from our previously remarked oobservations, above) equivalent to uniform Evans-Lopatinsky stability without constraint. On the other hand, the relation between uniform Evans-Lopatinsky stability and uniform stability of the associated inviscid shock has not been shown for slow MHD shocks; whether or not this holds has been cited in [30] as an important open problem in the nonlinear stability theory for MHD shocks. Thus, in the present setting, issues pertaining to uniform Evans-Lopatinski, nonlinear stability, and their relations, remain unclear. As our main object is instability, far from the origin $(\xi,\lambda)=(0,0)$, this is no obstruction in the present study. However, we mention it as an important aspect of the theory for further development.

3. Inviscid stability analysis: the Lopatinsky determinant

The main goal in this section is to introduce the Lopatinsky determinant, which will be the main tool used in the study of inviscid stability. Initially we derive some properties of the shock type. Recall that shocks are categorized as Lax, undercompressive and compressive depending on the number of characteristics entering the shock. Lax shocks are further categorized by their characteristic field: the unique family entering on both sides. To begin with, we obtain a useful parametrization of the shocks we study.

Recall that the Rankine-Hugoniot conditions (2.22) specialized to a stationary shock $\mathcal{W}(x,t)=\overline{\mathcal{W}}(x_{1})=\mathcal{W}_{\pm}$ in the $x_{1}$ direction of the MHD equations (1.1) are

[TABLE]

where $[f_{1}(\overline{\mathcal{W}})]:=f_{1}(\overline{\mathcal{W}}^{+})-f_{1}(\overline{\mathcal{W}}^{-})$.

Lemma 3.1** (Parametrization of MHD planar shocks with zero speed).**

Let

[TABLE]

$x_{1}\gtrless 0$, be a planar shock solution satisfying the Rankine-Hugoniot conditions in (3.6). Assume that the shock is parallel, i.e., $u_{2}^{+}=0$ and $h_{2}^{+}=0$.

1. (i)

One can parametrize the slow shocks connecting to the state $(\overline{\rho}^{+},u_{1}^{+},0,h_{1}^{+},0)$ to the right (i.e., $x_{1}>0$) using variables $R$ and $M$, defined as

[TABLE]

where $M^{2}=\frac{(u_{1}^{+})^{2}}{p_{\rho}(\rho^{+})}$ ($M$* is also known as the downstream Mach number), and $u_{2}^{-}=0$, $h_{2}^{-}=0$, $h_{1}^{+}=h_{1}^{-}$;* 2. (ii)

For large magnetic field, slow shocks, i.e., 2-shocks, are characterized by $R>1$ and $M<1$; 3. (iii)

There exist two numbers

[TABLE]

according to which three scenarios are possible: fast Lax shocks (i.e., extreme, or gas-dynamical type) for $0\leq|h_{1}|\leq H_{*}$, intermediate shocks for $H_{*}\leq|h_{1}|\leq H^{*}$ and slow shocks for $H^{*}\leq|h_{1}|$ (in particular $|h_{1}|\to\infty$).

We remark that in the context of 3D MHD, the 2-D slow shocks seen here become $3$-shocks

Proof.

Considering the jump conditions in (3.6) we readily observe that the ratios $\frac{\rho^{+}}{\rho^{-}}$ and $\frac{u_{1}^{-}}{u_{1}^{+}}$ are in fact equal, so that one can define the variable $R$ as in (3.7). Both properties $u_{2}^{-}=0$, $h_{2}^{-}=0$ follow from solving the third and fifth rows of the jump conditions in (3.6); the fourth equation in (3.6) implies that $h_{1}^{+}=h_{1}^{-}$. Setting $M:=\frac{u_{1}^{+}}{\sqrt{p_{\rho}(\rho^{+})}}$, we obtain the rightmost condition in (3.7) from the second relation in the jump condition (3.6). This establishes (i).

With regards to (ii), classical compressibility conditions in gas dynamics, namely $p_{\rho}(\cdot)>0$, gives that $R>1$; by inspection of (3.7), one can see that this condition implies that $0\leq M<1$, for $\gamma\geq 1$.

Last, we prove (iii). The type of the shock is determined by the number of incoming characteristics, or eigenvalues of $A_{0}^{-1}A_{1}$, in the $x_{1}$ direction on either side of the shock, where $A_{j}$ are as in (2.5). Computing $A_{0}^{-1}A_{1}$ and its spectrum on both sides of the shock, we have:

[TABLE]

In the large magnetic field scenario $h_{1}\to\infty$, $\beta>0$ and a zero speed shock we count 2 negative/3 positive (resp. 1 negative/4 positive) eigenvalues for $x_{1}>0$ (resp., $x_{1}<0$) whenever $u_{1}^{\pm}\pm\sqrt{p_{\rho}^{\pm}}\lessgtr 0$. Using the definitions in (3.7) we conclude that this is equivalent to $M<1$ when we use the constraints for $x_{1}>0$; similarly, we derive that $R>1$ using the constraints for $x_{1}<0$. Therefore, the discussion above says that we just need to analyze the signs of $1-\frac{h_{1}}{u_{1}^{+}\sqrt{\rho^{+}}}$ and $1-\frac{h_{1}}{u_{1}^{-}\sqrt{\rho^{-}}}=1-\frac{h_{1}}{u_{1}^{+}\sqrt{R\rho^{+}}}$. As $R>1$ we have

[TABLE]

Therefore, recalling $R=\frac{u_{1}^{-}}{u_{1}^{+}}$, we obtain the following table corresponding with (1.3):

∎

Notice that whenever $\beta>0$ the $\beta$-model preserves slow shocks, but this property can also be verified whenever a positive multiple of $\mathrm{div}(h)$ is added to any upstream/downstream side of the propagating shock. In particular, one may add different multiples to upstream and downstream sides, as for example in [45, 26], allowing one to symmetrize the equations while preserving 2-shock structure.444 Nonconstant $\beta$ is inconvenient however for the viscous case, destroying conservative structure.— Last, we remark that the case $\beta=0$ is degenerate, for the matrix $A_{1}$ is not invertible in this case; the invertibility of $A_{1}$ is necessary in the construction of the Lopatinsky determinant, which we discuss next.

3.1. Lopatinsky determinant: construction and asymptotic analysis

In this section we study the Lopatinsky determinant associated with inviscid parallel shocks and the onset of instability. Initially, we study the behavior of the system (2.5) by taking its Laplace-Fourier transform (in $t$ and $x_{2}$, respectively). We obtain

[TABLE]

where $v$ is the Laplace-Fourier transform of $u$ and

[TABLE]

where the variables $(u_{1},\rho,p_{\rho})$ should be read $(u_{1}^{+},\rho^{+},p_{\rho^{+}})$ (resp., $(u_{1}^{-},\rho^{-},p_{\rho^{-}})$) whenever $x_{1}>0$ (resp. $x_{1}<0$). Equation (3.8) consists of two systems of ODEs in the interior: one in $x_{1}>0$, one in $x_{1}<0$, coupled through Rankine-Hugoniot conditions at $x_{1}=0$ (Appendix C; see also [60] for further discussion on the technique). The Lopatinsky determinant is defined as

[TABLE]

The parameter $\lambda$ is a spectral parameter indicating solutions to the system (2.5) with growth $\sim e^{\lambda t}$ (thus, $\mathrm{Re}(\lambda)>0$ corresponding to instability); $\mathscr{E}^{\pm}(\lambda,\xi)$ denote manifolds of (spatial) decaying solutions in $x_{1}\gtrless 0$. In the rest of the paper we omit the dependence of these spaces on $\lambda$ and $\xi$, simply writing $\mathscr{E}^{\pm}$.

3.1.1. Large-$h_{1}$ asymptotics

Following [26, 15] we now study the large magnetic field $\overline{h}_{1}$ asymptotics. It is convenient to define the quantity $\displaystyle{\varepsilon:=\frac{1}{\overline{h}_{1}}}$, which parametrizes the underlying viscous profile $\overline{\mathcal{V}}=\overline{\mathcal{V}}^{(\varepsilon)}$. Our analysis consists of Taylor expanding the roots of the Lopatinsky determinant defined in (3.9) considered as a function in $\varepsilon$ with $\xi$ held fix at $1$, i.e., $\varepsilon\mapsto\lambda(\varepsilon)$. It is shown that

[TABLE]

where $\lambda_{2}>0$. One can conclude that an unstable regime occurs in the large magnetic field scenario, as verified in [26].

The study of the spaces $\mathscr{E}^{\pm}$ is equivalent to analyzing the eigenvalues of $-(A_{1}^{\pm})^{-1}(\lambda A_{0}^{\pm}+i\xi A_{2}^{\pm})$ and their associated eigenspaces. The approach we adopt relies on careful analytical estimates allied with the use of symbolic computations (carried out in SAGE [49] [8]). The main idea is the following: assume that the spectral parameter $\lambda$ can be expanded as

[TABLE]

Thus each element $\mu\in\sigma(-(A_{1}^{\pm})^{-1}\left(\lambda A_{0}^{\pm}+i\xi A_{2}^{\pm}\right))$ can be expanded as

[TABLE]

Afterwards we find a set of eigenvectors of $-(A_{1}^{\pm})^{-1}\left(\lambda A_{0}^{\pm}+i\xi A_{2}^{\pm}\right)$ spanning the spaces $\mathscr{E}^{\pm}$; namely, a mapping

[TABLE]

For the sake of convenience we drop the indexes “$\pm$” for now, since these formulas work for both cases if one use the notation in (3.7). If we also expand the eigenvector $X(\mu)$ in (3.13) in $\varepsilon$ terms we have

[TABLE]

Observation 3.2**.**

With regards to the symbolic computations, a few remarks are in order:

1. (i)

In order to find the expansion in (3.12) we find the characteristic polynomial $p(\cdot)$ of the matrix $-(A_{1}^{\pm})^{-1}\left(\lambda A_{0}^{\pm}+i\xi A_{2}^{\pm}\right)$ when (3.11) holds; we conclude by matching coefficients. The characteristic polynomial $p(\mu)$ can be shown to expand as $p(\mu)=p(\mu_{0},\mu_{1},\mu_{2},\ldots)=p_{0}(\mu_{0})+\varepsilon p_{1}(\mu_{0},\mu_{1})+\varepsilon^{2}p_{2}(\mu_{0},\mu_{1},\mu_{2})+\ldots$, with $p_{j}$ computed explicitly using SAGE; 2. (ii)

Higher order terms in the expansion (3.12) in terms of $\varepsilon$ can be easily obtained using SAGE, because once we have $\mu_{0}\,\ldots\,\mu_{i}$ the problem of finding $\mu_{i+1}$ is linear; these terms won’t be written explicitly here though; 3. (iii)

It is not hard to show that $-\frac{\lambda}{\beta}$ is an eigenvalue of $-(A_{1}^{\pm})^{-1}\left(\lambda A_{0}^{\pm}+i\xi A_{2}^{\pm}\right)$. This also holds true in the nonparallel case, as we will show later; 4. (iv)

Notice that, upon scaling, we can take $\xi=1$ (redefine $\lambda\to\lambda\xi$); 5. (v)

Let $X_{i}$ denote the $i^{th}$ column of a square matrix $A=[X_{1},X_{2},\ldots X_{n}]$ and $X_{i}(\varepsilon)=X_{i}^{(0)}+\varepsilon X_{i}^{(1)}+\varepsilon^{2}(X_{i}^{(2)})+\mathcal{O}(\varepsilon^{3})$. Then

[TABLE]

where $A_{k}=\underset{a_{1}+\ldots+a_{n}=k}{\sum}det(X_{1}^{(a_{1})},\ldots,X_{n}^{(a_{n})})$.

Taking into account these observations, we plug (3.11), (3.12) and (3.13) into (3.9), to obtain

[TABLE]

In what follows we shall exploit the multi-linearity of the determinant function in order to find the terms in this expansion; each term $\Delta_{i}$ is a function of $\left(\lambda_{i}\right)_{i\in\mathbb{N}},$ but this dependence will be most of the time omitted and made explicit as we carry out our computations. We make use of (3.14) to verify the condition $\Delta^{(\varepsilon)}(\lambda,1)\equiv 0$: we look for $(\lambda_{i})_{i\in\mathbb{N}}$ that gives $\Delta_{i}(\lambda)=0$ for all $i\in\mathbb{N}.$ An important step in the analysis consists of an explicit representation of the manifolds $\mathscr{E}^{\pm}$ referred to in (3.9), an investigation that we reformulate as the study of the eigenvalues of the matrices in (3.8). Indeed, one can observe by inspection that $-\frac{\lambda}{\beta}$ is an eigenvalue of $-(A_{1}^{\pm})^{-1}\left(\lambda A_{0}^{\pm}+iA_{2}^{\pm}\right)$, while an expansion of the eigenvalues of the form (3.12) readily shows that their zeroth order terms $\mu_{0}$ for $x_{1}\gtrless 0$ are

[TABLE]

where all the variables $(u_{1},\rho,p_{\rho})$ should be read $(u_{1}^{+},\rho^{+},p_{\rho^{+}})$ (resp., $(u_{1}^{-},\rho^{-},p_{\rho^{-}})$) whenever $x_{1}>0$ (resp. $x_{1}<0$). Therefore, since we are looking for “decaying” manifolds we must have, for $x_{1}>0$, eigenspaces associated to the following eigenvalues:

[TABLE]

Analogously, “decaying” manifolds in $x_{1}<0$ must be eigenspaces associated to the eigenvalue $\displaystyle{1+\mathcal{O}(\varepsilon^{2})}$. Notice that the number of eigenvalues in each interior $x_{1}\gtrless 0$ is consistent with the analysis derived from Hersh’s Lemma, hence it suffices to analyze the number of positive and negative eigenvalues of $-A_{1}^{-1}A_{0}$ .

3.1.2. Asymptotic instability

Choosing an appropriate parametrization of the decaying manifolds $\mathscr{E}^{\pm}$ one can show that $\Delta_{0}=\Delta_{1}=0$ in equation (3.14). A careful computation shows then that $\Delta_{2}=\Delta_{2}(\lambda_{0},\mu_{0})=det(X_{1}^{(1)},X_{2}^{(0)},X_{3}^{(0)},X_{4}^{(1)},X_{5}^{(0)}).$ We have

[TABLE]

where $\displaystyle{E_{*}=\frac{i\,\lambda_{0}(u_{1}^{+})^{2}}{{\left(u_{1}^{+}+\sqrt{p_{\rho^{+}}}\right)}p_{\rho^{+}}}+\frac{i\,\lambda_{0}^{3}u_{1}^{+}}{{\left(u_{1}^{+}+\sqrt{p_{\rho^{+}}}\right)}^{2}p_{\rho^{+}}}-\frac{i\,\lambda_{0}^{3}(u_{1}^{+})^{2}}{{\left(u_{1}^{+}+\sqrt{p_{\rho^{+}}}\right)}^{3}p_{\rho^{+}}}-\frac{i\,\lambda_{0}u_{1}^{+}}{p_{\rho^{+}}}.}$ The block structure of this matrix allow us to see right away that $\Delta_{2}=\mathcal{O}(\lambda_{0}^{3})$ and $\Delta_{2}=o(\lambda_{0})$. There exists a $\lambda_{0}=\mathcal{O}(1)$ satisfying $\Delta_{2}=0$ at which, however, the parametrization of the manifolds $\mathscr{E}^{\pm}$ through their eigenvectors is lost because two of those eigenvectors coincide. These points are called glancing modes; this issue was also pointed out in [26, §3.4]. Resorting to generalized eigenvectors one can show that the Lopatinsky determinant does not vanish for this value of $\lambda_{0}$, hence the only solution to $\Delta_{2}=0$ is $\lambda_{0}=0$; see also [8].

In the search for roots of $\Delta_{3}=0$ the now look for $\lambda_{1}$. Thanks to the multilinearity of the determinant function, we readily obtain that $\Delta_{3}=0$. We go to the next term in (3.14): thanks to Observation 3.2 it is not hard to see that $\Delta_{4}=det(X_{1}^{(1)},X_{2}^{(1)},X_{3}^{(0)},X_{4}^{(1)},X_{5}^{(1)}).$ A computation shows that the latter determinant has order $\mathcal{O}(\lambda_{1}^{2})$, but not order $o(\lambda_{1}^{2})$, which implies that $\lambda_{1}=0.$ A similar analysis leads to the expression of the next term, that is,

[TABLE]

Since the first two rows are linearly dependent, we can see that the determinant of the latter matrix is zero. A more involved analysis is necessary in dealing with the $\varepsilon^{6}$-order term in (3.14): $\Delta_{6}$ can be written as

[TABLE]

A simple analysis shows that the last determinant is zero due to the structure of columns 1, 3, 4 and 5. We end up with

[TABLE]

Setting $\Delta_{6}=0$ and solving for $\lambda_{2}$, we obtain

[TABLE]

This function is clearly positive for all $R>1$, from where we readily conclude inviscid instability for all values of $\gamma\geq 1$, a result that extends and improves the results of [26, Page 3036], which were limited to the case $\gamma\in[1,2]$. Furthermore, our result shows a higher level of accuracy when compared to the results of [26] (see also the Appendix C) and those of [50]; this improvement is clear once we compare the predicted analytically determined value for the instability of the Lopatinsky determinant with those values observed numerically, as discussed in Section 4.

3.1.3. Comparison with previous results: parallel case

We investigate the Lopatinsky determinant and its roots numerically for different values of $\gamma$, corresponding to monoatomic gas ($\gamma=5/3$) and diatomic gas ($\gamma=7/5$; for instance, $O^{2}$) , $\gamma=3$ (artificial gas) and compare those values to the analytically predicted result in equation (3.15); the comparison is shown in table 2.

Another representative description of the good agreement between the analytical result and the numerical study can be also seen in Figures 3 and 1. In particular, Fig. 1 points out the accuracy of our results when compared to those presented in [26]. However, it is worthwhile to stress that the analysis in the latter paper gives the correct order for the root of the Lopatinsky determinant, i.e., $\lambda=\mathcal{O}(\varepsilon^{2})$; apparently, this result was known at a formal level in the Astrophysics community since the late 80’s; see for example [50, §2.1].

3.2. Full inviscid stability diagram

In this section, following the scaling of [6], we fix $u_{1}^{-}=1$, $\rho^{-}=1$ so that $\rho^{-}u_{1}^{-}=\rho^{+}u_{1}^{+}=1$, where the latter is due to the Rankine-Hugoniot conditions (3.6). With regards to the parametrization of Lemma (3.1), it consists of $R=\frac{1}{u_{1}^{+}}$. Note for this choice of parameters that the slow shock classification of (1.3) simplifies to

[TABLE]

We complete our study of parallel inviscid shock stability by a numerical stability analysis over all parameters, complementing the asymptotic study of the previous subsections. Recall [26] that 1-D stability, or nonvanishing of $\Delta(\lambda,1)$, has previously been verified. Thus, without loss of generality, we may fix $\xi=1$ by homogeneity (see Observation 3.2-(iv)), reducing the question of stability to nonvanishing of $\Delta(\lambda,1)$ on $\mathrm{Re}\left(\lambda\right)\geq 0$. Noting that $\Delta(\lambda,1)$ is analytic in $\lambda$, this can be done by a winding number computation. Indeed, if we evaluate $\Delta(\lambda,1)$ along a contour in the complex plane that encloses any possible unstable roots of $\Delta(\lambda,1)$, and if the resulting image contour has winding number 0, then the associated shock is stable, and if the winding number is positive, then the shock is unstable. We numerically evaluate $\Delta(\lambda,1)$ along the contour $\partial(\{z\in\mathbb{C}:\mathrm{Re}(z)\geq 0\}\cap\{z\in\mathbb{C}:|z|=r\}$ where $R>0$ is is sufficiently large. In practice, we took $r=10$, which appears to be amply large. When $\mathrm{Re}(\lambda)=0$, the real part of the eigenvalues of (3.1) collapse to zero making it difficult numerically to detect the correct bases for evaluating $\Delta(\lambda,1)$. To get around this technical difficulty, we simply in practice shift to the right of the contour on which we compute $\Delta(\lambda,1)$ by $1e-4$. As displayed in Fig. 2 (b), we examine inviscid stability for $u_{1}^{+}\in\{0.05,0.1,...,0.9,0.95\}$ and $h_{1}\in\{1.1,1.2,...,3.9,4\}$.

The results displayed in Fig. 2-(b) indicate a single stability transition for each fixed $u_{1}^{+}$ sufficiently close to $u_{1}^{-}=1$ as $h_{1}$ is increased from $H^{*}=1$ (stability) to $\infty$ (instability). For smaller $u_{1}^{+}$, corresponding to larger-amplitude waves, all slow shocks appear to be multi-d unstable independent of the strength of the magnetic field $h_{1}$, hence there is no stability transition. Similarly as in the large-$h_{1}$ case, we observe through a winding number computation of the Evans function that in the unstable case the instability corresponds to a double real root, so that the stability transition as described in the introduction corresponds to passage of a double root through the origin.

3.2.1. The critical destabilization parameter.

Based on the above observations, to pinpoint the location of the stability transition $h_{1}$ for a given fixed $u_{1}^{+}$, we have only to numerically solve $\Delta(0,1)=0$, considered as an equation in $h_{1}$. (Note that, as a consequence of reflection symmetry, $\Delta(\lambda,\xi)$ may be normalized to be real valued for $\lambda\in\mathbb{R}$). Again, to avoid technical difficulties to do with pure imaginary $\lambda$, we solve the approximate equation $\Delta(\lambda_{0},1)=0$ where $\lambda_{0}=1e-5$ via the bisection method, where $\Delta(\lambda,1)$ is normalized by $\Delta(0.1,1)=1$. The result is displayed in Fig. 2(a)(b) together with the results of our more complete coarse-mesh computations. The thick dashed line in Fig. 2 marks the critical destabilization parameter. Apparently the stability region is exclusively determined by this critical destabilization parameter curve.

3.3. Remarks on the nonparallel case

Our approach to the nonparallel case relies on numerical winding number computations, using similar techniques to those used in the parallel case: we first compute an eigenbasis associated to decaying manifolds on both $x_{1}\gtrless 0$ sides, and later expand their entries in $\varepsilon$ in order to obtain an expression for the Lopatinsky determinant $\Delta$ in terms of $\varepsilon$ as in (3.14). One obtains an expansion of the type

[TABLE]

with $\lambda=\lambda_{0}+\lambda_{1}\varepsilon+\mathcal{O}(\varepsilon^{2})$. It is not hard to see that $\Delta_{-6}=\Delta_{-5}=0$, since both corresponding matrices have rows of zeros. As before, in order to find instability we need to show that there exists a solution to $\Delta_{-4}=\Delta_{-4}(\lambda_{0})=0$ such that $Re(\lambda_{0})>0$. The analytical study of this determinant is very complicated, even if we use symbolic computations, so we approach this part numerically: we verify numerically using a winding number computations that the Lopatinsky determinant (3.16) has a simple root $\lambda=\lambda_{0}+\lambda_{1}\varepsilon+\mathcal{O}(\varepsilon^{2})$ such that $\lambda_{0}$ has a positive real part, so we do have instability. This result confirms and elucidates the assertions in [26, Remark 3.5] regarding the nonparallel case. In particular, that the critical eigenvalue is order $\mathcal{O}(1)$ rather than order $\mathcal{O}(\varepsilon^{2})$ as in the parallel case. However, the root we found does not agree with the explicit formula given in [26, Remark 3.5]; as the proof of this formula is not given in [26] we are unable to determine the reason for this discrepancy.

Besides these asymptotics, for some specific shocks we also carried out a numerical Lopatinsky analysis, as one can see from the results displayed in Fig. 9 that in the nonparallel case the roots are not real valued. This appearance of complex roots corresponds to a break of $\mathcal{O}(2)$ symmetry upon linearization (see also [46, 47, 23]).

4. Viscous stability analysis: Evans function

We also study the viscous linear stability analysis using Evans function techniques. Roughly speaking, the study of the Evans function proceeds as follows: consider the system,

[TABLE]

and make the change of coordinates $x_{1}\to x_{1}-st$ to obtain,

[TABLE]

We linearize about a planar traveling wave solution $\bar{u}$, traveling in the direction $x_{1}$ to obtain,

[TABLE]

where $\bar{A}_{j}:=A_{j}(\bar{u})$ and the linearization of $A_{1}(u)\partial_{x_{1}}u$ is $\bar{C}u+\bar{A}_{1}\partial_{x_{1}}u$. Next we take the Fourier transform in the variables $\xi_{2},\ldots,\,\xi_{d}$, and the Laplace transform in $t$ to obtain the eigenvalue problem,

[TABLE]

where $\widetilde{A}_{j}:=\hat{A}_{j}+(\bar{B}_{j1})^{\prime}$, $\bar{B}^{j}:=\bar{B}_{j1}+\bar{B}_{1j}$, $\hat{A}_{1}u:=\bar{A}_{1}u-s\bar{A}^{0}u-dB_{11}(\bar{u})(u,\bar{\partial}_{x_{1}}u)$, and $\hat{A}_{j}u:=\bar{A}_{j}u-dB_{j1}(\bar{u})(u,\bar{\partial}_{x_{1}}u)$ for $j\geq 2$. We obtain our Evans function coefficient matrix by using the flux form, only wherever we see $\widetilde{A}^{\xi}:=\sum_{j=2}^{d}\xi_{j}\widetilde{A}_{j}$, we instead use $i(\bar{A}_{1})^{\prime}-i\bar{C}+\sum_{j=2}^{d}\xi_{j}\widetilde{A}_{j}$ (for further discussion, see [3]).

The Evans function $(\lambda,\xi)\mapsto D(\lambda,\xi)$ consists of a measurement at $x_{1}=0$ of the transversality between the decaying manifolds of the ODE (4.1) below when restricted to the spaces $x_{1}\gtrless 0.$ It is an analytic function of both its parameters whenever $\lambda$ is in the domain of consistent splitting (see Section 2.2.5), which in the present case includes $\{(\xi,\lambda):\Re\lambda\geq 0\}\setminus\{(0,0)\}$. For our study of stability in a channel the relevant values of $\xi$ are $\xi\in\frac{2\pi}{L}\mathbb{Z}$, $L$ being the width of the channel we are studying.

According to the results in [60], the Lopatinsky determinant is (in most cases) a first order approximation of the Evans function in the low frequency regime; consequently, inviscid instability implies viscous instability.555 Indeed this implication holds also in cases such as slow MHD shocks for which the approximation property is not known; see [30, Theorem 2.30] and [30, Lemma 8.3 and Prop. 83]. We note that even for slow MHD shocks, the first-order approximation property holds for generic frequency angles $(\xi,\lambda)$ [30, Lemma 8.3]. The latter implication is our main motivation in the search for zeros of the Evans function, i.e., values of the spectral parameter $\lambda$ such that $D(\lambda,\xi)=0$ for some $\xi\in\frac{2\pi}{L}\mathbb{Z}$.

Apart from constraint issues, substantial new difficulties in going from 1-D to multi-D Evans function computations arise:

1. (i)

Number of equations/parameters: we deal with $5\times 5$ system of equations with downstream/upstream shock conditions. The complexity, both mathematical and numerical, is enormous. Symbolic computations are necessary to create the code without risk of human error. Fortunately, the $\beta$-model allows use of previously tested code (STABLAB; see more in [5]); 2. (ii)

Unexpected issues with Evans function asymptotics related to Eulerian vs Lagrangian coordinates make computation of the multi-D numerical Evans function practically impossible; issues were only recently resolved in this project and, to knowledge of the authors, in only one other project ([33]).

4.1. Computing the profile

To solve for the viscous profile numerically, we cut the domain in half and use a coordinate change to reflect the interval $(-\infty,0]$ to $[0,\infty)$. We use matching conditions at $x=0$ and projective boundary conditions at $x=\infty$ that select the decaying solution. To solve the resulting three point boundary value problem, we use MATLAB’s bvp5c solver with relative and absolute error tolerances set respectively to $10^{-6}$ and $10^{-8}$.

4.2. Computation of decaying manifolds and eigenfunctions

We recall that the Evans function takes the form

[TABLE]

where

[TABLE]

and $W_{1}^{-\infty},...,W_{k}^{-\infty}$ and $W_{k+1}^{+\infty},...,W_{n}^{+\infty}$ form a basis for the solution space of (4.2) that decays as $x\to-\infty$ and as $x\to+\infty$, respectively. If $(\lambda_{0},\xi_{0},v_{0})$ is an eigenvalue, Fourier mode, eigenfunction triple, then $D(\lambda_{0},\xi_{0})=0$ and $v_{0}$ can be expressed as a linear combination of $W_{1}^{-\infty},...,W_{k}^{-\infty}$ when $x\in(-\infty,0]$ and as a linear combination of $W_{k+1}^{+\infty},...,W_{n}^{+\infty}$ when $x\in[0,+\infty)$. Hence, to solve for an eigenvalue $\lambda_{0}$ and eigenfunction $v_{0}$ corresponding to a fixed $\xi_{0}$, we may do as in [34], and solve (4.2) as a boundary value problem with $\lambda$ as a free parameter. At $x=\pm\infty$, we use projective boundary conditions, $P^{\pm\infty}v_{0}(\pm\infty)=0$, which force the projection of $v_{0}$ onto the unstable and stable manifolds at $x=\pm\infty$, respectively, to be zero. The projective boundary conditions at $x=\pm\infty$ provide $k$ and $n-k$ boundary conditions, which leaves one additional boundary condition corresponding to the free parameter. We provide a phase condition, such as $\|v_{0}(0)\|=1$ or a component of $v_{0}(0)$ is unity, which selects an eigenfunction from the family $\{cv_{0}:c\in\mathbb{C},c\neq 0\}$. In practice, to numerically approximate $v_{0}$ we divide the domain into two parts, $(-\infty,0]$ and $[0,+\infty)$, and then perform the change of coordinates $x\to-x$ on $(-\infty,0]$, thus doubling the dimension of the system (4.2) now posed on $[0,+\infty)$. We then pose the boundary value problem on the finite interval $[0,L]$ where $L$ is the truncation value approximating infinity as determined in solving the traveling wave profile. In the end, there is one phase condition given at $x=0$ and $n$ projective boundary conditions given at $x=L$.

To obtain an initial guess for the boundary value problem, we apply STABLAB’s built in root finding capabilities, such as the method of moments or a two-dimensional bisection method using squares in the complex plane, to the Evans function to find a $\tilde{\lambda}_{0}$ which approximates the eigenvalue $\lambda_{0}$ of interest. To approximate the eigenfunction $v_{0}$, we set

[TABLE]

and $W^{R}(x):=[W_{k+1}^{+\infty}(x;\tilde{\lambda}_{0},\xi_{0}),...,W_{n}^{+\infty}(x;\tilde{\lambda}_{0},\xi_{0})]$ and then find $C=(c_{L},c_{R})^{T}$ that minimizes $\|[W^{L}(0),-W^{R}(0)]C\|$ in the least squares sense subject to $\|C\|=1$.

In practice, solving for $W^{L}(x)$ and $W^{R}(x)$ is difficult because of competing modes of $A(x;\lambda,\xi)$ as $x\to\pm\infty$. Thus, we compute $W^{L}(x)$ and $W^{R}(x)$ using the method of continuous orthogonalization of [36]. In this method, we set $W^{L}=\Omega_{L}\alpha_{L}$ and $W^{R}=\Omega_{R}\alpha_{R}$ where $\Omega_{L}$ and $\Omega_{R}$ are orthonormal basis of $W^{L}$ and $W^{R}$ respectively. In particular, as detailed in [36], $\Omega_{L}$ and $\alpha_{L}$ satisfy the well-conditioned ODEs

[TABLE]

Thus, we can solve the ODE for $\Omega_{L}$ and $\Omega_{R}$ and then minimize $\|[\Omega_{L}(0),-\Omega_{R}(0)]C\|$ subject to $\|C\|=1$, then afterward solve for $\alpha_{L}$ and $\alpha_{R}$ by initializing the associated ODE at $x=0$. We note that solving for $\alpha_{L}$ from $x=0$ to $x=-L$, or for $\alpha_{R}$ from $x=0$ to $x=L$, is numerically well posed as error decays in this direction of integration. We then recover $W^{L}$ and $W^{R}$, which provides an initial guess for the boundary value problem described previously.

Fig. 4 exemplifies the applicability of the numerical construction here described: it shows a graph of the real part of the variable $u_{1}$ of the eigenfunction associated to the bifurcating eigenvalue. Notice the loss in the planar structure, which is also pointed out in [46] in the strictly parabolic case for $\mathcal{O}(2)$ steady bifurcations.

4.3. Viscous stability diagram

To determine stability of the viscous shock waves, we compute the Evans function, similar to the Lopatinsky determinant, on a contour $\Omega_{r}:=\partial(\{z\in\mathbb{C}:\mathrm{Re}(z)\geq 0\}\cap\{z\in\mathbb{C}:|z|=r\}$, where $r$ is now chosen by curve fitting the Evans function to within 0.2 relative tolerance of its asymptotic behavior $D(\lambda)\sim C_{1}e^{C_{2}\sqrt{\lambda}}$, indicating that any zeros of the Evans function that may exist lie within $\Omega_{r}$; see Appendix A for more details. However, we do limit $R\leq 128$ for practical reasons since the time to compute $D(\lambda)$ becomes unreasonable for $|\lambda|$ too large. To compute the Evans function, we use the method of continuous orthogonalization ([37]) described in Section 4.2, computed in pseudo-Lagrangian coordinates for better conditioning as described in [4]. We initialize the Evans ODE with a basis that varies analytically in $\lambda$ via the method of Kato [40] as described in [35]. All of these methods are built into the STABLAB platform with which we perform our computations [5].

When $\gamma=5/3$, we sample the winding number for the Fourier coefficient $\xi\in[0.001,0.004,0.007,$ $0.01,0.04,0.07,0.1,0.14,0.17,0.2]$ for various values of $u_{1}^{+}$ and $h_{1}$ and plot the resulting stability diagram in Fig. 2-(a). When $\gamma=7/5$, we obtain the corresponding stability diagram plotted in Fig. 6. We note that the z-axis in Fig. 6 indicates the value of $\xi$ to give a sense of which modes are unstable. In Table 3 we indicate for various parameters the radius needed to enclose any potentially unstable eigenvalues and we indicate the unstable root when it exists.

4.4. On the symmetry of eigenfunctions and equivariance of the Evans function

It is not possible to conclude from our analytical results that the bifurcating eigenvalues have associated dimension 2 rather than $2n$, $n\in\{2,3,\ldots\}.$ On the other hand, as we discuss next, an interesting conclusion can be derived with regards to the symmetry of the decaying manifolds discussed in this section.

By definition, in $\mathcal{O}(2)$ symmetric systems, if $v(x,y)$ is a solution then $\mathcal{R}v(x,-y)$ is also a solution, for $\mathcal{R}$ an orthogonal matrix. Fourier transforming in the y-direction, we observe that any real eigenvalue $\lambda$ has an eigenfunction $e^{iky}w(x)$ and also an eigenfunction $e^{-iky}\mathcal{R}w(x)$. However, by complex symmetry so is $e^{iky}\overline{w}(x)$. Likewise, $e^{iky}\mathcal{R}\overline{w}(x)$. This suggests that an eigenvalue $\lambda=0$ associated with $k_{*}\neq 0$ should have total multiplicity 4, 2 for each of $\pm k_{*}$, unless the apparently non-generic situation occurs that $\mathcal{R}w$ and $\overline{w}$ are constant multiples of one another.

The analogy to complex conjugation is apparent: given any eigenfunction $w$, form $(w+\mathcal{R}\overline{w})$, the real part, and $(1/i)(w-\mathcal{R}\overline{w})$, the imaginary part of the eigenfunction $w$. It is easily seen that both of these are invariant under

[TABLE]

and span the space $\mathrm{span}\{w,\mathcal{R}\overline{w}\}$ contained in the eigenspace of $\mathscr{L}_{k_{*}}$. The same reasoning gives a symmetric basis of the subspaces of decaying solutions of $(\lambda-\mathscr{L})w=0$ at $\pm\infty$. So, we can construct an Evans function from these eigenfunctions, and whenever there is a zero, we can find an eigenfunction given by a real linear combination of them, which is thus itself symmetric under the mapping $\mathcal{T}$, that is, invariance under $\mathcal{T}$ is in fact generic. For there to be non-symmetric eigenfunctions, there would have to be a higher multiplicity of linear dependence. Moreover, MHD gives an explicit example where the multiplicity is in fact 2. Indeed, if there are only 2 eigenvalues, then, choosing the representatives of eigenfunctions having symmetry, we see that this uses up all the dimensions and there cannot be more. Consequently, eigenfunctions may always be chosen with $\mathcal{O}(2)$ symmetry. Furthermore, one can build an Evans function with $\mathcal{O}(2)$ symmetry and this detects “nice” eigenfunctions having the desired symmetry; there may well be others, but this would be “extra”, and there is no reason they would need to be there, that is, they are not generic.

4.5. Finding the critical destabilization parameter $h_{1}$

In the following discussion we assume $\gamma=5/3$ and $u_{1}^{+}=0.86$. Determining exactly where the stability transition occurs in $h_{1}$ for the Evans function is a little difficult because the contour on which we compute the Evans function comes close to a root near the stability transition. A closer approximation was tried, but we could only confirm that, for $\xi=0.005$, the Evans function has a root when $h_{1}=2$, but it does not when $h_{1}=1.999$. Several values of $\xi$ were tested when $h_{1}=1.999$, and no zeros of the Evans function were found; on the other hand, a root is found when $h_{1}=2$, being approximately $2.86\times 10^{-7}$. The Lopatinsky determinant for $h_{1}=2$ and $\xi=1$ has a root at $\lambda_{0}=4.26\times 10^{-4}$, which for $\xi=0.005$, corresponds to $\lambda_{0}=2.13\times 10^{-6}$. It was verified that the Lopatinsky determinant has no root to the right of a the vertical line $\lambda=10^{-4}$ for $h_{1}=1.995$, but it does for $h_{1}=1.996$. The contour cannot be taken much closer to the imaginary axis than $10^{-4}$ because of the essential spectrum.

In summary, we can estimate that the stability transition for the Lopatinsky determinant occurs at $h_{1}\approx 1.995$ and for the Evans function at approximately $h_{1}=2$. As discussed in the introduction, this is due to discretization of Fourier modes. The whole-space transition values agree, as a consequence of concavity of the associated spectral curves, illustrated in Fig. 5.

4.6. Verifying concavity

In the previous subsection, we have verified concavity of the critical spectral curve/agreement of (whole space) viscous and inviscid transition values for one (typical) choice of parameters, essentially by force, by computing the critical spectral curve $\lambda_{*}(\xi)$ and approximating the second derivative. In this subsection, we check concavity more efficiently using the implicit function theorem.

To verify concavity of the spectral curves $\lambda(\xi)=0$ at the critical transition, we approximate the quantity

[TABLE]

using finite difference quotients, where $\tilde{D}(\rho,\lambda_{0},\xi_{0})=D(\rho\xi_{0},\rho\lambda_{0})$ is the Evans function in polar coordinates. This may be recognized as the negative of the “effective viscosity coefficient” of [60, 57, 58, 12], with $\lambda(\xi)=\sigma\xi^{2}+O(\xi^{3})$. Negativity of $\sigma$ corresponds to the “refined stability condition” of the references.

To approximate $\sigma$, we obtain an initial basis for the Evans function ODE by finding a basis $B_{\pm}$ at $(\lambda,\xi)=(0,\varepsilon)$, and then multiplying $B_{\pm}$ on the left by an analytic projection onto the desired subspace, thus creating a locally analytically varying basis determining $D$. We then compute the difference quotient approximation

[TABLE]

for various values of $\varepsilon>0$. We perform a convergence study to verify the correctness of our approximation of $\sigma$, using five evenly spaced values of $\varepsilon$ between 1e-3 and 1e-6. We plot $\sigma$ against $u_{1}^{-}-u_{1}^{+}$ in Fig. 11(a), demonstrating that $\sigma$ is always negative. Recall that $\sigma>0$ indicates that instability may occur in the viscous system before it does in the inviscid system as the bifurcation parameter $h_{1}$ is increased. We also did a spot check to verify that we get the same value for $\sigma$ when we interpolate the curves $\lambda(\xi)=0$ with quadratic polynomials for various values of $h_{1}$, and then interpolate the second derivative of these quadratic polynomials in the variable $h_{1}$ with a quartic polynomial, which we evaluate at $h_{1}=H_{*}$. In Fig. 11(b), we demonstrate the quadratic interpolation of a typical curve $\lambda(\xi)=0$. We note that we compute $\sigma$ only for $u_{1}^{+}$ as small as $0.7005\approx U_{*}$, where $U_{*}$ is the value of $u_{1}^{+}$ at which shocks become inviscid stable at the minimum $h_{1}=H^{*}$ value for which they are 3-shocks.

The clear conclusion from Fig. 11(a) is that the spectral curve is indeed concave at transition to instability, for all relevant values of physical parameters, in the case of a monatomic gas $\gamma=5/3$.

4.6.1. Finer points: analyticity and glancing

The implicit function computation $\sigma=-\tilde{D}_{\rho}/\tilde{D}_{\lambda_{0}}$ of the second-order coefficient in $\lambda(\xi)=\sigma\xi^{2}+\dots$, or, equivalently

[TABLE]

defined by $\tilde{D}(\rho,\lambda_{0}(\rho),1)=0$, presupposes analyticity of $\tilde{D}$ in $(\rho,\lambda_{0})$. As pointed out in [60, 56], analyticity holds away from “glancing points”, defined as frequencies $\lambda_{0}$ for which $A_{1}^{-1}(\lambda_{0}A_{0}+i\xi_{0}A_{2})$ has neutral (i.e., zero real part) eigenvalues possessing a nontrivial Jordan block.

However, for the parallel MHD equations, [45, Lemma 7.2(ii)] specialized to the 2-D case considered here yields for $(\lambda_{0},\xi_{0})=(0,1)$ that there is always a Jordan block of dimension 2, hence an associated square-root singularity in the initializing decaying eigenspaces at both $x\to\pm\infty$, inherited by the manifolds of decaying solutions. This is readily verified by direct computation of the zero-eigenspace of $iA_{1}^{-1}A_{2}$, which may be seen to have geometric multiplicity 2 but algebraic multiplicity 3. Thus, we cannot simply appeal to nonglancing to conclude analyticity: there is always glancing! On the other hand, the fact that both decaying manifolds at $\pm\infty$ have a square-root singularity at $\lambda=0$ implies that the Evans determinant obtained as their exterior product, by a monodromy argument, or simply by composing the two square roots to obtain a linear factor, is analytic, despite the presence of glancing modes. This justifies our computations above.

5. Additional description of numerics

To create Fig. 8 (a), we used the method of moments described in [16] to determine the roots of the Evans function. The moments were computed on the contour $\partial((\{z\in\mathbb{C}:\mathrm{Re}(z)\geq 0\}\cap\{z\in\mathbb{C}:|z|\leq 0.01\})/\{z\in\mathbb{C}:|z|\geq 10^{-4}\})$. To create Fig. 8 (b), we used a forward finite difference scheme to approximate the derivative using the data given in Fig. 8 (a).

Appendix A On winding number computations

The study of unstable modes reduces to studying an eigenvalue equation: at the inviscid level through Lopatinsky determinant, $\Delta(\lambda_{\mbox{lop}},\xi)$; at the viscous level through an Evans function, $D(\lambda_{\mbox{ev}},\xi)$. Due to analyticity of these two objects in their parameters, the search for growing modes corresponds to verifying if, for a fixed $\xi^{*}$ there exists a root $\lambda_{\mbox{lop}}$ of the Lopatinsky determinant and a root $\lambda_{\mbox{ev}}$ of the Evans function in the half space $\mathrm{Re}(z)>0$ of the complex space $\mathbb{C}$. This computation relies then on winding number computations, based upon the argument principle. An example of these computations can be seen in Figures 12(a) and 12(b).

Appendix B Guide to [26] and [15] # 1: persistence of constraint condition, a second proof

We give here a second proof of Proposition 2.7 for the linearized MHD $\beta$-model,666Inherited from the nonlinear version, by linearity of the constraint. following [26, Remark 3.2]. In what follows, we write $\left(\mathcal{Q}_{1}\right)_{x_{i}}$ or $\mathcal{Q}_{1,x_{i}}$ to denote the partial derivative $\partial_{x_{i}}\left(\mathcal{Q}_{1}\right)$ of a quantity $\mathcal{Q}_{1}$ with respect to $x_{i}$. To begin with, we linearize (1.1c) about a shock profile $\overline{\mathcal{V}}=(\overline{\rho},\overline{u}_{1},0,\overline{h}_{1},0)$ (thus, $\overline{u}_{2}=0$, $\overline{h}_{2}=0$), obtaining

[TABLE]

In this section we will only make use of the first of these equations.

[TABLE]

[TABLE]

The linearized Rankine-Hugoniot conditions across a shock at $x_{1}=0$ yield

[TABLE]

We subtract equations (B.2) and (B.3),

[TABLE]

and we differentiate the equation (B.4a) (respectively, (B.4b)) with respect to $t$ (respectively, $x_{2}$), which is legitimate since we are not differentiating in any direction perpendicular to the shock front. We end up with

[TABLE]

where $\left(\overline{u}_{1}^{+}(h_{2}^{+}-Rh_{2}^{-})-h_{2}^{+}\,\overline{u}_{1}^{+}+h_{2}^{-}\,\overline{u}_{1}^{-}\right)_{x_{2}}=0$; thus,

[TABLE]

Now we differentiate the equation (B.1a) (respectively, (B.1b)) with respect to $x_{1}$ (respectively, $x_{2}$), obtaining

[TABLE]

Using that \mathrm{div}(h)\Big{|}_{t=0}=0 and the jump condition at $x_{1}=0$ provided by the Rankine-Hugoniot condition (B.7), we see by the characteristic method that the only possible solution is $\mathrm{div}(h)=0$ for $t>0$.

Appendix C Guide to [26] and [15] # 2: a different approach towards the derivation of Rankine-Hugoniot conditions

In this appendix we consider a different way to derive the Rankine-Hugoniot conditions upon linearization about a planar shock wave. There are several ways to do that, and in the context of MHD equations we refer to [26] and [45]. As pointed out by Majda in [44], the study of stability of planar shocks reduces to a free boundary problem, in which a parameter measuring the deformation of the planar structure of the shock is introduced. The main point of the analysis is to “trade” the latter unknown deformation by introducing a dynamic boundary condition at the linearized shock front. Now that we have given the rough idea we can put this heuristic on solid mathematical ground: assume initially that we are in a 2-D spatial domain. Throughout this section, the usual $l^{2}(\mathbb{R}^{n})$ inner product is written $\langle\cdot,\cdot\rangle$; a subindex $(\cdot)_{a}$ denotes the partial derivative with respect to the variable $a$, for $a\in\{t,x_{1},x_{2}\}$. Let $u(\cdot)$ be a planar traveling wave with speed $s$ solving the following system of equations

[TABLE]

where $u(x_{1}-st)=u^{\pm}$ ($u^{\pm}$ are constants) whenever $x_{1}-st\gtrless 0$ and $f_{0}$, $f_{1}$, $f_{2}\in\mathscr{C}^{\infty}(\mathbb{R};\mathbb{R}^{n})$. The Rankine-Hugoniot conditions are given by

[TABLE]

where $[\cdot]$ denotes the jump across the shock. We consider a perturbation of this system given by a function $u(\cdot)+v(\cdot)$, taking into account also perturbations in the shock front of the form $x_{1}-\phi(x_{2},t)=0$ for $\phi$ sufficiently smooth. In this case the Rankine-Hugoniot conditions are:

[TABLE]

As pointed out in Section 3, we can take s=0. Linearizing the Rankine-Hugoniot conditions about the shock profile $u(\cdot)$ and the shock $x_{1}-st=0$ we have

[TABLE]

where $Df_{1}(\cdot)$ denotes the Jacobian of the mapping $f_{1}$. We start with a trivial linear algebra result:

Claim C.1**.**

We can choose $n-2$ vectors $\{v_{1},\ldots v_{n-2}\}$ in $\mathbb{R}^{n}$ so that

[TABLE]

for all $i\in{1,\ldots,n-2}.$ Further, thanks to (C.1), we must have that $\langle v_{i},[Df_{1}(u)v]\rangle=0$.

Each one of these vectors provide a “static” constraint. It turns out that we can actually find another vector $\widetilde{w}$ - independent of $v_{1},\ldots,v_{n-2}$ - that is also orthogonal to $-\phi_{t}[f_{0}(u)]-\phi_{x_{2}}[f_{2}(u)]$. In this case, however, since $\phi_{t}$ and $\phi_{x_{2}}$ have a dynamic behavior (i.e., time dependence), we must expect the same for $\widetilde{w}$. The idea consists of looking for a vector of the form $\widetilde{w}=r\partial_{t}+s\partial_{x_{2}}$, where $r,s\in\mathbb{R}^{n}$ are unknowns still to be found. Recall that $[f_{0}(u)]$ and $[f_{1}(u)]$ are constant vectors, since $u$ is constant on both sides of the shock. Our aim is to satisfy $\langle r\partial_{t}+s\partial_{x_{2}},-\phi_{t}[f_{0}(u)]-\phi_{x_{2}}[f_{2}(u)]\rangle=0$; expanding the latter, we obtain the equivalent expression

[TABLE]

In order to verify this formula, it suffices to find $r$ and $s$ in $\mathbb{R}^{n}$ such that

[TABLE]

Our next step consists of defining a projection $u\mapsto\mathbb{P}_{v}(u)$ that projects the vector $u\in\mathbb{R}^{n}$ in the space spanned by the space $v\in\mathbb{R}^{n}$. The following properties of this operator, whose proofs we omit, are standard results in linear algebra:

Claim C.2**.**

Given a vector $v\in\mathbb{R}^{n}$, $v\neq 0$ and defining the operator $\mathbb{R}^{n}\ni u\mapsto\mathbb{P}_{v}(u)$, then the following properties are satisfied:

1. i.

$\mathbb{P}_{v}(0)=0$; 2. ii.

$\mathbb{P}_{v}(v)=v$; 3. iii.

$\mathbb{P}_{v}(w)=0\quad\Leftrightarrow\quad v\perp w\quad\Leftrightarrow\quad\mathbb{P}_{w}(v)=0.$

Let $[f_{0}]^{\perp}\in\mathbb{R}^{n}$ (respectively, $[f_{1}]^{\perp}\in\mathbb{R}^{n})$ be defined so that $\{[f_{0}]^{\perp},[f_{0}],v_{1},\ldots,v_{n-2}\}$ (respectively, $\{[f_{1}]^{\perp},[f_{1}],v_{1},\ldots,v_{n-2}\}$) spans $\mathbb{R}^{n}$. Define

[TABLE]

where $\widetilde{\mu}\in\mathbb{R}$ will be defined later. Without loss of generality, we assume that $\mathbb{P}_{[f_{0}]^{\perp}}([f_{2}])\neq 0$ and $\mathbb{P}_{[f_{2}]^{\perp}}([f_{0}])\neq 0$; indeed, using Claim C.2(iii) and the fact that $\{v_{1},\ldots v_{n-2}\}\perp\mathrm{span}\{[f_{0}],[f_{1}]\}$, if one of these projections is vanishes then the $[f_{0}]$ and $[f_{1}]$ are linearly dependent, thus the choice of a $(n-1)^{th}$ orthogonal vector is reduced to a trivial problem. On the other hand, it is easy to see that the last two equations in (C.2) are satisfied. So we proceed as follows: we plug $r$ and $s$ as defined in (C.3) in the first equation of (C.2) to find that $\widetilde{\mu}$ should be

[TABLE]

Applying the result to (C.1) we see that $\displaystyle{\langle r\partial_{t}+s\partial_{x_{2}},-\phi_{t}[f_{0}(u)]-\phi_{x_{2}}[f_{2}(u)]\rangle=0}$, which implies that

[TABLE]

The latter equation is called a dynamic Rankine-Hugoniot condition; notice that it is independent of the perturbation $\phi$ of the shock front.

Remark C.3*.*

A generalization of the results in this appendix to cases with a higher number of vectors is a bit tricky. Indeed, one needs to find more vectors $r_{j,k},1\leq j\leq d-1,1\leq k\leq d-1,$ corresponding to $d-1$ dynamic conditions (this in the generic case that $[f_{0}],[f_{2}],...,[f_{d}]$ are independent, so that there are $n-d$ static conditions and in total there will be the needed $n-1$ total conditions for extreme shock). These must have similar properties to those in (C.2), that is: a) $r_{k}(\lambda,\eta):=\lambda r_{1,i}+i\eta r_{2,k}+$ … $i\eta_{d}r_{d,k},1\leq k\leq d-1$ be independent for each $(\lambda,\eta)\neq 0$, and b) $\langle r_{k}(\lambda,\eta),\lambda[f_{0}]+\eta_{2}[f_{2}]+...\rangle=0.$ Now, the case $d=2$ is easy, because we only need find one of these vectors, and there are $d+d(d-1)/2$ homogeneous constraints, so we can always find a nontrivial solution; to see that it is non-vanishing for all $(\lambda,\eta)$ one just looks and sees a contradiction if one entry but not the other is vanishing. However, when $d\geq 3$, this is not a simple task. Furthermore, there are degenerate cases where $[f_{0}],[f_{2}],...,[f_{d}]$ are not independent, which need to be treated slightly differently also it seems. In conclusion: the method works for the current purposes. In spite of its limitations, it illustrates a case where one can find Rankine-Hugoniot conditions explicitly, even when the shock front is an unknown.

C.1. The dynamic jump condition of [26]

Recall that $u_{2}^{\pm}=0$, $h_{2}^{\pm}=0.$ We make use of equations (1.1) in order to derive the vectors $f_{0}$, $f_{1}$ and $f_{2}$.

[TABLE]

The jumps across the shock are:

[TABLE]

It is easy to find the vectors mentioned in claim C.1: $v_{1}=e_{2}$, $v_{2}=e_{4}$ and $v_{3}=e_{5}$. To calculate the static Rankine-Hugoniot conditions we need the Jacobian of $f_{1}$, $Df_{1}((\rho,u_{1},0,h_{1},0))$:

[TABLE]

Using the notation in [26, §3.2], the result in C.1 provides the static Rankine-Hugoniot conditions. The persistence of the divergence free condition (Section 2 and Appendix B) implies that

[TABLE]

Let $\mathcal{V}=(\rho,u_{1},u_{2},h_{1},h_{2})$. We use $\displaystyle{\langle v_{1},[Df_{1}(\overline{\rho},\overline{u}_{1},0,\overline{h}_{1},0)\mathcal{V})]\rangle=0}$ to get

[TABLE]

The condition

[TABLE]

does not provide anything relevant (since the $4^{th}$ row of $Df_{1}$ is zero), whereas the last static condition, $\langle v_{3},[Df_{1}(\overline{\rho},\overline{u}_{1},0,\overline{h}_{1},0)\mathcal{V})]\rangle=0$ gives

[TABLE]

Remark C.4*.*

If we use the preserved constraint and appropriate normalization (see [26, §3.1]), equation (C.9) yields the first Rankine-Hugoniot condition in [26, Equation (44)], while equation C.11 corresponds to the last equation in [26, Equation (44)].

Now we derive the dynamic Rankine-Hugoniot condition. We begin by calculating the projections mentioned at claim C.2:

[TABLE]

Set $r=[f_{2}]$ and $s=\lambda[f_{0}]$, where $\lambda$ is defined by $\displaystyle{\lambda=-\frac{\langle[f_{2}],[f_{2}]\rangle}{\langle[f_{0}],[f_{0}]\rangle}}.$ The dynamic Rankine-Hugoniot condition will be given by equation (C.4), $\displaystyle{\langle[f_{2}]\partial_{t}+\lambda[f_{0}]\partial_{x_{2}},[Df_{1}(\overline{\rho},\overline{u}_{1},0,\overline{h}_{1},0)\mathcal{V})]\rangle=0}$, i.e.,

[TABLE]

Define the scalars $\mathcal{A}_{1}$ and $\mathcal{A}_{3}$ so that $[f_{0}]=\mathcal{A}_{1}e_{1}$, $[f_{2}]=\mathcal{A}_{3}e_{3}$; it follows that $\mathcal{A}_{1}=\overline{\rho}^{+}(\frac{R-1}{R})$ and, by equation (3.6), $\mathcal{A}_{3}=-\overline{\rho}^{+}(\overline{u}_{1}^{+})^{2}(1-R)$. Clearly, $\lambda=-(\mathcal{A}_{3}/\mathcal{A}_{1})^{2}$. After some computations, the above equation is reduced to

[TABLE]

Following [26, §3.2], we use the normalization $\widetilde{t}\simeq t\overline{u}_{1}^{+}$ to rewrite the previous equation as

[TABLE]

or, equivalently,

[TABLE]

At this point we make some substitutions. First, we solve the Rankine-Hugoniot condition (C.11) for $(u_{2}^{+}-u_{2}^{-})$ and substitute the result in the first term of (C.12). Then we solve (C.9) for ${u}_{1}^{+}$, giving

[TABLE]

and plug the result in the second term of (C.12); note that $\mathcal{G}$ is independent of $u_{1}^{\pm}$ and $\rho^{\pm}$. This gives

[TABLE]

where $\displaystyle{b_{1}=\frac{\overline{u}_{1}^{+}}{2\overline{\rho}^{+}}\left(1+\frac{1}{M^{2}}\right)}$ and $\displaystyle{b_{2}=\frac{\overline{u}_{1}^{+}}{2\overline{\rho}^{+}}R^{2}\left(1+\frac{1}{M^{2}R^{\gamma+1}}\right)}$. Multiplying by $\displaystyle{\frac{M^{2}}{M^{2}-q^{2}}}$, we obtain the dynamic Rankine-Hugoniot condition

[TABLE]

where $\displaystyle{b_{3}=-R(\frac{(M^{2}-1)}{(M^{2}-q^{2})}\rho^{+}}$, $\displaystyle{b_{4}=-\frac{RM^{2}}{(M^{2}-q^{2})}(b_{2}-R)}$, and $\displaystyle{b_{5}=-\frac{M^{2}(R-1)}{M^{2}-q^{2}}}$.

This final equation does not agree with the dynamic Rankine-Hugoniot condition [26, Equation (44)(ii)] of Freistühler & Trakhinin. For, according to their calculations,

[TABLE]

which implies that

[TABLE]

A more careful analysis shows that Equation (C.14) is not a linear combination of (C.13) and the other jump conditions in (C.8)-(C.11) and the interior equations evaluated on $x_{1}=0$. Indeed, if that were the case then $\partial_{\widetilde{t}}\{h_{2}^{-}\}=0$. However there is no way of obtaining $h_{2}^{-}$ as a linear combination of the jump conditions (because it would introduce a linear term in $h_{2}^{+}$) or interior equations evaluated on $x_{1}=0$ (because it would introduce variables that are spatial derivatives in $x_{1}$). One can conclude that the the dynamic Rankine-Hugoniot condition in [26, Equation (44)] is incorrect.

This small error did not do much harm to the main calculations of Freistühler & Trakhinin. In particular they obtained they obtained the correct order for the root of the Lopatinsky determinant (see [26, Equation (61)]), since in the large magnetic field we have $q\to\infty$ and the coefficients of (C.13) and (C.14) only differ by an $\mathcal{O}(\varepsilon^{3})$ term.

Remark C.5*.*

In the $\beta$-model the derivation of the Rankine-Hugoniot conditions using the reasoning presented in this section follows similar lines; the main differences are that

[TABLE]

The jump across the shock are:

[TABLE]

and finally,

[TABLE]

Furthermore, unlike the case presented before, in equation C.10 we would obtain $h_{1}^{+}=h_{1}^{-}$, which we already know for MHD due to persistence of the constraint (1.2) (see Section 2.1.2).