Blake Barker

Benjamin Melinand
email: meli-nand@ceremade.dauphine.fr.

Kevin Zumbrun
email: kzumbrun@indiana.edu

Kevin Zumbrun Existence

Existence and stability of steady noncharacteristic solutions on a finite interval of full compressible Navier-Stokes equations

Keywords: Steady solutions, gas dynamics, Evans function

Existence and stability of steady noncharacteristic solutions on a finite interval of full compressible Navier-Stokes equations

Introduction

In this paper, continuing investigation in [START_REF] Melinand | Existence and stability of steady compressible Navier-Stokes solutions on a finite interval with noncharacteristic boundary conditions[END_REF] of the isentropic case, we study by a combination of analytical and numerical techniques the existence, uniqueness, and stability of steady solutions of the full (nonisentropic) 1D compressible Navier-Stokes equations on a bounded interval, with noncharacteristic inflow-outflow boundary conditions, and more generally of hyperbolic-parabolic systems of conservation laws of similar abstract type.

This corresponds to the 1D version of the "shock tube" problem of describing flow in a finite length and width channel, with prescribed boundary conditions at the left and right ends. Our main interest is in large-amplitude data. Small-amplitude 1D existence, uniqueness, and spectral stability are shown for general symmetrizable systems in [START_REF] Melinand | Existence and behavior of steady solutions on an interval for general hyperbolic-parabolic systems of conservation laws[END_REF].

As developed in the viscous shock case [BHZ10, BHLZ18a, BHLZ18b, HLZ09, HLZ17], a convenient method to study spectral stability is via numerical Evans function investigations. A useful necessary condition, also based on Evans function considerations, is positivity of the stability index, a mod two count of the Morse index of the linearized operator about the wave. This was trivially evaluable in the isentropic case [START_REF] Melinand | Existence and stability of steady compressible Navier-Stokes solutions on a finite interval with noncharacteristic boundary conditions[END_REF], but is complicated in general. In particular, it does not seem to be analytically evaluable for the nonisentropic case considered here. We carry out here both full Evans function and stability index calculations at the same time, both using the numerical code STABLAB [START_REF] Barker | A MATLAB-based numerical library for evans function computation[END_REF].

Description of main results

Our main analytical result is the global existence of steady solutions of the full polytropic gas equations (2.1) (Theorem 4.7), proved by a Brouwer degree argument using detailed and special ODE estimates, applied to the "Cauchy-to-boundary value" map Ψ defined in Section 2.2. We show, moreover, that global uniqueness of solutions of (2.1) is roughly equivalent to transversality of steady profiles as solutions of the ODE connection problem (2.5)-(2.9). This is equivalent to the nonvanishing of the Jacobian det(dΨ) of (2.9) (Proposition 5.1).

Nonvanishing of det(dΨ) is also seen to be equivalent to nonvanishing of the stability index (Lemma 6.1). Hence a change in sign implies appearance of both nonuniqueness and instability: the usual "exchange of stability" scenario familiar from finite-dimensional ODE. Thus we may study uniqueness in passing, in the course of a larger study of spectral stability.

Augmenting our analytical results for the full polytropic gas equations, we carry out such a study in Section 7 by a systematic numerical Evans function investigation of the "feasible set" C u 0 ,e 0 of profiles realizable by numerical shootings. Our numerical findings are that, on the feasible set C u 0 ,e 0 , the stability index is uniformly positive, indicating uniqueness of large-amplitude solutions, and that steady solutions exhibit uniform spectral stability. We note that nonlinear stability can be shown to follow from spectral stability by similar considerations to those of [MZ19, Section 6]; see [START_REF] Melinand | Existence and behavior of steady solutions on an interval for general hyperbolic-parabolic systems of conservation laws[END_REF].

On the other hand, we show numerically in Section 8.1 that both uniqueness and stability can fail for gas dynamics with an artificial convex equation of state.

Discussion and open problems

The first local existence/uniqueness result for small-amplitude data was established in [START_REF] Kweon | Compressible Navier-Stokes equations in a bounded domain with inflow boundary condition[END_REF], in multi-D. We extend that result here to large-amplitude data in the 1-D case by a combination of analytical (existence) and numerical (uniqueness) investigation, establishing (numerically) time-evolutionary stability as well.

Our findings of global existence and uniqueness for the noncharacteristic problem parallel those of Lions [START_REF] Lions | Mathematical topics in fluid mechanics[END_REF] in the characteristic case u = 0 on the boundary, for which he shows global existence and uniqueness of solutions for arbitrary prescribed average density, in 1-and multi-D. However, they are obtained by quite different techniques, which, moreover, are special to 1D. Indeed, though perhaps intuitively expectable, especially given the uniform shock stability results of [START_REF] Humpherys | Spectral stability of ideal-gas shock layers[END_REF][START_REF] Humpherys | Multidimensional stability of large-amplitude Navier-Stokes shocks[END_REF] for the compressible Navier-Stokes equations, our results of large-amplitude existence, uniqueness and stability are obtained by a combination of exhaustive numerical investigations, and rather delicate degree-theoretic arguments specific to the equations of 1D polytropic gas dynamics under study.

Our investigations of stability belong, rather, to a newer family of investigations blending numerical and analytical techniques to study dynamics and bifurcation of shock waves and related solutions of hyperbolic-parabolic conservation and balance laws, cf. [HLZ09, BFZ15, BHLZ15a, BMZ21, SZ01, Zum10].

In [START_REF] Serre | Boundary layer stability in real vanishing viscosity limit[END_REF][START_REF] Zumbrun | Stability of noncharacteristic boundary layers in the standing-shock limit[END_REF] the case of steady solutions on a half-line was investigated and it was shown that instability of steady solutions can occur, even for the most standard ideal polytropic gas law. This suggests that the question of stability at least is not a foregone conclusion for steady solutions on the interval. Moreover, the nature of instablity found in [START_REF] Serre | Boundary layer stability in real vanishing viscosity limit[END_REF][START_REF] Zumbrun | Stability of noncharacteristic boundary layers in the standing-shock limit[END_REF] involved change of sign in the stability index, which in the present case would signal nonuniqueness as well. On the other hand, our numerical findings (Section 7) indicate that neither of these phenomena in fact occur for polytropic gas dynamics on the interval.

This begs the question whether such detailed and special arguments are necessary, or whether there might instead exist some more straightforward argument for all or part of our results via general principles, such as, e.g., existence of convex entropy as used in Section 4. We give a partial answer to this question in Section 8, exhibiting a counterexample involving an equation of state presented in [START_REF] Barker | Convex entropy, Hopf bifurcation, and viscous and inviscid shock stability[END_REF] for which the equations of compressible gas dynamics possess a convex entropy, but global stability and uniqueness are violated. It is seen that the associated transition to instability can involve either steady bifurcation to multiple solutions, or Hopf bifurcation to time-periodic solutions. The latter phenomenon is significant as the first example of Hopf bifurcation for stationary solutions of compressible gas dynamics, similar to "galloping" or "cellular" instabilities in detonation [START_REF] Texier | Transition to longitudinal instability of detonation waves is generically associated with Hopf bifurcation to time-periodic galloping solutions[END_REF].

It is an interesting question whether our existence result extends to general equations of state considered in [START_REF] Barker | Convex entropy, Hopf bifurcation, and viscous and inviscid shock stability[END_REF]. Note that we obtain nonuniqueness results for a particular equation of states in Section 8. A further very interesting open problem is the extension of our large-amplitude existence results to the true multi-D shock tube problem, generalizing the small-amplitude existence-uniqueness results of [START_REF] Kweon | Compressible Navier-Stokes equations in a bounded domain with inflow boundary condition[END_REF], and the determination of stability of steady multi-D solutions even in the small-amplitude case.

Preliminaries

Equations of motion

The 1D compressible Navier-Stokes equations in Eulerian coordinates are (2.1)

ρ t + (ρu) x = 0 , (ρu) t + (ρu 2 + p) x = αu xx , (ρE) t + (ρuE + pu) x = κT xx + (αuu x ) x where E = e + u 2 2 , p = Γρe, e = c v T ,
and ν = κ cv . Here Γ, c v , ν and α are fixed positive constants; see [Bat99, HLZ09, HLZ17]. As described in [START_REF] Melinand | Existence and stability of steady compressible Navier-Stokes solutions on a finite interval with noncharacteristic boundary conditions[END_REF] in the isentropic case, we seek steady solutions on the interval [0, 1], with noncharacteristic inflow-outflow boundary conditions (2.2) (ρ, u, e)(0) = (ρ 0 , u 0 , e 0 ), (u, e)(1) = (u 1 , e 1 ).

By changing ρ by ρ 0 ρ, u by 1 ρ 0 u, t by ρ 0 t and e by 1 ρ 2 0 e (notice that we can not change x without changing the length of the interval), we assume in the following that (2.3) ρ 0 = 1, u 0 , e 0 , u 1 , e 1 > 0.

Profile equations and formulation as mapping problem

Our main interest is the study of steady solutions, i.e. solutions of (2.4)

(ρu) x = 0 , (ρu 2 + Γρe) x = αu xx , ρu e + u 2 2 + Γρeu x = νe xx + (αuu x ) x
together with (2.2). In order to find these steady solutions, we use a shooting method. Integrating (2.4) from 0 to x and rearranging using (2.3), there exists constants of integration c = (c 1 , c 2 ) to be determined so that we obtain similarly as in [START_REF] Humpherys | Multidimensional stability of large-amplitude Navier-Stokes shocks[END_REF] the profile ODE (2.5)

α u 0 u = c 1 + u + Γ e u , ν u 0 e = c 2 -c 1 u - 1 2 u 2 + e,
together with ρ = u 0 u and with the initial data

(2.6) u(0) = u 0 > 0, e(0) = e 0 > 0.
In this setting

(2.7)

c 1 = α u 0 u (0) -u 0 -Γ e 0 u 0 , c 2 = ν u 0 e (0) + αu (0) -e 0 - 1 2 u 2 0 -Γe 0 ,
where u 0 and e 0 are given and (u (0), e (0)) has to be determined in order to satisfies (u(1), e(1)) = (u 1 , e 1 ).

The domain of the ODE is the set (u, e) ∈ R 2 , u > 0 for which the right hand side of (2.5) is well-defined and from which we can reconstruct ρ. Indeed, we remark that u > 0 is imposed by ρu = constant and ρ > 0. Hence we may ignore the variable ρ in the following. Note also that the physical solutions are the ones for which e is also positive. For a fixed choice of left data (ρ 0 , u 0 , e 0 ) (meaning, by our previous normalization, just a fixed choice of u 0 and e 0 ), we define now the "Cauchy-to-boundary value" mapping

(2.8) Ψ : (c 1 , c 2 ) → (u, e)(1),
where (u, e) denotes the maximal solution of (2.5)-(2.6) for the given value of c = (c 1 , c 2 ). Evidently, solutions of (2.2)-(2.4) thus correspond to solutions of the mapping problem (2.9) Ψ(c) = (u 1 , e 1 ).

The feasible set

In (2.8), we did not specify the domain of c. It is indeed our first order of business to determine it. For a fixed choice of left data (u 0 , e 0 ), we define the feasible set C u 0 ,e 0 as the set of all c for which (2.5)-(2.6) has a continuous solution (u, e) on [0, 1] where u and e are both positive on [0, 1]. Note that C u 0 ,e 0 is not empty since (-u 0 -Γ e 0 u 0 , -(1 + Γ)e 0 -1 2 u 2 0 ) ∈ C u 0 ,e 0 (that corresponds to the constant solution of problem (2.5)-(2.6)). Then, we have the following crucial observation.

Proposition 3.1. The set C u 0 ,e 0 is open and its boundary consists of c for which there exists continuous functions (u, e) on [0, 1], solution of Problem (2.5)-(2.6) on [0, 1), with u, e both positive on [0, 1) and such that e(1) = 0.

Before proving Proposition 3.1, we establish a preliminary result. Lemma 3.2. Let c = (c 1 , c 2 ) ∈ R 2 . Let x * ∈ (0, 1] such that a solution (u, e) of (2.5)-(2.6) is defined on [0, x * ). We have the following statements :

(i) If e > 0 on [0, x * ), then (u, e) is bounded on [0, x * ) uniformly with respect x * and one can extend continuously (u, e) to x * .

(ii) If there exists a constant ẽ > 0, e ≥ ẽ > 0 on [0, x * ), then there exists a constant ũ > 0, u ≥ ũ on [0, x * ).

(iii) If (u, e) → 0 simultaneously as x → x * with (u, e) both positive on [0, x * ), then c 1 < 0 and c 2 < 0.

(iv) Assume c 1 < 0, c 2 < 0, (u, e) is a solution of (2.5)-(2.6) on [0, x * ] and u, e > 0 on [0, x * ]. For any ε > 0, there exists δ > 0 depending only on c 1 , c 2 , ε such that if u(x * ), e(x * ) ≤ δ, there exists x ∈ [x * , x * + ε] such that (u, e) extends continuously as a solution of (2.5)-(2.6) on [0, x) with u, e > 0 on [0, x) and e(x) → 0 as x → x.

Remark 3.3. As we will see in the proof of (iv), one can prove that if c 1 < 0, c 2 < 0 and (u, e) is a solution of (2.5)-(2.6) on [0, 1) with u and e both positive on [0, 1), there exists a constant M > 0 depending only on c 1 , c 2 , Γ, such that for any δ > 0 small enough, for any x * ∈ (0, 1), if 0 < u(x * ), e(x * ) ≤ δ, then 0 < e ≤ δ and 0 < u ≤ M δ on [x * , 1).

Proof of Lemma 3.2. (i) Since e > 0 on [0, x * ), we get 1 2

α u 0 u 2 + ν u 0 e 2 = c 1 u + u 2 + Γe + c 2 e -c 1 eu -eu 2 /2 + e 2 ≤ c 1 u + u 2 + Γe + c 2 e -c 1 eu + e 2 ≤ 1 2 c 2 1 + u 2 + 2u 2 + Γ 2 + e 2 + c 2 2 + e 2 + c 2 1 e 2 + u 2 + e 2 ≤ A α u 0 u 2 + ν u 0 e 2 + B
for some constants A, B > 0 depending only on c 1 , c 2 , Γ, α, ν, u 0 . Hence, |(u, e)| grows at most exponentially, in particular remaining bounded on [0, x * ). Furthermore, u and e can be continuously extended to x * since (u 2 ) and (e 2 ) are bounded and then integrable on [0, x * ).

(ii) The term Γ e u in the u-equation serves as a barrier meaning that there exists u * > 0 such that for any u ∈ (0, u * ], c 1 + u + Γ ẽ u ≥ 0. (iii) Evidently, c 1 < 0, or else u > 0 for u, e > 0, contradicting the assumed convergence to 0. Then, for u > 0 sufficiently small, this implies that -c 1 u -1 2 u 2 > 0 and hence ν u 0 e > c 2 + e. Therefore, c 2 < 0 or else e > 0 for e > 0 and u > 0 sufficiently small, again contradicting convergence. This proves (iii).

(iv) We assume now that c 1 , c 2 < 0, that (u, e) are both positive on [0, x * ] and that u(x * ), e(x * ) ≤ δ. By point (i), one can extend (u, e) as a solution of (2.5)-(2.6) on a interval that strictly contains [0, x * ]. We introduce x = sup{x ≤ x * + ε , (u, e) extends as a solution of (2.5)-(2.6) and are positive on [0, x)} and we keep the notation (u, e) for the solution on [0, x). We first note that so long as u and e remain less than c 2 c 1 -1 on [x * , x) we have ν u 0 e < -1 2 u 2 ≤ 0 and thus e is decreasing on [x * , x). Next, based on the u-equation, several situations can happen around x ∈ [x * , x) :

(a) If u(x) > -c 1 - √ c 2 1 -4Γe(x) 2
, u (x) < 0 and u is decreasing around x.

(b) If u(x) = -c 1 - √ c 2 1 -4Γe(x) 2 , u (x) = 0, u (x) = Γ e (x) u(x) < 0, u is decreasing around x. (c) If u(x) < -c 1 - √ c 2 1 -4Γe(x) 2 , then u < -c 1 - √ c 2 1 -4Γe 2 ≤ 2Γ |e| |c 1 | around x. Therefore, for δ small enough, u ≤ max 1, 2Γ |c 1 | δ and ν u 0 e < c 2 2 on [x * , x
) and e goes to zero as

x → x with |x -x * | ≤ 2ν u 0 |c 2 | δ.
This proves assertion (iv). Thanks to this lemma we can assert that

C u 0 ,e 0 = c ∈ R 2
, where e > 0 on [0, 1] , (u, e) the maximal solution of (2.5)-(2.6) .

Proof of Proposition 3.1. For c = (c 1 , c 2 ), we denote by (u, e) the maximal solution of Problem (2.5)-(2.6). Note that the following map is locally Lipschitz

(3.1) Φ : (u, e) ∈ (u, e) ∈ R 2 , u > 0 → c 1 + u + Γe u , c 2 -c 1 u - 1 2 u 2 + e .
If c ∈ C u 0 ,e 0 , then u and e are defined and positive on [0, 1] and by continuous dependence on parameters of solutions of an ODE, c lies in the interior of C u 0 ,e 0 . In particular C u 0 ,e 0 is open.

We assume in the following that c ∈ C c u 0 ,e 0 . Thanks to Lemma 3.2(i)-(ii), there exists x * ∈ (0, 1] such that u and e are defined and continuous on [0, x * ), u, e > 0 on [0, x * ) and (u, e) can be extended to x * with e(x * ) = 0 and u(x * ) ≥ 0. Our goal is to show that c ∈ ∂C u 0 ,e 0 if and only if x * = 1. Three different situations can then occur. Case (i) : x * = 1 and then e(1) = 0. Case (ii) : u(x * ) > 0 and x * < 1. Then (u, e) is defined on an interval that strictly contains [0, x * ] and e (x * ) ≤ 0. In that case let us show that e must actually cross 0 and must become negative as x crosses x * . Since u(x * ) > 0, (u, e) is defined on an interval that strictly contains [0, x * ] and e (x * ) ≤ 0. Several subcases occur.

Subcase (ii)(a) : If e (x * ) < 0, e crosses 0 and becomes negative as

x crosses x * . Subcase (ii)(b) : If e (x * ) = 0 and u(x * ) = -c 1 , then ν u 0 e (x * ) = -(c 1 + u(x * ))u (x * ) = - u 0 α (c 1 + u(x * )) 2 < 0
and e crosses 0 and becomes negative as x crosses x * . Subcase (ii)(c) : If e (x * ) = 0 and u(x * ) = -c 1 , repeated differentiation shows that derivatives of e and u at x * vanish to all orders. By analyticity of solutions of an analytic ODE (note that u > 0), e ≡ 0 and u ≡ -c 1 , contradicting e(0) = e 0 > 0 so that this subcase can not occur. With such a fact in hand, there exists ε > 0 small enough such that (u, e) is defined on [0, x * + ε], e negative on ]x * , x * + ε] and u positive on [0, x * + ε]. By continuous dependence on parameters that c lies in the interior of C c u 0 ,e 0 and c ∈ ∂C u 0 ,e 0 . Case (iii) : u(x * ) = 0 and x * < 1. In this case, Lemma 3.2(iii) shows that c 1 , c 2 < 0. Let (c 1 , c2 ) close enough to (c 1 , c 2 ) and denote by (ũ, ẽ) the maximal solution of Problem (2.5)-(2.6) associated to (c 1 , c2 ). We then take a δ associated to ε = 1 -x * in Lemma 3.2(iv) that works for any (c 1 , c2 ) close enough to (c 1 , c 2 ). By continuity of u and e, there exists a number µ > 0 small enough such that 0 < u(x * -µ), e(x * -µ) ≤ δ 2 . Then, by continuous dependence on parameters, for any (c

1 , c2 ) close enough to (c 1 , c 2 ), (ũ, ẽ) is defined on [0, x * -µ] and 0 < ũ(x * -µ), ẽ(x * -µ) ≤ δ. Lemma 3.2(iv) shows that there exists x ∈ [x * -µ, 1 -µ], ẽ(x) = 0. In particular, in that case c ∈ ∂C u 0 ,e 0 .
We can now show that Ψ defined in (2.8) is continuous.

Proposition 3.4. The map Ψ is continuous on C u 0 ,e 0 and can be extended to C u 0 ,e 0 as a continuous map denoted again Ψ.

Proof. The fact that Ψ is continuous on C u 0 ,e 0 follows from continuous dependence on parameters of solutions of an ODE (and the fact that the map Φ defined in (3.1) is locally Lipschitz). We consider now c = (c 1 , c 2 ) ∈ ∂C u 0 ,e 0 . Proposition 3.1 shows the maximal solution (u, e) of (2.5)-(2.6) is defined and continuous on [0, 1) and can be extended continuously to 1 with e(1) = 0 and u(1) ≥ 0. Therefore, we can define Ψ(c) = (u(1), 0). If u(1) > 0, (u, e) is defined on an interval that strictly contains [0, 1] and by continuous dependence on parameters, Ψ is continuous at c. We now have to deal with the case u(1) = 0. Lemma 3.2(iii) shows that c 1 , c 2 < 0. Consider ε > 0. Let c ∈ C u 0 ,e 0 close enough to c and denote by (ũ, ẽ) the maximal solution of Problem (2.5)-(2.6) associated to c. By continuity of (u, e) there exists x * ∈ (0, 1) such that 0 < u(x * ), e(x * ) ≤ ε 2 . Then, by continuous dependence on parameters, for any c close enough to c, we have 0 < ũ(x * ), ẽ(x * ) ≤ ε. Finally, by taking ε small enough, Remark 3.3 and the fact that ũ is continuous at 1 give 0 < ũ(1) ≤ M ε (where M depends only on c and Γ). Hence, Ψ is continuous at c.

Existence

We are now ready to study existence. We first show that Ψ is "proper" in the following sense.

Proposition 4.1. Assume that u 0 > 0, e 0 > 0 are fixed. Let c = (c 1 , c 2 ) ∈ R 2 , such that |c| 1 and denote by (u, e) the maximal solution of (2.5)-(2.6). Then, if c ∈ C u 0 ,e 0 , either u(1)

1, e(1)

1 or 0 < u(1) 1.
Proof. Several situations can happen.

Case (i) (c 1 1). If c ∈ C u 0 ,e 0 , α u 0 u ≥ c 1 + u 0 and u(1) ≥ u 0 α (c 1 + u 0 ) + u 0 1. Case (ii) (c 2
-1). We consider the energy y = α 2u 0 u 2 + ν u 0 e. Then,

y = c 2 + (Γ + 1)e + u 2 2 ≤ c 2 + M y
where M is a positive constant depending only on Γ, u 0 , α, ν. Then, for any x in the domain of definition of y, y(x) ≤ y(0)e M x + c 2 M (e M x -1), so that if c 2 < -M e M -1 y(0)e M , e must vanish at a point between 0 and 1 and c ∈ C u 0 ,e 0 . Case (iii) (c 2 1). Let c ∈ C u 0 ,e 0 . We consider again the energy y = α 2u 0 u 2 + ν u 0 e. Then,

y = c 2 + (Γ + 1)e + u 2 2 ≥ c 2 + my
where m is a constant depending only on Γ, u 0 , α, ν. Therefore either e(1)

1 or u(1) 1. Case (iv) (c 1 -1 and c 2 ≤ √ -c 1 ) Let c ∈ C u 0 ,e 0 .
Using again y = α 2u 0 u 2 + ν u 0 e and following case (ii), we get y ≤ -c 1 + M y on [0, 1]. Therefore, for c 1 -1, there exists a constant B depending only on u 0 , e 0 , α, ν, Γ such that e ≤ B √ -c 1 on [0, 1]. Using this fact on the u-equation of System (2.5), we get

α u 0 u ≤ c 1 + u + ΓB √ -c 1 u . Note that for x ∈ [0, 1], if u(x) ∈ -c 1 4 -1 2 1 4 c 2 1 -4ΓB √ -c 1 , -c 1 4 + 1 2 1 4 c 2 1 -4ΓB √ -c 1 , then α u 0 u (x) ≤ 1 2 c 1 and denote b = -c 1 4 -1 2 1 4 c 2 1 -4ΓB √ -c 1 .
We then notice that, for Remark 4.2. We proved in the previous proposition that there exists a constant A > 0 depending only on u 0 , e 0 , α, ν, Γ such that for any

c 2 ≤ -A, c ∈ C u 0 ,e 0 (see case(ii)). Note also that if ν u 0 e 0 + c 2 -c 1 u 0 -1 2 u 2 0 + e 0 ≤ 0 and c 1 + u 0 ≥ 0, then c / ∈ C u 0 ,e 0 . Indeed, in this case, u increases and -c 1 u -1 2 u 2 ≤ -c 1 u 0 -1 2 u 2
0 so that e is decreasing, e ≤ -e 0 and then e crosses 0 in (0, 1].

The previous proposition is not empty in the sense that C u 0 ,e 0 is not bounded.

Lemma 4.3. Assume u 0 > 0, e 0 > 0 are fixed. There exists a positive number A depending only on Γ, u 0 , e 0 , α, ν such that if

c 1 +u 0 + ΓA(|c 2 |+1) u 0 < 0 and c 2 +e 0 > 0, then (c 1 , c 2 ) ∈ C u 0 ,e 0 .
Proof. For c = (c 1 , c 2 ), we denote by (u, e) the maximal solution of Problem (2.5)-(2.6) and by I its interval of definition. Following case (ii) in the previous proposition there exists a constant A > 0 depending only on Γ, u 0 , e 0 , α, ν such that e ≤ A(|c 2 | + 1) on [0, 1] ∩ I. Then we note that since

c 1 < c 1 + u 0 < 0, the map y ∈ [0, u 0 ] → -c 1 y -1 2 y 2 is nonnegative so that for any x ∈ [0, 1] ∩ I such that 0 < u(x) ≤ u 0 , we have ν u 0 e (x) ≥ c 2 + e(x). Note also that for x ∈ [0, 1] ∩ I such that u(x) = u 0 , α u 0 u (x) ≤ c 1 + u 0 + Γ A(|c 2 | + 1) u 0 < 0.
In particular u (0) < 0, e (0) > c 2 + e 0 > 0 and {x ∈ (0, 1] ∩ I , u(x) = u 0 } is empty. By Lemma 3.2(i)-(ii), [0, 1] ⊂ I, e is increasing and 0 < u ≤ u 0 on [0, 1]so that c ∈ C u 0 ,e 0 .

We now define for ε > 0,

E ε = {(x, y) ∈ R 2 , ε < x, y < 1 ε } and Ω ε = Ψ -1 (E ε )
. By continuity of Ψ (Proposition 3.4) and Proposition 4.1, Ω ε is open, bounded and Ω ε ⊂ C u 0 ,e 0 . We denote by Ψ ε the restriction of Ψ to Ω ε , and by d(Ψ ε , Ω ε , (e 1 , u 1 )) the Brouwer degree [Bro11, Mil65, Hir94, Pra06, DM21] of Ψ ε in Ω ε with respect to the target (e 1 , u 1 )).

Recall, for regular values (e 1 , u 1 ), defined as values for which Ψ ε is differentiable and full rank on Ψ -1 ε (e 1 , u 1 ),

(4.1) d(Ψ ε , Ω ε , (e 1 , u 1 )) := c∈Ψ -1 ε (e 1 ,u 1 ) sgn det dΨ ε (c),
that is, the degree counts roots with sign depending on orientation. For arbitrary (not necessarily regular) values (e 1 , u 1 ) for which Ψ -1 ε (u 1 , e 1 ) = ∅, d(Ψ ε , Ω ε , (e 1 , u 1 )) = 0. Thus, nonzero Brouwer degree implies existence of a solution. Finally, recall that d(Ψ ε , Ω ε , (e 1 , u 1 )) is homotopy invariant, so long as Ψ -1

ε (e 1 , u 1 ) remains disjoint from ∂Ω ε . Typically, degree is evaluated at a regular value, then deduced for other values by homotopy invariance.

Corollary 4.4. Assume that u 0 > 0 and e 0 > 0 are fixed. Let u 1 > 0, e 1 > 0. Then for ε > 0 small enough, (u 1 , e 1 ) ∈ Ψ(∂Ω ε ) and the Brouwer degree d(Ψ ε , Ω ε , (e 1 , u 1 )) is independent of (u 1 , e 1 ) and ε.

Proof. Let u 1 > 0, e 1 > 0. First, Proposition 4.1 shows that Ψ -1 (u 1 , e 1 ) is bounded and included in the open set Ω ε for ε small enough. In particular, d(Ψ ε , Ω ε , (e 1 , u 1 )) is independent of ε small enough. Furthermore, we also get from Proposition 4.1 that for any

t ∈ [0, 1], (1 + (1 -t)u 1 , 1 + (1 -t)e 1 ) ∈ Ψ(∂Ω ε ) if ε is small enough. Hence, by homotopy invariance, d(Ψ ε , Ω ε , (e 1 , u 1 )) and d(Ψ ε , Ω ε , (1, 1)) are equal.
At this point we make use of the fundamental property that gas dynamics has an associated convex entropy η(ρ, ρu, ρE) in the sense of [START_REF] Lax | Hyperbolic systems of conservation laws and the mathematical theory of shock waves[END_REF][START_REF] Kawashima | On the normal form of the symmetric hyperbolic-parabolic systems associated with the conservation laws[END_REF], namely η = -ρS, where S(ρ -1 , e) is thermodynamic entropy; see [START_REF] Kawashima | On the normal form of the symmetric hyperbolic-parabolic systems associated with the conservation laws[END_REF]§4]. That is, writing (2.1) as

U t + f (U ) x = (B(U )U x )
x , where U = (ρ, ρu, ρE) t , there hold: for any U ∈ R 3 , the Hessian matrix d 2 η U > 0; dη U • df U = dq U for some flux q; and d 2 η U B(U ) is symmetric positive semidefinite, hence V, d 2 η U B(U )V ≤ 0, with equality if and only if d 2 η U B(U )V = 0, or equivalently B(U )V = 0. Composing the equations on the left by dη, we have

(4.2) η(U ) t + q(U ) x = dη U (B(U )U x ) x = (dη U (B(U )U x )) x -d 2 η U U x , B(U )U x = (dη U (B(U )U x )) x -U x , d 2 η U B(U )U x , giving η(U ) t + q(U ) x -(dη U (B(U )U x ))
x ≤ 0 with equality if and only if B(U )U x = 0. Using the definition B(U )U x = (αu x , κT x + αuu x ) t given by (2.1), we find that B(U )U x = 0 is equivalent to (u, T ) x = 0, and therefore to (u, e) x = 0. Integrating the time-independent profile equation from x = 0 to 1, we thus obtain (4.3) q( Û ) -dη Û (B( Û ) Ûx ) 1 0 ≤ 0, with equality if and only if (û, ê) ≡ constant.

Lemma 4.5. For steady gas dynamics on an interval, (2.1)-(2.2), with constant boundary conditions (u 0 , e 0 ) = (u 1 , e 1 ), the unique global solution is given by the constant solution (ρ, û, ê) ≡ (ρ 0 , u 0 , e 0 ). Equivalently c * = (-u 0 -Γe 0 /u 0 , -(1 + Γ)e 0 -u 2 0 /2) is the unique global solution of Ψ(c) = (u 1 , e 1 ); moreover, it is nondegenerate, with sgn det dΨ(c * ) = +1.

Proof. From ρû =: m ≡ constant, we obtain (ρ, u, e)(0) = (ρ, u, e)(1), or U (0) = U (1) in the vectorial notation above. Note that η may be modified by the addition of any affine function while preserving its properties as a convex entropy (by changing q accordingly). Thus, by an appropriate affine shift, we may arrange that η(U (0)) = 0 and dη U (0) = 0, so that the left-hand side vanishes in (4.3) (since η(U (0)) = η(U (1)) and dη U (0) = dη U (1) ), and therefore (û, ê) ≡ constant. But, then, ρ = m/û ≡ constant as well, and so (ρ, û, ê) ≡ (ρ 0 , u 0 , e 0 ) as claimed. The computation of dΨ(c * ) amounts to integration of a 2 × 2 constant-coefficient linearized equations about this constant solution, hence may be carried out explicitly to find that sgn det dΨ(c * ) = +1. We omit this calculation as we will show it in a simpler and more general way later on. See Subsection 6.1.

Remark 4.6. Lemma 4.5 may be recognized as a particularly concrete instance of results stated for general systems in [MZ21, Thm 2.10] and [MZ21, Prop. 2.9]. In particular, sgn det dΨ(c * ) = +1 is seen by abstract considerations to hold for constant solutions of general symmetrizable systems, without explicit calculation. We note for gas dynamics that the key identity (4.2) may be obtained readily from the thermodynamic relation de = T dS + p dv defining S, where V = 1/ρ, or S t = T -1 (e t -pv t ), together with (2.1), without verifying convexity or symmetrizability, with no need to invoke general theory.

Theorem 4.7 (Large-data existence). For steady gas dynamics on an interval, (2.1)-(2.2), there is at least one steady solution for every choice of left and right data.

Proof. Applying Corollary 4.4, we find that the Brouwer degree is independent of the target (u 1 , e 1 ). Thus we may compute the degree at the constant data (u 1 , e 1 ) = (u 0 , e 0 ). By Lemma 4.5, Ψ -1 (u 1 , e 1 ) consists of the single point c * = (-u 0 -Γe 0 /u 0 , -(1 + Γ)e 0 -u 2 0 /2), at which sgn det dΨ(c * ) = +1. Thus, by (4.1), the degree at (u 1 , e 1 ) is +1. This implies that the Brouwer degree is +1 for all values of the target, implying existence of a solution.

Uniqueness

We next characterize uniqueness, by a global version of the local [BFZ15, Lemma 3.10].

Proposition 5.1. If γ := det dΨ(c) does not vanish on the feasible set C u 0 ,e 0 , then solutions of (2.9) are globally unique for each choice of data (ρ 0 , u 0 , e 0 , u 1 , e 1 ). If on the other hand γ changes sign on the feasible set C u 0 ,e 0 , then even local uniqueness is violated; in particular, there is at least one choice of data possessing multiple solutions.

Proof. Nonvanishing of γ implies nonvanishing of the Jacobian determinant det dΨ(c), which implies dΨ(c) full rank and sgn det dΨ(c) ≡ +1 for all c. In particular all values are regular and it follows that the degree of Ψ with respect to a target (u 1 , e 1 ) is equal to +n, where n is the number of solutions for that data. Since we have already shown that degree is identically equal to +1, this is a contradiction unless roots are unique i.e., n = 1. This proves the first assertion. For the second assertion, just notice that uniqueness implies that degree is equal to the sign of γ at the unique solution and therefore a change of sign in γ implies a change in degree. Hence, by contradiction, uniqueness is impossible when γ changes sign.

Conclusion:

Uniqueness or nonuniqueness hinges on nonvanishing of det dΨ(•) on C u 0 ,e 0 .

Spectral stability and the Evans function

We can reduce Problem (2.1) to Note that for λ = 0 the previous system can be written in the alternative form (6.3) ρu + ûρ = 0 ,

ρ t + (ρu) x = 0 , ρu t + ρuu x + (Γρe) x =
(ρûu + Γρe + Γêρ) x = αu xx , (1 + Γ)ρûe + ρû 2 u x = νe xx + α (ûu x + ûx u) x .
We are using here the standard approach [AGJ90, GZ98] of rewriting (6.1) as a firstorder system and a Cauchy problem. Note that, after eliminating ρ, (6.1) may be rewritten as a first-order system in (u, e, u , e ), following the standard approach of [AGJ90, GZ98], with homogeneous data prescribed on (u, e) at both ends. The Evans function may thus be defined via a "shooting" construction, similarly as in [START_REF] Rousset | Inviscid boundary conditions and stability of viscous boundary layers[END_REF][START_REF] Serre | Boundary layer stability in real vanishing viscosity limit[END_REF] for the half-line case, as (6.4) D(λ) := det u 1 (1) u 2 (1) e 1 (1) e 2 (1) ,

where (ρ j , u j , e j ) are solutions of (6.1) with initial conditions (ρ 1 , u 1 , e 1 , u 1 , e 1 )(0) = (0, 0, 0, 1, 0), (ρ 2 , u 2 , e 2 , u 2 , e 2 )(0) = (0, 0, 0, 0, 1); that is, as the Wronskian at x = 1 of a basis of solutions satisfying the boundary conditions at x = 0. This Wronskian vanishes precisely when there exists a solution vanishing in (e, u) at both x = 0, 1, i.e., an eigenfunction. Evidently, D(•) is analytic in λ on all of C, and real-valued for λ in R, with zeros corresponding to eigenvalues of the linearized operator about the associated steady solution.1

Conclusion: Spectral stability is equivalent to nonvanishing of D on { λ ≥ 0}.

The stability index

Clearly D is real-valued for real λ. It is readily seen (see, e.g. [START_REF] Melinand | Existence and stability of steady compressible Navier-Stokes solutions on a finite interval with noncharacteristic boundary conditions[END_REF]) that D(λ) = 0 for λ real and sufficiently large, hence we may define as in [START_REF] Gardner | The gap lemma and geometric criteria for instability of viscous shock profiles[END_REF] the Stability index (6.5) µ := sgnD(0) lim

λ→+∞ real sgnD(λ)
as a nonvanishing multiple ±sgnD(0) of sgnD(0). Evidently, µ determines the parity of the number of roots of the Evans function with positive real part, or, equivalently (since complex roots occur in conjugate pairs), the number of positive real roots, with +1 corresponding to "even" and -1 to "odd". As such, it is often useful in obtaining instability information. Moreover, we have the following key observation relating the low-frequency stability and the stability index information to transversality of the steady profile solution of the standing-wave ODE.

Lemma 6.1. The zero-frequency limit D(0) is equal to αν u 2 0 multiplied by the Jacobian determinant det dΨ(c) associated with problem (2.9) evaluated at any root c; in particular, (6.6) sgnD(0) = sgn det dΨ(c).

Proof. The proof amounts to the observation that the operations of linearization and integration of the standing-wave ODE commute. Taking the variation of the profile equation (2.5) with respect to c gives (6.7)

α u 0 u = ċ1 + u + Γ ė û - ê û2 u , ν u 0 ė = ċ2 -( ċ1 û + c 1 u) -û u + ė, ( u, ė)(0) = (0, 0),
where ˙denotes variation. Furthermore, we deduce from relations (2.7) that

( ċ1 , ċ2 ) = α u 0 u (0), ν u 0 ė (0) + α u (0) .
It is readily verified for λ = 0 that the eigenvalue equations (6.3) can be integrated from 0 to x to yield the same system (6.7) (note that ρû = u 0 ). Therefore, keeping the notations of (6.4), for ( ċ1 , ċ2 ) = (1, 0), ( u, ė) = u 0 α (u 1 , e 1 ) -u 2 0 ν (u 2 , e 2 ), whereas for ( ċ1 , ċ2 ) = (0, 1), ( u, ė) = u 0 ν (u 2 , e 2 ). The result follows.

Remark 6.2. The previous lemma gives us another way to compute D(0). Considering the problem (6.8)

α u 0 u = d 1 + 1 -Γ ê û2 u + Γ û e , ν u 0 e = d 2 -d 1 û - α u 0 û u + e + Γ ê û u, u(0) = 0, e(0) = 0, we have D(0) = det u 1 (1) u 2 (1) e 1 (1) e 2 (1)
where (u 1 , e 1 ) solves (6.8) for (d 1 , d 2 ) = α u 0 -u 0 -Γ e 0 u 0 , α -e 0 -1 2 u 2 0 -Γe 0 and (u 2 , e 2 ) solves (6.8) for (d 1 , d 2 ) = -u 0 -Γ e 0 u 0 , ν u 0 -e 0 -1 2 u 2 0 -Γe 0 (see (2.7) for the link between (c 1 , c 2 ) and (u (0), e (0)). Remark 6.3. Lemma 6.1 is analogous to the Zumbrun-Serre/Rousset lemmas of [START_REF] Zumbrun | Viscous and inviscid stability of multidimensional planar shock fronts[END_REF][START_REF] Rousset | Inviscid boundary conditions and stability of viscous boundary layers[END_REF] in the whole-and half-line case, which say D(λ) ∼ γδ(λ) for |λ| 1, where γ is a Wronskian encoding transversality of the associated standing-wave ODE and δ is a Lopatinski determinant for the inviscid stability problem (here trivially nonvanishing).

Conclusion:

Both Brouwer degree γ = sgn det dΨ(•) and stability index µ are determined by sgn(D(0)), hence (by Proposition 5.1 and the discussion just above) uniqueness and topological stability information may be obtained by evaluation of D(0) on the feasible set C u 0 ,e 0 . In particular, differently from the cases of the whole-or half-line (see, e.g., the discussion of [Zum01, §6.2]), changes in stability/Morse index associated with passage of a single eigenvalue through λ = 0 are necessarily associated with bifurcation/nonuniqueness.

Numerical investigations

For simple gases, the ratio ν α follows closely to the prediction (7.1)

ν α = 27Γ + 12 16
of statistical mechanics [START_REF] Humpherys | Spectral stability of ideal-gas shock layers[END_REF][START_REF] Humpherys | Multidimensional stability of large-amplitude Navier-Stokes shocks[END_REF]. 2 In our numerics, we will assume, further, (7.1).

Feasible set

For our numerical studies, we rescale equation (2.1) through the following change of coordinates, ρ = ρ 0 ρ, u = u 0 ū, e = u 2 0 ē, t = t u 0 , ᾱ := α ρ 0 u 0 , ν := ν ρ 0 u 0 , which allows us to always fix ρ 0 = u 0 = 1. We note that the assumption concerning the ratio of viscosities for simple gases still holds under this change of coordinates, 16ν = ᾱ(27Γ+12). Hereafter, we drop the bar notation. To map out the feasible set, we solve the profile equation (2.5) (with u 0 = 1) as an initial value problem on the interval [0, 1] with initial conditions (u, e)(0) = (1, e 0 ) for various values of the integration constants c 1 , c 2 . We center the map about the integration constants corresponding to the fixed point,

c 1 = -1 -Γe 0 , c 2 = -1 2 -(1 + Γ)e 0 .
We check to ensure that u and e remain positive throughout the unit interval and that finite blowup does not occur; see Appendix A.1 for details about computational algorithms. In Figure 1, we plot some examples of the feasible set. Note that the feasible set is unbounded (see Lemma 4.3).

We tested the following parameters to see if they lie in the feasibility set,

16

, and

c 1 = -1 -Γe 0 + ∆c 1 , c 2 = -1 2 -(1 + Γ)e 0 + ∆c 2 .

Evans function computations

To numerically compute the Evans function, we use the package STABLAB [START_REF] Barker | A MATLAB-based numerical library for evans function computation[END_REF], which is well tested by this point; for example see [BHLZ15b, BJN + 17, BLZ11]. We provide details about numerical conditioning and algorithm choices we use in STABLAB in Appendix A.2.

Winding number computations

To test for the existence of unstable eigenvalues, we compute the Evans function on a contour consisting of the boundary ∂S of the set S := {z ∈ B(0, 100) : (z) ≥ 0}. We use the functionality built into STABLAB [START_REF] Barker | A MATLAB-based numerical library for evans function computation[END_REF] that adaptively chooses the mesh along ∂S so that the relative error between any two consecutive points on the image of ∂S under the Evans function, C S , varies by no more than 0.2. We then compute the winding number of C S , which is the number of eigenvalues of (6.1) inside S. In Figure 2, we demonstrate the profile and corresponding Evans function computation for representative parameters.

We compute the Evans function on the contour ∂S for the parameters, if they are in the feasible set, given by (Γ, α, e 0 , ∆c 1 , ∆c 2 ) ∈ {2/3, 2/5, 1} × lin(0.1, 2, 10) × lin(0.001, 10, 30) × lin(-50, 50, 50) × lin(-50, 50, 50), where lin(a, b, c) indicates the set containing c evenly spaced points in the interval [a, b], ν = α(27Γ+12)

16

, c 1 = -1 -Γe 0 + ∆c 1 , and c 2 = -1 2 -(1 + Γ)e 0 + ∆c 2 . In all, we computed the Evans function on 670,926 contours, and in all cases found the winding number to be zero. These computations took the equivalent of approximately 83.8 computation days on a desktop with 10 duo cores.

Computations in original coordinates

As mentioned at the beginning of Section 7, we use a convenient scaling for the numerics. To give an idea of the region studied in the original coordinates corresponding to the analytical results, we provide the following plots. In Figure 3 (a), we plot the coordinates for the initial data, and in Figure 3 (b) we plot the coordinates for the final data, for the profiles for which existence is shown. In Figures 3 (c) and (d), we plot the same initial and final profile data for which we were able to compute the Evans function and show stability. We note that computing the Evans function, though straightforward on much of the feasibility region, is very challenging in some parts due to stiffness of the associated ODEs. This stiffness can come from rapid transition of the profile solution associated with differences in initial and final data that are 1; see Figures 3 (a)-(b) for examples of such data. It can also come from points near the feasibility boundary, where the profile nearly blows up to infinity or down to zero on the interior of the computational interval; as one can see from the linearized profile equations (6.8), either of these events leads to blowup of coefficients and associated stiffness of this system of ODE. 17

Global uniqueness/stability index

For the parameters in the feasible set described in Section 7.1, we computed the Evans function at the origin, D(0). We found that the smallest value of D(0) for the computed parameters was 4.19e-4, far greater than the absolute error tolerance 1e-10 of the computation method, as described further in Appendix A.3. Thus, D(0) appears not to vanish on the feasible set; See Figure 4 for a demonstration of how D(0) varies as c 1 and c 2 vary in the feasible set. By Proposition 5.1, uniqueness of profiles is equivalent to nonvanishing of the Jacobian determinant γ associated with the profile ODE, while Lemma 6.1 shows that γ is a real positive multiple of D(0). Thus, the numerically observed nonvanishing of D(0) indicates global uniqueness on the feasible set, as announced in the introduction. Moreover, nonvanishing of D(0) implies further that the sign of the stability index µ, defined in (6.5) as the sign of a nonvanishing real multiple of D(0). is constant. By numerical evaluation at a single choice of profile-alternatively, by continuation to the constant-profile limit, for which the limiting sign of D(λ) at positive real infinity may be computed to be µ = +1-we find that µ ≡ +1 on the feasible set, consistent with (though not implying) spectral stability.

A numerical counterexample

We now consider equations (2.1) subject to the equation of state ē(τ, S) = e S τ + S + τ 2 2 considered in [START_REF] Barker | Convex entropy, Hopf bifurcation, and viscous and inviscid shock stability[END_REF], where τ corresponds to specific volume and S corresponds to entropy.

Specific density is given by ρ = 1 τ , and T = ēS = e S τ + 1, so e S = T -1 ρ , or (8.1)

S = Ŝ(ρ, T ) = ln T -1 ρ .
From this, we obtain

p = p(ρ, T ) = -ē τ = ρ(T -1) - 1 ρ , e = ê(ρ, T ) = T -1 + ln T -1 ρ + 1 2 ρ 2 , (8.2)
closing the system, together with the energy relation E = e + 1 2 u 2 , in terms of variables (ρ, u, T ), T > 1. Alternatively, inverting the relation e = ê(T, ρ) using êT > 0 for T > 1, we may consider it as implicitly determining a system in the usual variables (ρ, u, e), with e > 0. This is the system referred to as the local model in [START_REF] Barker | Convex entropy, Hopf bifurcation, and viscous and inviscid shock stability[END_REF]. Notably, the function η := -ρ Ŝ(ρ, T ), with Ŝ as in (8.1) considered as a function of the conservative variables (ρ, ρu, E) is a convex entropy for system (2.1) in the sense of [START_REF] Lax | Hyperbolic systems of conservation laws and the mathematical theory of shock waves[END_REF][START_REF] Kawashima | On the normal form of the symmetric hyperbolic-parabolic systems associated with the conservation laws[END_REF]; see [START_REF] Barker | Convex entropy, Hopf bifurcation, and viscous and inviscid shock stability[END_REF].

In [START_REF] Barker | Convex entropy, Hopf bifurcation, and viscous and inviscid shock stability[END_REF] it was shown that the local model considered on the whole line has unstable shock waves for parameters for which the inviscid system has stable waves. Here, we demonstrate that the local model considered on a finite interval has parameters for which uniqueness of solutions fails, and also other, nearby parameters for which a Hopf-bifurcation occurs.

These results are guided by the general principles of [START_REF] Zumbrun | Stability of noncharacteristic boundary layers in the standing-shock limit[END_REF] relating spectra of standing shocks on the whole line to spectra of pieces thereof, considered as solutions on a truncated domain. See [MZ21, §3.2] for further discussion in the specific case of a finite interval. The first relevant principle is that spectra on the interval are, for λ ≥ 0 and λ = 0 given in the limit as interval length goes to infinity-equivalently, as viscosity goes to zero-by the direct sum of spectra on the whole line together with spectra of constant boundary layers on the half-line with data corresponding to that on the left (resp. right) endpoint of the interval. This implies that strict instability on the whole line implies strict instability on the interval with associated stability transition as amplitude is increased from a (presumably stable; see [MZ21, Proposition 3.2]) constant steady solution to an unstable one.

The second principle is that in the same large interval length/small viscosity standingshock limit, the stability index does not vanish ([MZ21, Prop. 3.3]), or equivalently D(0) = 0. Thus, if a homotopy is taken from stable constant solutions to unstable standing shock solutions, entirely within the class of standing shocks with sufficiently large interval/small viscosity, then the associated stability transition cannot correspond to a simple crossing of an eigenvalue through the origin λ = 0, as D(0) = 0, and must therefore involve the crossing of one or more pairs of complex conjugate roots, i.e., a Hopf-type scenario.

On the other hand, the first cited principle implies that two of these roots must be near the pair of roots at the origin of the whole-line shock as it undergoes transition to instability: one "translational" eigenvalue fixed at λ = 0 and the crossing eigenvalue corresponding to instability. Thus, we have the picture of a Hopf bifurcation with very nearby roots, i.e., with associated time-period going to infinity, a quite delicate scenario. This makes numerical verification somewhat sensitive; however, it also aids us in finding a more standard bifurcation in the form of a single crossing eigenvalue through λ = 0, as we are able to find by playing with the left and right boundaries of the interval for a given, sufficiently largeamplitude standing shock on the whole line.

Nonuniqueness

Abstract bifurcation result. We first demonstrate (numerically) a bifurcation implying nonuniqueness. Namely, considering the restriction to a finite interval [x L , x R ], of an appropriate standing-shock solution of the local model on the whole line (described in detail below), we show that D(0) changes sign as x L and x R vary; see Figure 7 (a)-(c). Defining by c * (x L , x R ) the value of c corresponding to the shock profile on [x L , x R ], define the map

Φ(c; x L , x R ) := ψ(c * (x L , x R ) + c; x L , x R ) -ψ(c * (x L , x R ); x L , x R ),
where ψ(x; x L , x R ) is the solution map ψ associated with the interval [x L , x R ]. Then Φ(0; x L , x R ) ≡ 0, reflecting the fact that the shock profile restricted to [x L , x R ] solves its own data. Existence of additional roots c = 0 for some x L , x R implies nonuniqueness for the same data. Nonuniqueness may be detected, therefore, using the following abstract bifurcation result, in the spirit of Proposition 5. Proof. Arguing by contradiction, suppose that c = 0 is the unique root of Φ(c; p) = 0 for p in a neighborhood of p * . Thus, Φ does not vanish on the boundary of a small ball B(0, r), hence the topological degree of Φ(•; p) is independent of p. However, at p for which det(dΦ(0; p)) > 0, the degree is by the assumed uniqueness of roots equal to +1, while at points p for which det(dΦ(0; p)) < 0, the degree is -1, a contradiction.

Numerical investigation: methods. To demonstrate non-uniqueness numerically, we first solve approximately for the profile corresponding to the whole-line shock. We then take the piece of that solution on [x L , x R ] as the profile for the finite boundary problem posed on the same interval. Next, following the intuition described at the end of the previous subsection, we seek nonuniqueness by appropriately varying parameters of this "truncated shock" and nearby steady solutions.

The computations showing non-uniqueness are relatively difficult. In the following discussion, S -:= lim x→-∞ S(x), is the left end state value of entropy in the whole-line shock wave solution of the local model. To solve for the profile, we fix the parameters α = κ = 1 and take S -= 1. From the Rankine-Hugoniot conditions, we obtain the other parameters. We then use a boundary value solver to obtain the viscous shock solution.

this computationally favorable regime where M 1 and M 2 are substantially transverse. We plot the resulting profiles corresponding to (ĉ 1 , ĉ2 ) and (c 1 , c2 ) in Figures 6(a)-(b). We note that there is approximately a 20% difference between the upper and lower curves over the interval [-3, 3] depicted, in terms of the ratio of the ≈ 0.2 maximum difference between the two curves to the ≈ 1.0 total variation of each curve, far more than can be attributed to numerical error.

Thus, we conclude that these profiles indeed give an explicit example of nonuniqueness. It may be checked, further, that the two associated parameter pairs as expected lie on opposite sides of the fold curve D(0) ≡ 0, having value D(0) of opposite signs. Indeed, the nullclines of M 1 and M 2 are approximately parallel, and transverse to the fold curve, giving further validation of the numerically observed fold bifurcation scenario. the Evans function computations, we use MATLAB's ode15s with the requested relative and absolute error tolerance set to 1e-10 and 1e-12 respectively. The ode15s routine is a variable-step solver based on variable differentiation formulas of first through fifth orders [START_REF] Shampine | The MATLAB ODE suite[END_REF].

A.4 Computational effort

Computations were done on a desktop with 128GB Ram and a 4.0GHz i7-6950X Intel processor with 25 MB Cache and 10 cores with 20 threads. Computations were done in Matlab using parallel processing. It took 1.37 days of computation to create the data for the final feasibility study figure, Figure 1. It took 83.8 days of computation on all 10 cores to compute the data for the final Evans function figure, Figure 3. Cumulative computations took longer, exceeding five months. Each of the computations in Section 8 took a substantial part of a day to compute. One of the main reasons the computations were time consuming is stiffness of the associated ODE systems. For instance, for the ounterexample of Section 8, continuation of the profile was necessary in order to achieve required accuracy, simple shooting being prohibitively ill-conditioned. Indeed, this project is similar numerically in scope and delicacy to those described in [START_REF] Barker | Convex entropy, Hopf bifurcation, and viscous and inviscid shock stability[END_REF] and [START_REF] Barker | Transverse bifurcation of viscous slow MHD shocks[END_REF], which together represent a new level of computational challenge in numerical Evans function studies.

c

  ∈ C u 0 ,e 0 and c 1 -1, u rapidly goes under b and stays under b. Therefore u(1) ≤ b ≤ 4ΓB √ -c 1 and u(1) 1.

(2

  Γ, α, e 0 , ∆c 1 , ∆c 2 ) ∈ {2/3, 2/5, 1} × lin(0.1, 2, 10) × lin(0.001, 10, 30) × lin(-50, 50, 50) × lin(-50, 50, 50), In the notation of[START_REF] Humpherys | Multidimensional stability of large-amplitude Navier-Stokes shocks[END_REF], α = 2µ + η = 4 3 µ, γ = Γ + 1, and κ cv µ = 9γ-5 4 , giving the result.

Figure 1 :

 1 Figure 1: (a) Plot of the feasible set as e 0 varies when Γ = 1, α = 0.1, and ν = 0.244. (b) Plot of the feasible set with black dots, the set where u goes negative on [0, 1] with blue circles, the set where e goes negative on [0, 1] with green stars, and the set where there is finite time blowup on [0, 1] with red + signs for α = 2, ν = 3.75, Γ = 2/3, e 0 = 2. (c) Plot of the feasible set with black dots, the set where u goes negative on [0, 1] with blue circles, the set where e goes negative on [0, 1] with green stars, and the set where there is finite time blowup on [0, 1] with red + signs for α = 0.2, ν = 1, Γ = 2/3, e 0 = 2. A bold magenta dot marks the constant solution on plots (b) and (c).

Figure 2 :

 2 Figure 2: For the parameters α = 0.1, Γ = 1, ν = 0.2438, e 0 = 0.001, c 1 = -18.35, and c 2 = 0.5184, we plot (a) the boundary layer profile, and (b) the image of S := {z ∈ B(0, 100) : (z) ≥ 0} under the Evans function. The winding number is zero indicating spectral stability of the boundary layer profile.

Figure 3 :

 3 Figure 3: (a) Plot of initial conditions (u 0 , e 0 ) considered. (b) Plot of resulting (u 1 , e 1 ) end conditions when Γ = 2/3 and -50 ≤ c 1 ≤ 50 (approximately 50 points) and -50 ≤ c 2 ≤ 50 (approximately 200 points). In Figures(c) and (d) we plot the profile data for which we were able to compute the Evans function.

  3

Figure 4 :

 4 Figure 4: In these figures, α = 0.7 3, Γ = 2/3, ν = 1.375, and e 0 = 0.001. (a) Plot of D(0) against c 1 and c 2 . (b) Plot of the feasible set corresponding to Figure (a) with black dots, the set where u goes negative on [0, 1] with blue circles, the set where e goes negative on [0, 1] with green stars, and the set where there is finite time blowup on [0, 1] with red + signs.

  1 and [BFZ15, Lemma 3.10]. Proposition 8.1. Let Φ(c; p) : R m × R satisfy Φ(0; p) ≡ 0. If γ := det(dΦ(0; p)) changes sign as p crosses a particular bifurcation value p = p * , then Φ(•; p) has a nontrivial root c = 0 for p arbitrarily close to p * .

Figure 5 :

 5 Figure 5: Figures (a)-(b) demonstrate that D(0) changes sign as c 1 and c 2 vary. Figures (c)-(d) indicate that there are distinct profiles that solve the same data since there are nullclines of M 1 and M 2 that intersect twice. (a) Plot of D(0) against c 1 and c 2 . (b) Plot of sign(D(0)) against c 1 and c 2 . (c) Plot of the nullclines of M 1 and M 2 . Dots indicate intersections of the nullclines. (d) Plot of only the two intersecting nullclines seen in (c).

Figure 6 :Figure 7 :

 67 Figure 6: (a) Plot of the two profiles, solving the same data, against x. The solid blue curves and dashed red curves correspond to the profiles with c 1 and c 2 values plotted as dots with the same colors in Figure 5(c)-(d). (b) Zoomed-in picture of (a) near x = x R .

  αu xx , ρe t + ρue x + Γρeu x = νe xx + αu 2

		12
	from which we obtain the eigenvalue problem around a steady state (ρ, û, ê)
		λρ + (ρu + ûρ) x = 0 ,
	(6.1)	λρu + (ρûu + Γρe + Γêρ) x + ûx (ρu + ûρ) = αu xx ,
		λρe + (ρûe) x + êx (ρu + ûρ) + Γρêu x + Γû x (ρe + êρ) = νe xx + 2αû x u x ,
	with boundary conditions
	(6.2)	(ρ, u, e)(0) = 0, (u, e)(1) = 0.
		x ,

Indeed, as standard in Evans function theory, zeros correspond in both location and multiplicity to eigenvalues of the linearized operator about the wave; see, e.g.,[START_REF] Alexander | A topological invariant arising in the stability analysis of travelling waves[END_REF][START_REF] Gardner | The gap lemma and geometric criteria for instability of viscous shock profiles[END_REF][START_REF] Zumbrun | Errata to: "Pointwise semigroup methods, and stability of viscous shock waves[END_REF] in the whole-line case.

That µ = +1 for constant profiles may be deduced also by the fact [MZ21, Proposition 3.2] that constant steady solutions of general systems with convex entropy are stable.

Research of B.B. was partially supported under NSF grant no. DMS-140087. Research of K.Z. was partially supported under NSF grant no. DMS-1700279.

1

Next, we use continuation with 30 evenly spaced steps in S -to obtain the solution at S -= -5. That is, we change the parameter S -by a small amount and solve for the other parameters given by the Rankine-Hugoniot conditions, then use the profile solution corresponding to the previous value of S -as an initial guess in the boundary value solver to solve for the profile for the new parameters. In solving for the whole-line profile, we use STABLAB which adaptively increases the spatial domain [-L, L], L 1, until the profile converges to the fixed-point end states, corresponding to the shock at x = ±∞, to within requested tolerance.

To compute the Evans function, we use the same procedure as described in Section 7.2, except that we evaluate the Wronskian to obtain the Evans function at x = 0 instead of (x L + x R )/2, and we use "pseudo-Lagrangian coordinates" as described in [START_REF] Barker | Euler versus Lagrange: the role of coordinates in practical Evans-function computations[END_REF] to reduce winding in our winding number studies without changing the zeros of the Evans function. For algorithm details, see Appendix A.3.

Nonuniqueness. To demonstrate abstract non-uniqueness of profile solutions, we take a piece of the whole-line shock for an unstable wave in the local model, and truncate it to a finite interval as described just above. By varying the boundary on the left of this finite interval, we are able to observe a change of sign of the Evans function evaluated at the origin, D(0), indicating by Proposition 8.1 that non-uniqueness of solutions occurs.

Multiple solutions. To find explicit nonunique profiles satisfying the same data, we investigate further, considering not only translated pieces of the unstable whole line shock, but also other steady solutions nearby. Namely, fixing the interval to be [x L , x R ] = [-33.17, 2.9], we compute the Evans function at the origin for profiles with varying c 1 and c 2 to find regions in c 1 and c 2 for which D(0) has opposite sign; see Figures 5(a)-(b). As clearly evident in Figure 5(b), the results are consistent with a fold bifurcation of the map ψ, with orientation sgn det dψ changing across a smooth curve in (c 1 , c 2 ), in which case we may expect an open set of distinct parameter pairs (ĉ 1 , ĉ2 ) and (c 1 , c2 ) corresponding to different profiles solving the same data, lying in regions for which D(0) has opposite sign.

A nice way to numerically find parameter pairs corresponding to two distinct profiles solving the same data is to look at nullclines of the mappings

Here (c * 1 , c * 2 ) are fixed constants of integration that correspond to the whole line shock, which constants of integration we find by solving for them in the Rankine-Hugoniot equation. The other terms used in defining M 1 and M 2 , that is u L (c 1 , c 2 ) and T L (c 1 , c 2 ), are the components of the profiles evaluated at x = x L . We note that these profiles have the same data at x = x R as the the profile corresponding to

. The level sets of M 1 and M 2 intersect in two locations, which we name (ĉ 1 , ĉ2 ) and (c 1 , c2 ), along the same curves, indicating that these constants of integration correspond to two distinct profiles solving the same data; see Figure 5 (c)-(d).

Near the fold curve along which D(0) = 0, the nullclines of M 1 and M 2 are by necessity parallel, making computations of their intersection delicate. Our strategy therefore is to move well away from the D(0) = 0 curve into negative and positive orientation parts of the plane by varying (c * 1 , c * 2 ), then to determine the M 1 , M 2 intersections numerically in

A Documentation of STABLAB

In this appendix, we describe additional computational details geared toward the reader interested in reproducing results. In particular, we provide some references regarding the MATLAB-based package STABLAB that we used extensively throughout this paper. STA-BLAB [START_REF] Barker | A MATLAB-based numerical library for evans function computation[END_REF] is a well-tested package for studying stability of traveling waves using the Evans function. This package has been successfully used in a variety of studies; for example see [Bar20, BHLZ15b, BHLL18, BJN + 17, BLZ11, GLM20, Lyt17, FS18]. For an overview of the methods used in STABLAB, please see [START_REF] Barker | Evans function computation[END_REF][START_REF] Barker | Evans function computation for the stability of travelling waves[END_REF].

A.1 Details of feasibility study

To verify the correctness of our code when computing the feasibility set, we independently coded by hand a constant step-size Runge-Kutta-Fehlberg four to fifth order scheme and compared it to the solution we obtained using standard suite software in MATLAB. For improved accuracy, for the large scale study we use MATLAB's ode15s [START_REF] Shampine | The MATLAB ODE suite[END_REF] routine which is an adaptive step, stiff ODE solver. The solver warnings alert us to finite blowup, and testing the sign of a solution tells us whether or not u and e remain positive throughout the unit interval. In Figure 1, we plot some examples of the feasible set. Note that the feasible set is unbounded (see Lemma 4.3).

A.2 Details of the Evans function computations

We now provide details about numerical conditioning and algorithm choices for the Evans function computations. For background regarding the methods mentioned, please see [START_REF] Barker | Evans function computation for the stability of travelling waves[END_REF]. To compute the Evans function, we use the the flux coordinates described in Section 3.1 of [START_REF] Barker | Balanced flux formulations for multidimensional Evans-function computations for viscous shocks[END_REF], which is equivalent to computing with coordinates (ρ, u, e, u , e ) as described in Section 6. These coordinates are important to use in practice in order to reduce the variation in the image of the Evans function. To improve numerical conditioning of the computation, we evaluate the Evans function wronskian at x = 1/2 with ODE solutions given in the definition of the Evans function initialized at x = 0 with {(0, 1, 0, 0, 0) T , (0, 0, 1, 0, 0) T } and at x = 1 with {(1, 0, 0, 0, 0) T , (0, 1, 0, 0, 0) T , (0, 0, 1, 0, 0)} to recover, as given by Abel's Theorem, a non-vanishing multiple of the Evans function.

We also use the method of continuous orthogonalization [START_REF] Humpherys | An efficient shooting algorithm for Evans function calculations in large systems[END_REF] without the radial equation using Drury's method [START_REF] Drury | Numerical solution of Orr-Sommerfeld-type equations[END_REF] in order to compute the ODE solution, which resolves computational challenges due to differing growth modes. To verify the correctness of our code, we compute D(0) with the radial equation by initializing the ODE solutions at x = 0 only and evolving them to take the determinant at x = 1 with the initializing basis there, and check that this matches the value of D(0) computed with the definition given in (6.4).

A.3 Details of the non-uniqueness study

The boundary value solver we refer to in Section 8.1 is MATLAB's routine bvp5c, which uses a four-stage Lobatto IIIa formula [START_REF] Kierzenka | A bvp solver that controls residual and error[END_REF]. We set the tolerance in bvp5c to 1e-6. For