Balanced flux formulations for multidimensional Evans function computations for viscous shocks

TL;DR

This paper introduces a new coordinate formulation for computing the Evans function in multidimensional viscous shock stability analysis, extending the benefits of integrated coordinates from 1D to higher dimensions.

Contribution

It proposes a novel coordinate choice that enables effective Evans function computations in multidimensional settings, where integrated coordinates are not available.

Findings

New coordinate formulation facilitates Evans function computation in multiple dimensions

Extends benefits of integrated coordinates from 1D to multidimensional cases

Enhances stability analysis of viscous shock profiles in higher dimensions

Abstract

The Evans function is a powerful tool for the stability analysis of viscous shock profiles; zeros of this function carry stability information. In the one-dimensional case, it is typical to compute the Evans function using Goodman's integrated coordinates [G1]; this device facilitates the search for zeros of the Evans function by winding number arguments. Although integrated coordinates are not available in the multidimensional case, we show here that there is a choice of coordinates which gives similar advantages.

Equations

$$f^0(U{)}_t+\sum _{j=1}^d{f}^j(U{)}_{x_j}=\sum _{j,k=1}^d(B^{jk}(U)U_{x_k}{)}_{x_j}.$$

$$f^j:{\mathbb{R}}^n\to {\mathbb{R}}^n,j=0,1,\dots ,d.$$

$$B^{jk}(U)=\left(\begin{array}{cc}0& 0\\ 0& b^{jk}(U)\end{array}\right).$$

$$\left(\begin{array}[]{c|c}r\times r&r\times(n-r)\\ \hline\cr(n-r)\times r&(n-r)\times(n-r)\end{array}\right)\,,$$

$$U=\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right),u_1\in {\mathbb{R}}^r,u_2\in {\mathbb{R}}^{(n-r)},$$

$$M=\left(\begin{array}{cc}{M}_{11}& M_{12}\\ M_{21}& M_{22}\end{array}\right),$$

$$U(x,t)=\bar{U}(x_1-st),\underset{z\to \pm \infty }{\lim }\bar{U}(z)=U_{\pm },$$

$$f^{0}(U)_{t}+\big{(}f^{1}(U)-sf^{0}(U)\big{)}_{x_{1}}+\sum_{j=2}^{d}f^{j}(U)_{x_{j}}=\sum_{j,k=1}^{d}(B^{jk}(U)U_{x_{k}})_{x_{j}}\,.$$

$$\det (A_{11}^1(\bar{U})-s{A}_{11}^0(\bar{U}))\ne 0\text{(hyperbolic noncharacteristicity)}$$

$$\sigma \left(\sum {\eta }_j{\eta }_k{b}^{jk}(\bar{U})\right)\ge \theta |\eta {|}^2,\theta >0,\text{for all}\eta =({\eta }_1,\dots ,{\eta }_d)\in {\mathbb{R}}^d\text{(parabolicity).}$$

$$B^{jk}(U)=\left(\begin{array}{cc}0& 0\\ b_{21}^{jk}(U)& b_{22}^{jk}(U)\end{array}\right),$$

$$\det (A_{11}^1-A_{12}^1(b_{22}^{11}{)}^{-1}{b}_{21}^{11}-s{A}_{11}^0)(\bar{U})\ne 0,$$

$${\tilde{f}}^1(\bar{U}{)}^{\prime }=(B^{11}(\bar{U}){\bar{U}}^{\prime }{)}^{\prime },{\tilde{f}}^1(\bar{U}):=f^1(\bar{U})-s{f}^0(\bar{U}).$$

$$B^{11}{U}^{\prime }={\tilde{f}}^1(\bar{U})-{\tilde{f}}^1(U_{-}).$$

$$0$$

$$b^{11}(\bar{U}){\bar{u}}_2^{\prime }$$

$${\tilde{f}}^1(U_{+})-{\tilde{f}}^1(U_{-})=0.$$

$${\bar{A}}^0{U}_t+\sum _{j=1}^d({\bar{A}}^{j}U{)}_{x_j}=\sum _{j,k=1}^d({\bar{B}}^{jk}{U}_{x_k}{)}_{x_j},$$

$${\bar{A}}^{1}U:=A^1(\bar{U})U-s{A}^0(\bar{U})-\mathrm{d}{B}^{11}(\bar{U})(U,{\bar{U}}_{x_1})$$

$${\bar{A}}^0:=A^0(\bar{U}),{\bar{A}}^{j}U:=A^j(\bar{U})U-\mathrm{d}{B}^{j1}(\bar{U})(U,{\bar{U}}_{x_1}),(j\ne 1),{\bar{B}}^{jk}:=B^{jk}(\bar{U}).$$

$$\begin{gathered} \lambda {\bar{A}}^{0}U+({\bar{A}}^{1}U{)}^{\prime }+\sum _{j=2}^d\mathrm{i}{\xi }_j{\bar{A}}^{j}U=({\bar{B}}^{11}{U}^{\prime }{)}^{\prime }+\sum _{k=2}^d(\mathrm{i}{\xi }_k{\bar{B}}^{1k}U{)}^{\prime } \\ +\sum _{j=2}^d\mathrm{i}{\xi }_j{\bar{B}}^{j1}{U}^{\prime }-\sum _{j,k=2}^d{\xi }_j{\xi }_k{\bar{B}}^{jk}U. \end{gathered}$$

$$W^{\prime }=\mathbf{A}(x_1;\lambda ,\xi )W.$$

$$D(\lambda ,\xi ):=\det (W_1^{+},\dots ,W_k^{+},W_{k+1}^{-},\dots W_N^{-}){|}_{x_1=0}.$$

$$\lambda {\bar{A}}^{0}U+({\bar{A}}^{1}U{)}^{\prime }=({\bar{B}}^{11}{U}^{\prime }{)}^{\prime }.$$

$$w:={\bar{A}}^{0}U,W^{\prime }=w,$$

$$\lambda W+{\bar{A}}^1({\bar{A}}^0{)}^{-1}{W}^{\prime }={\bar{B}}^{11}(({\bar{A}}^0{)}^{-1}{W}^{\prime }{)}^{\prime }.$$

$$Z:=\left(\begin{array}{c}W\\ (0,I_{n-r})({\bar{A}}^0{)}^{-1}{W}^{\prime }\end{array}\right).$$

$${\mathbf{A}}_{\mathrm{int}}(x;\lambda )=\left(\begin{array}{ccc}-\lambda {\bar{A}}_{11}^0\mathbf{a}& 0& {\bar{A}}_{12}^0-{\bar{A}}_{11}^0\mathbf{a}{\bar{A}}_{12}^1\\ -\lambda {\bar{A}}_{21}^0\mathbf{a}& 0& {\bar{A}}_{22}^0-{\bar{A}}_{21}^0\mathbf{a}{\bar{A}}_{12}^1\\ -\lambda ({\bar{b}}^{11}{)}^{-1}{\bar{A}}_{21}^1\mathbf{a}& \lambda ({\bar{b}}^{11}{)}^{-1}& ({\bar{b}}^{11}{)}^{-1}({\bar{A}}_{22}^1-{\bar{A}}_{21}^1\mathbf{a}{\bar{A}}_{12}^1)\end{array}\right)$$

$$\lambda W_1=-({\bar{A}}_{11}^1,{\bar{A}}_{12}^1)({\bar{A}}^0{)}^{-1}{W}^{\prime },$$

$$\lambda {\bar{A}}^{0}U=({\bar{B}}^{11}{U}^{\prime }-{\bar{A}}^{1}U{)}^{\prime },$$

$$\mathsf{f}:={\bar{B}}^{11}{U}^{\prime }-{\bar{A}}^{1}U,$$

Full text

B.B. was partially supported by NSF grant DMS-0801745

J.H. was partially supported by NSF grant DMS-0847074 (CAREER)

G.L. was partially supported by NSF grants DMS-0845127 (CAREER) and DMS-1413273

K.Z. was partially supported by NSF grant DMS-0801745

1. 1 Introduction

2. 1.1 Overview

3. 1.2 Plan

4. 2 Preliminaries

5. 2.1 Conservation laws with viscosity

6. 2.2 Standing-wave profiles & the eigenvalue problem

7. 2.2.1 Profile solutions

8. 2.2.2 Linearization, eigenvalue problem

9. 2.2.3 First-order systems

10. 3 Formulating the Evans function

11. 3.1 Flux variables and integrated coordinates ($d=1$)

12. 3.1.1 Integrated coordinates

13. 3.1.2 Flux variables

14. 3.1.3 Balanced flux variables

15. 3.2 Flux, balanced flux, and modified balanced flux variables $(d>1)$

16. 3.2.1 Flux Variables $(d>1)$

17. 3.2.2 Balanced flux variables $(d>1)$

18. 3.2.3 Modified balanced flux variables $(d>1)$

19. 4 Low-frequency behavior of the balanced flux forms

20. 5 Practical Considerations

21. 5.1 Kato bases

22. 5.2 Application in different frequency regimes

23. A Integrated coordinates: $b_{21}\neq 0$ ($d=1$)

1. Introduction

1.1. Overview

The Evans function has proven to be a potent theoretical and numerical tool for the stability analysis of viscous shock profiles; see, e.g., [HLZ1, Z1]. In the multidimensional setting, the Evans function $D$ is a function of frequencies $(\lambda,\xi)\in\{\operatorname{Re}\lambda\geq 0\}\times\mathbb{R}^{d-1}$ where the complex spectral parameter $\lambda$ is dual to time and the vector $\xi$ is dual to the transverse spatial directions. Zeros of $D$ with $\operatorname{Re}\lambda>0$ correspond to perturbations that grow exponentially in time. A central task in the stability analysis of viscous shock profiles is therefore the determination of the number and location of zeros (if any) of $D$ in the unstable half space. Indeed, generalized spectral stability—roughly, the absence of such zeros—is a sufficient condition for nonlinear asymptotic stability with explicit algebraic-in-time rates of decay in $L^{p}$, $p\geq 2$; see [Z1]. The precise statement of generalized spectral stability is formulated in terms of the Evans function itself. For important physical problems, e.g., gas dynamics or magnetohydrodynamics, locating zeros of the Evans function and verifying this Evans-function condition is a task that requires the numerical approximation of $D$.

We discuss here one important practical aspect of computing Evans functions associated with viscous shock profiles for multidimensional systems of conservation laws with physically appropriate “real” or partially parabolic viscosity. In particular, we develop the basic properties of various formulations of the Evans function based on particualr choices of the phase variables in the first-order formulation of the associated eigenvalue problem. We call these the flux, the balanced flux, and the modified balanced flux formulations, and we show that these formulations have concrete benefits for the numerical computation of $D$. The flux coordinates have their origins in the work of Goodman [G1, G2], and the balanced flux formulation was originally introduced by Plaza & Zumbrun [PZ] for the purpose of analyzing the spectral stability of small-amplitude multidimensional relaxation shocks. Here, we propose a further modification of these coordinates that preserves the desirable property of analyticity with respect the complex eigenvalue parameter.

To put these developments in context, we recall that in one space dimension it is a standard practice to use Goodman’s tactic of integrated coordinates [G1]. Importantly, this maneuver removes the translational eigenvalue at the origin and is advantageous both for energy estimates [HZ] and for computation of the Evans function [BHRZ, HLZ1]. We show that our balanced flux coordinates and modified balanced flux coordinates also have this desirable property. Indeed, we give a new, transparent proof of this fact which recovers Zumbrun & Serre’s fundamental link [ZS] between low-frequency behavior of the (viscous) Evans function and the (inviscid) Lopatinskiĭ determinant (see Proposition 4.1 below). Moreover, unlike the balanced flux formulation, our modified version retains analyticity with respect to $\lambda$ while still reducing to the usual integrated Evans function when $\xi=0$ (the one-dimensional case). For each constant $\xi$-slice of frequency space, this is perhaps the truest generalization of the integrated Evans function to the multidimensional setting. However, some of the radial uniformity is lost.

Both of our balanced flux formulations accommodate multidimensional systems with real viscosity and are therefore applicable to physical systems such as gas dynamics or magnetohydrodynamics. Indeed, one significant benefit of the flux framework presented here is that it provides a systematic choice of “good” coordinates for Evans-function computations for the stability of viscous shock profiles. We note that, for example, the coordinates used for the Evans-function computations for one-dimensional gas dynamics in [HLZ1] were created on an ad hoc basis. That is, though they were based on Goodman’s integrated coordinates, the actual construction relied heavily on the precise form of the equations of compressible gas dynamics. Here, under minimal hypotheses, we show—for a broad class of equations in one and several space dimensions—that there is a choice of coordinates that accommodates real viscosity, maintains analyticity with respect to the spectral paramter $\lambda$, and removes the translational eigenvalue at the origin. These features are all important for the practical treatment of shock stability by numerical computation of the Evans function. Finally, the utility of our proposed formulation is demonstrated by numerical computations of the Evans function for multi-dimensions in the fundamental settings of isentropic [BHLZ] and nonisentropic [HLZ2] gas dynamics and MHD [BMZ]; collectively, these are the first successful multi-dimensional Evans-function computations for viscous shock waves.

1.2. Plan

In §2 we establish the setting of our analysis. Namely, we outline the rather general framework of hyperbolic–parabolic systems of conservation laws to which our flux formulations apply. In §3 we describe the balanced flux formulation and its modification. For the convenience of the reader, we describe the full development in both the important special case $d=1$ and in the general case $d>1$. This slight repetition allows us to highlight the connection between the balanced flux form and the oft-used integrated coordinates in one dimension. We establish the main result of the paper—a proof detailing the low-frequency behavior of the flux forms—in §4. Finally, in conclusion, we discuss some practical consequences in §5. Appendix A outlines the generalization of integrated coordinates to the case that the lower-left-hand block of the viscosity matrices does not vanish; see Remark 1.

2. Preliminaries

2.1. Conservation laws with viscosity

A number of physical systems take the form of conservation laws with partially parabolic or “real” viscosity. That is, they are partial differential equations of block hyperbolic–parabolic type with form

[TABLE]

Here, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$, $t\in\mathbb{R}$, and $U\in\mathbb{R}^{n}$ with

[TABLE]

The $d^{2}$ viscosity matrices $B^{jk}\in\mathbb{R}^{n\times n}$ are each assumed to have the block structure

[TABLE]

The blocks in (2.1) have sizes

[TABLE]

and we write $U$ as

[TABLE]

to respect this block structure. We write $A^{j}(U):=\mathrm{d}f^{j}(U)$ for $j=0,1,\ldots,d$, and, when necessary, we write any $n\times n$ matrix $M$ in block form

[TABLE]

with block sizes as in (2.2). We also write $f^{j}_{\ell}$ with $\ell=1$ (or $2$) to denote the first $r$ (or the last $n-r$) component functions of the flux $f^{j}$.

Our interest is in the stability of planar viscous shock profiles. Thus, we consider traveling-wave solutions of the form

[TABLE]

and, given our interest in the stability of these waves, our first step is to transform to moving coordinates $\tilde{x}_{1}=x_{1}-st$ in which $\bar{U}$ becomes stationary. This gives (dropping tildes) the modified system of equations

[TABLE]

We make the structural assumptions

[TABLE]

and

[TABLE]

Remark 1*.*

Our structural conditions apply in complete generality to the principal equations of continuum mechanics: compressible gas dynamics, MHD, and viscoelasticity. The methods described here can be extended to the case that the viscosity matrices have nonzero lower left-hand blocks, i.e.,

[TABLE]

under (H2) and the modified hyperbolic condition

[TABLE]

introducing an “approximate parabolic coordinate” $\check{u}_{2}=u_{2}+b_{22}^{-1}b_{21}u_{2}$ similarly as in [MZ1, MZ2]. This is essential, for example, if there does not exist a true parabolic variable, i.e., $b_{12}\partial_{u_{1}}+b_{22}\partial_{u_{2}}$ is not a matrix multiple of $\nabla_{U}\tilde{u}_{2}(U)$, $\tilde{u}_{2}\in\mathbb{R}^{n-r}$ for some “exact parabolic coordinate” $\tilde{u}_{2}$. It can be useful also if it is more convenient to work with a coordinate other than the true parabolic one. However, in practice we find it more convenient to work with the actual parabolic variable, and so, to simplify the presentation, we will restrict our attention to the main case (2.1), (H1), (H2). We briefly treat the more general case (in one space dimension) in Appendix A.

2.2. Standing-wave profiles & the eigenvalue problem

2.2.1. Profile solutions

Examining (2.6), we see that the standing wave $\bar{U}$ must satisfy the ordinary differential equation (${}^{\prime}=\mathrm{d}/\mathrm{d}x_{1}$)

[TABLE]

Evidently, equation (2.7) can be integrated once to

[TABLE]

Note that, using block structure, we may rewrite (2.8) as

[TABLE]

We expect that the algebraic equation (2.9a) defines a submanifold of $\mathbb{R}^{n}$ on which (2.9b) defines a flow. To solve for $u_{1}$ in terms of $u_{2}$, locally at least, the implicit function theorem requires that $\det\tilde{A}^{1}_{11}(\bar{U})\neq 0$, or, equivalently, (H1). This motivates the introduction of assumption (H1).

An obvious necessary condition for the existence of a traveling-wave connection is that the end states $U_{\pm}$ must be equilibria of (2.8). Therefore, from (2.9), we obtain immediately the Rankine–Hugoniot condition

[TABLE]

Remark 2* (Hyperbolic classification).*

We denote by $i_{{\scriptscriptstyle+}}$ the number of characteristics incoming to the shock from the right and by $i_{{\scriptscriptstyle-}}$ the number of characteristics incoming from the left. We write $i:=i_{{\scriptscriptstyle+}}+i_{{\scriptscriptstyle-}}$ for the total number of incoming characteristics. Then, the hyperbolic classification of $\bar{U}(\cdot)$, i.e., the classification of the associated hyperbolic shock $(U_{{\scriptscriptstyle-}},U_{{\scriptscriptstyle+}})$, is given in the table below.

[TABLE]

2.2.2. Linearization, eigenvalue problem

Supposing, as above, that $U(x,t)=\bar{U}(x_{1})$ is a steady solution of (2.6), we linearize about $\bar{U}$ to obtain an equation describing the approximate evolution of a perturbation also called $U=U(x,t)$. The linearized equations for $U$ read

[TABLE]

where

[TABLE]

and

[TABLE]

Taking the Laplace transform in time (dual variable $\lambda$) and Fourier transform (dual variable $\xi=(\xi_{2},\ldots,\xi_{d})$) in the transverse spatial directions $(x_{2},\dots,x_{d})$, finally, we obtain the generalized eigenvalue equation

[TABLE]

In (2.11) we have now used $U=U(x_{1},\lambda,\xi)$ to represent the transformed perturbation. Now, a basic criterion for stability of the viscous profile $\bar{U}$ is that the eigenvalue equation (2.11) should have no solutions which decay at $x_{1}=\pm\infty$ with $(\lambda,\xi)\in\{\operatorname{Re}\lambda>0\}\times\mathbb{R}^{d-1}$. Searching for such values of $(\lambda,\xi)$ is the spectral stability problem, and the Evans function $D=D(\lambda,\xi)$ vanishes precisely at such values. Our focus, then, is on locating zeros (if any) of $D$. It is clearly advantageous to design an Evans function with as much structure as possible to aid the search for unstable zeros. For example, analyticity is valuable; it allows the search for zeros to proceed using the argument principle.

2.2.3. First-order systems

One may visualize the the construction of the Evans as follows. The basic set-up is based on reformulating the eigenvalue problem (2.11) as a first-order system of differential equations

[TABLE]

We note that our block structure hypotheses imply that $\mathbf{A}$ is an $N\times N$ matrix with $N=(2n-r)$. If $\bar{U}$ decays rapidly to its limiting values $U_{{\scriptscriptstyle\pm}}$ as $x_{1}\to\pm\infty$, then the coefficient matrix $\mathbf{A}$ should also have constant (with respect to $x_{1}$) limiting values. We denote these by $\mathbf{A}^{{\scriptscriptstyle\pm}}(\lambda,\xi)$.

Then, the Evans function is built out of the subspaces of solutions of (2.12) which grow at $-\infty$ and decay at $+\infty$; the construction of these subspaces starts with an analysis of the constant-coefficient limiting system $W^{\prime}=\mathbf{A}^{{\scriptscriptstyle\pm}}(\lambda,\xi)W$. That is, if the collection $\{W_{1}^{{\scriptscriptstyle+}},\ldots,W_{k}^{{\scriptscriptstyle+}}\}$ forms a basis for the solutions of (2.12) that decay at $+\infty$ and, similarly, $\{W_{k+1}^{{\scriptscriptstyle-}},\ldots W_{N}^{{\scriptscriptstyle-}}\}$ spans the solutions that grow at $-\infty$, the Evans function can be written as

[TABLE]

Evidently, if $D(\lambda_{\circ},\xi_{\circ})=0$, then (2.13) shows that there is a linear dependence between these two subspaces. But then there must be a solution which decays at both $\pm\infty$, an eigenfunction.

Remark 3*.*

Clearly, different choices of bases lead to distinct Evans functions, and the Evans function is highly non unique. However, the construction guarantees that each representative chosen from the family of Evans functions has the fundamental property that it vanishes at eigenvalues of (2.11). Indeed, in a companion paper, we discuss how differing coordinate systems at the level of original partial differential equation (2.6) influence the character of the resulting Evans function(s) [BHLZ]. A related issue is the previously mentioned use of integrated coordinates and the ability to manipulate the character of the Evans function through the formulation of the first-order system (2.12). For example, there are several choices of the phase variable $W$. Even though all of these Evans functions carry the same stability information, different versions may be more amenable to analysis or computation in various regimes/settings. For example, when counting zeros by the argument principle, it may be useful to limit excessive winding and unwinding.

3. Formulating the Evans function

3.1. Flux variables and integrated coordinates ($d=1$)

For its independent interest and to showcase the relationship between the flux and balanced flux variables we introduce below and the integrated coordinates that are commonly used in one-dimensional Evans-function calculations, in this subsection we specialize to a single space dimension ($d=1$). We recall that, when formulated in terms of integrated coordinates, the Evans function has the useful property that, for Lax or overcompressive shocks, it does not vanish at the origin. Among other benefits, this feature is useful for the practical computation of the Evans function. While integrated coordinates do not naturally generalize to the multidimensional setting, the flux and balanced flux forms do; we describe this generalization to the case $d>1$ below in Section 3.2.

3.1.1. Integrated coordinates

In the case $d=1$, the linearized equation (2.11) collapses ($\xi=0$), and we may write the associated eigenvalue equation as

[TABLE]

To obtain the integrated Evans function, we define

[TABLE]

and we find, integrating (3.1),

[TABLE]

To write (3.3) as a first-order system, we set

[TABLE]

We thus have $Z^{\prime}=\mathbf{A}_{\mathrm{int}}(x;\lambda)Z$, where (denoting by $\mathbf{a}$ the inverse of the matrix $\bar{A}^{1}_{11}$, recall (H1))

[TABLE]

is obtained by solving for $(I,0)(\bar{A}^{0})^{-1}W^{\prime}$, whence, together with the coordinate $(0,I)(\bar{A}^{0})^{-1}W^{\prime}$, we obtain $(\bar{A}^{0})^{-1}W^{\prime}$ and thus $W^{\prime}$. For this step, multiply (3.3) by $(I,0)$ to obtain

[TABLE]

from which we see that, provided $\bar{A}^{1}_{11}$ is invertible (same assumption needed for flux variables, and indeed even for framing via implicit function theorem of the profile equation; see discussion below (2.9)), we can solve for $(I,0)(\bar{A}^{0})^{-1}W^{\prime}$ in terms of the known coordinates $W$ and $(0,I)(\bar{A}^{0})^{-1}W^{\prime}$ of $Z$. We shall not carry out this computation in detail, as we shall reproduce it by an equivalent and somewhat simpler derivation below.

3.1.2. Flux variables

We now describe an alternative way to write the eigenvalue equation as a first-order system. To write the eigenvalue equation in flux variables, we observe that (3.1) can be rewritten as

[TABLE]

which motivates the definition of the flux variable

[TABLE]

or

[TABLE]

with

[TABLE]

Provided that the $r\times r$ matrix $\bar{A}^{1}_{11}$ is invertible—as assumed in (H1), we see immediately from (3.7a) that

[TABLE]

Thus, using (3.9), we may write (A.3) as a first-order system in flux variables as

[TABLE]

with

[TABLE]

and the coefficient matrix $\mathbf{A}_{\mathsf{f}}$ given by

[TABLE]

Remark 4*.*

Perhaps the quickest route to the form of $\mathbf{A}_{\mathsf{f}}$ comes from multiplying (3.6) from the left by the matrix

[TABLE]

from which we immediately obtain

[TABLE]

Observe that the first row of (3.14) gives (3.9), and the row operation in (3.13) has eliminated $u_{1}$ from the second row.

3.1.3. Balanced flux variables

Introducing the balanced flux variable

[TABLE]

and

[TABLE]

effects a scaling transformation on (3.10); we find in this case that the eigenvalue ODE can be written as

[TABLE]

with

[TABLE]

identical to $\mathbf{A}_{\mathrm{int}}$. That is, the balanced flux and integrated formulations exactly agree.

Remark 5*.*

Noting that the balanced flux variable $\mathsf{f}^{\sharp}$ satisfies

[TABLE]

we see that $(\mathsf{f}^{\sharp})^{\prime}$ and $W^{\prime}$ agree, so that the ODEs must be equivalent. Moreover, keeping in mind the relation (3.19), one may check directly that the two described derivations coincide.

3.2. Flux, balanced flux, and modified balanced flux variables $(d>1)$

3.2.1. Flux Variables $(d>1)$

We now proceed to describe the flux formulation for the multidimensional system (2.6). The starting point is the eigenvalue equation (2.11). First, we note that

[TABLE]

so that we may rearrange (2.11) to

[TABLE]

where $\tilde{A}^{j}:=\bar{A}^{j}+(\bar{B}^{j1})^{\prime}$ and $\bar{B}^{j}:=\bar{B}^{j1}+\bar{B}^{1j}$. Thus, we may define the flux variable $\mathsf{f}$ by

[TABLE]

and the goal is to recast (3.20) as a first-order system. $W^{\prime}=\mathbf{A}_{\mathsf{f}}(x_{1};\lambda,\xi)W$ with

[TABLE]

We write $\bar{B}^{\xi}:=\sum_{j\neq 1}\xi_{j}\bar{B}^{j}$ (and similarly, $\tilde{A}^{\xi}:=\sum_{j\neq 1}\xi_{j}\tilde{A}^{j}$, $\bar{B}^{\xi\xi}$, and so on), and we note that our block structure assumption implies that $\bar{B}^{\xi}$ has the form

[TABLE]

Thus, we may perform a simplifying row operation on (3.21); we multiply on the left by

[TABLE]

Equation (3.21) then becomes

[TABLE]

Evidently, the third row of $\mathbf{A}_{\mathsf{f}}$ can be read off from (3.23). In addition, (3.23) contains the fundamental identity

[TABLE]

which allows us to eliminate $u_{1}$ in favor of $u_{2}$ and $\mathsf{f}_{1}$. To obtain the first two rows of $\mathbf{A}_{\mathsf{f}}$, we write

[TABLE]

in terms of components, so that

[TABLE]

Thus, we may use (3.24) to eliminate $u_{1}$ from (3.26) and obtain $\mathbf{A}_{\mathsf{f}}$. Thus, continuing to denote $(A_{11}^{1})^{-1}$ by $\mathbf{a}$, we find that

[TABLE]

Remark 6*.*

This may readily be seen to be equivalent to the usual Evans function of [Z1, Eq. (3.1), p. 356], based on variable $(u,b_{11}u_{1}^{\prime}+b_{22}u_{2}^{\prime})$, since the phase variables for the two Evans functions are conjugate by a frequency-independent coordinate transformation.

3.2.2. Balanced flux variables $(d>1)$

Proceeding from the flux form, we may, similarly as in §3.1.3 above, define a multi-dimensional Evans function analogous to the integrated Evans function used in one dimension. We define

[TABLE]

with

[TABLE]

to obtain an alternative Evans function. In this case, the first-order system takes the form $\tilde{W}^{\prime}=\mathbf{A}^{\sharp}(x_{1};\lambda,\xi)\tilde{W}$ with $\tilde{W}=(\mathsf{f}^{\sharp},u_{2})^{\top}$ and

[TABLE]

This determines an Evans function $D^{\sharp}(r,\xi^{\sharp},\lambda^{\sharp})$ in the variables $(r,\xi^{\sharp},\lambda^{\sharp})$, from which we may then extract an Evans function

[TABLE]

3.2.3. Modified balanced flux variables $(d>1)$

Alternatively, we may replace $r$ in (3.29) with

[TABLE]

in the above derivation, to obtain an Evans function

[TABLE]

that is analytic in $\lambda$, reducing to the usual (1D) integrated Evans function for $\xi=0$, and still has the desirable property that it is nonvanishing at the origin (where it is now multi-valued, depending on limiting angle). This is perhaps the truest generalization of the integrated Evans function to multi-dimensions, considered $\xi$-slice by $\xi$-slice. However, it loses some uniformity in replacing $|\lambda,\xi|$ by the norm-equivalent (for $\operatorname{Re}\lambda\geq 0$) quantity $|\xi|+\lambda$.

4. Low-frequency behavior of the balanced flux forms

The Evans function $D_{b\mathsf{f}}$ in the phase variables $\tilde{W}=\begin{pmatrix}\mathsf{f}^{\sharp}\\ u_{2}\end{pmatrix}$ is approximately homogeneous near $(0,0)$ and analytic along rays through the origin (equivalently, when written in polar coordinates). As noted above, in one dimension, $(\mathsf{f}^{\sharp})^{\prime}=U$, and so the Evans function $\tilde{D}(\lambda,\xi)$ determined by phase variables $\begin{pmatrix}\mathsf{f}^{\sharp}\\ u_{2}\end{pmatrix}=\begin{pmatrix}\mathsf{f}^{\sharp}\\ (\mathsf{f}^{\sharp}_{2})^{\prime}\end{pmatrix}$ is exactly the usual integrated Evans function, with the desirable property that, for Lax or overcompressive shocks, it does not vanish at the origin. We now show that this desirable property persists also for multi-dimensions. Recall [M, ZS] that inviscid stability of multidimensional shock waves is determined by a Lopatinski determinant $\Delta(\lambda,\xi)$ analogous to the Evans function, defined on $\xi\in\mathbb{R}$, $\operatorname{Re}\lambda\geq 0$; a shock is uniformly inviscid stable if $\Delta\neq 0$ on $\{\operatorname{Re}\lambda\geq 0\}\setminus\{(0,0)\}$.

Proposition 4.1**.**

With appropriately chosen bases at $x_{1}=\pm\infty$,

[TABLE]

for $r:=|\lambda,\xi|$ sufficiently small, where $(\check{\lambda},\check{\xi}):=r^{-1}(\lambda,\xi)$, $\Delta(\lambda,\xi)$ is the inviscid Lopatinski determinant, and $\gamma$ is a transversality coefficient for the traveling-wave ODE that is a constant independent of angle $(\check{\lambda},\check{\xi})$. In particular, for a uniformly inviscid stable shock, $D_{b\mathsf{f}}(\lambda,\xi)$ has a nonvanishing limit $\Delta(\check{\lambda},\check{\xi})$ as $(\lambda,\xi)\to 0$ with $(\check{\lambda},\check{\xi})$ held fixed.

Proof.

This follows from the fact that the component bases of decaying solutions for the standard (“unintegrated”) flux system, when multiplied by ${\rm diag}\{I,|\lambda,\xi|\}$, yield a basis for the balanced flux system, except for “fast” decay elements (corresponding to incoming inviscid modes) which vanish in the $\mathsf{f}$ coordinate and must be treated instead by multiplying by ${\rm diag}\{|\lambda,\xi|^{-1},I\}$ using L’Hôpital’s rule. But, this means the determinant will change by factor of

[TABLE]

where $i$ is the number of incoming and $o$ the number of outgoing hyperbolic modes, and $c=i-o$ is the “degree of compressivity” [ZH, ZS], equal for Lax or overcompressive shocks to the number of zeros at the origin.

Here, we are using the important fact that fast modes may be chosen analytically in $r$ and (locally) in $(\check{\lambda},\check{\xi})$, and in the usual flux variables are independent of angle $(\check{\lambda},\check{\xi})$, being of form

[TABLE]

Thus when multiplied by ${\rm diag}\{|\lambda,\xi|^{-1},I\}$, they transform to form

[TABLE]

which are legitimate basis elements that are still analytic in polar coordinates $(r,\check{\xi},\check{\lambda})$. This validates the choice of fast bases. Likewise, we may check that the first coordinates $\mathsf{f}$ of slow modes are chosen as bases for the associated Lopatinski determinant, so that $\begin{pmatrix}\mathsf{f}\\ 0\end{pmatrix}$ are still independent, and clearly independent of fast modes, so we still have a basis. For details on construction of “fast” and “slow” basis elements in the vicinity of the origin, see [ZH, ZS, Z1].

We thus have $\tilde{D}(\lambda,\xi)\equiv r^{-c}D(\lambda,\xi)$, where $D$ is the usual Evans function as defined, e.g., in [Z1], whereupon (4.1) follows from the fundamental property of the standard Evans function [ZS, Z1], valid for Lax and overcompressive shocks, that

[TABLE]

(For undercompressive shocks $D(\lambda,\xi)=\Delta(\lambda,\xi)+o(r^{c})=r\Delta(\check{\lambda},\check{\xi})+o(r)$, so that the number of zeros at the origin in balanced flux coordinates equals or exceeds the number of zeros in standard flux coordinates.) ∎

Remark 7*.*

The modified balanced flux formulation also removes zeros at the origin, by essentially the same argument, substituting for $|\lambda,\xi|$ the commensurate (for $\operatorname{Re}\lambda\geq 0$) quantity $|\xi|+\lambda$.

5. Practical Considerations

5.1. Kato bases

Proposition 4.1 concerns the Evans function induced by a particular choice of (local) bases near $|\lambda,\xi|=0$, or, equivalently [GZ, BrZ, HuZ, Z1], by the choice of initializing stable/unstable eigenbases $R_{\pm}$ of the limiting coefficient matrices at $x\to\pm\infty$, where $R_{\pm}$ are matrices with columns corresponding to basis elements. In standard practice, this is done using the Kato ODE

[TABLE]

where $P$ is the corresponding stable/unstable eigenprojector [BrZ, HuZ, BHLytZ], and $\dot{}$ denotes variation along a given path in frequency space $(\lambda,\xi)$. It can be shown that this is the unique choice such that $P\dot{R}=0$. When there exist locally (jointly) analytic bases and projectors $V_{\pm}$, $P_{\pm}$ with respect to $(r,\hat{\xi},\hat{\lambda})$, we can write $R_{\pm}=V_{\pm}\alpha_{\pm}$ and use (5.1) to derive a linear analytic-coefficient ODE for coefficients $\alpha_{\pm}$, which are therefore locally (jointly) analytic in $(r,\hat{\xi},\hat{\lambda})$ as well. The change from bases $R_{\pm}$ to $V_{\pm}$ changes the value of the resulting Evans function by a nonvanishing analytic factor, hence the conclusions of Proposition 4.1 remain valid for the standard Kato basis as well. Points where analyticity in $(\hat{\xi},\hat{\lambda})$ fails at $r=0$ correspond (see [ZS, Z1]) to glancing modes for the associated inviscid problem, in which the coefficient matrix possesses a Jordan block. It is shown in [Z1, MeZ] that, at such a point, variations in $r$ and $\lambda$ enter “together,” to lowest order as a linear combination in the lower lefthand corner of the standard Jordan form. For example, a model for a glancing mode/Jordan block of order $2$ is

[TABLE]

where $\tau(\hat{\xi})$ is an analytic function of $\xi$. In this case, writing $\delta:=\hat{\lambda}-i\tau(\hat{\xi})+r$, and making a similar computation, we see that variations enter via a Puissieux series, through $\sqrt{\delta}$. In particular, we find that, both in the inviscid and the viscous problem, the Evans function exhibits a square root-type singularity in $\hat{\lambda}$ at $\hat{\lambda}=i\tau(\hat{\xi})$. This gives a useful check for multi-dimensional viscous computations; see [HLZ2].

5.2. Application in different frequency regimes

By (4.1), the balanced flux form of the Evans function is, with uniform stability, nonvanishing at the origin. This is useful for numerical conditioning in the delicate low-frequency regime. To conveniently check intermediate frequencies by a robust winding number computation, we may instead use the modified balanced flux formulation $D_{mb}$, recovering the desirable property of analyticity in $\lambda$. This reduces to the usual integrated Evans function for $\xi=0$, and still has the desirable property that it is nonvanishing at the origin (now multi-valued, depending on the limiting angle, as is the balanced flux version). For an example of how this works in practice,

Appendix A Integrated coordinates: $b_{21}\neq 0$ ($d=1$)

In the case $d=1$, we indicate the changes incurred by dropping the condition $b^{jk}_{21}=0$ in (2.1). Suppose

[TABLE]

Linearizing about the steady solution $\bar{U}$, we obtain, as before, the linearized system

[TABLE]

and we may write the associated eigenvalue equation as

[TABLE]

In (A.3), we have written

[TABLE]

Defining $w:=\bar{A}^{0}U$, $W^{\prime}=w$, as in Section 3.1.1, and integrating (A.3), we have again (3.3). Setting now

[TABLE]

we obtain by a similar, but more involved, computation to that in the case $b_{21}=0$, $Z^{\prime}=\mathbf{A}_{\mathrm{int}}Z$, where

[TABLE]

Here,

[TABLE]

and

[TABLE]

with

[TABLE]

Note that ($\text{H1}^{\prime}$) implies that $\mathcal{B}$ is invertible. The system is closed by solving for $(I,0)(\bar{A}^{0})^{-1}W^{\prime}$, whence, together with the coordinate

[TABLE]

we obtain $(\bar{A}^{0})^{-1}W^{\prime}$ and thus $W^{\prime}$. For this step, multiply (3.3) by $(I,0)$ to obtain

[TABLE]

from which we see that, provided the modified condition ($\text{H1}^{\prime}$) holds, we can solve for $(I,0)(\bar{A}^{0})^{-1}W^{\prime}$ in terms of the known coordinates $W$ and $(b_{22}^{-1}b_{21},I_{r})(\bar{A}^{0})^{-1}W^{\prime}$.