M C Roubaud

J M Lina

J Carrier

B Torrésani

SPACE-TIME EXTENSION OF THE MEM APPROACH FOR ELECTROMAGNETIC NEUROIMAGING

Keywords: MEG inverse problem, maximum entropy on the mean, wavelet decomposition, spatial filters, Kronecker covariance factorization, sleep slow waves

The wavelet Maximum Entropy on the Mean (wMEM) approach to the MEG inverse problem is revisited and extended to infer brain activity from full space-time data. The resulting dimensionality increase is tackled using a collection of techniques, that includes time and space dimension reduction (using respectively wavelet and spatial filter based reductions), Kronecker product modeling for covariance matrices, and numerical manipulation of the free energy directly in matrix form. This leads to a smooth numerical optimization problem of reasonable dimension, solved using standard approaches.

The method is applied to the MEG inverse problem. Results of a simulation study in the context of slow wave localization from sleep MEG data are presented and discussed.

INTRODUCTION

EEG and MEG inverse problems are notoriously difficult ill posed inverse problems. The corresponding electric potentials or magnetic fields measurements are far from sufficient to yield a unique solution, which can only be obtained by adding constraints or regularization. Most classical approaches (MNE, wMNE, LORETA, MCE, see [START_REF] Baillet | Encyclopedia of Computational Neuroscience, chapter Forward and Inverse Problems of MEG/EEG[END_REF] and references therein) do not explicitly use time dependence: source parameters are estimated at each time point and time correlations are not directly exploited (unless some pre-processing techniques used prior to inverse problem resolution). All those solutions can be interpreted as maximum a posteriori estimates, with gaussian noise distribution and various prior choices (often Gaussian or Laplacian).

An alternative approach was proposed in [START_REF] Lina | Waveletbased localization of oscillatory sources from magnetoencephalography data[END_REF], based on the Maximum Entropy on the Mean (MEM) principle combined with wavelet representation of time courses. This Bayesian technique introduces more freedom in the modeling, and yields smooth optimization problems of dimension much smaller than the number of sources. Nevertheless, while time correlations are to some extent captured by wavelet coefficients, inversion is still performed coefficientwise.

We describe here a space-time extension of the wMEM approach of [START_REF] Lina | Waveletbased localization of oscillatory sources from magnetoencephalography data[END_REF] in which time dependence is explicitely modeled and accounted for. To overcome the curse of dimensionality, we rely on three main ingredients: the Kronecker product factorization of covariance matrices (the noise covariance, in the spirit of [START_REF] Bijma | A mathematical approach to the temporal stationarity of background noise in MEG/EEG measurements[END_REF], and the source covariances), the matrix formulation of the MEM optimization problem, which generates considerable savings, and space and time dimension reduction techniques. The resulting (smooth and concave) optimization problem is solved using standard tools.

We apply this approach to the MEG inverse problem, on a dataset originating from a study of slow waves in deep sleep MEG signals [START_REF] Marc | Electromagnetic neuroimaging of sleep slow waves generators from MEG-EEG data fusion[END_REF]. Results on real data will be discussed in a forthcoming publication, we present here results of a simulation study that shows the ability of the approach to localize sources with such a given time course.

The paper is organized as follows. We describe in Sections 2 and 3 the main aspects of our approach, and provide numerical results in Section 4. Section 5 is devoted to conclusions and discussion.

PROBLEM STATEMENT, MEM APPROACH

Problem statement, notations

Observed signals are modeled as multivariate time series ℓ → z(ℓ) ∈ R J0 (J 0 being the number of sensors), and similarly the observation noise writes ℓ → n 0 (ℓ) ∈ R J0 . In the framework of the distributed sources model, we denote by K the number of mesh points on the cortex, and by ℓ → x(ℓ) ∈ R K the corresponding time courses. Denoting by G 0 ∈ R J0×K the lead-field matrix that summarizes the forward problem (propagation from the cortical surface to the sensors), this leads to the time domain observation equation

z(ℓ) = G 0 x(ℓ) + n 0 (ℓ) , ℓ = 1, . . . L 0 , (1)
L 0 being the number of time samples. These multivariate time series are reshaped as matrices Z, N 0 ∈ R L0×J0 , X ∈ 978-1-5386-5477-4/18/$31.00 c 2018 IEEE R L0×K . This yields the matrix formulation:

Z = XG T 0 + N 0 , (2)
where z = vec(Z T ). Here • T denotes matrix transposition and vec(•) vertical concatenation of columns. After dimension reduction (wavelet transform and projection onto a L dimensional wavelet subspace in time domain, and space domain reduction to J spatial filters, see Section 2.3 below), we obtain in matrix form

D = W G T + N = Y H T . (3)
Here,

D, N ∈ R L×J , W ∈ R L×K , and Y = [W, N ] (resp. H = [G; I J ]
) denotes the horizontal (resp. vertical) concatenation of matrices W and N (resp. G and I J ). This is the inverse problem that will be of interest to us here.

The vector MEM approach

The Maximum Entropy on the Mean (MEM, see [START_REF] Gamboa | Bayesian methods and maximum entropy for ill-posed inverse problems[END_REF][START_REF] Le Besnerais | A new look at entropy for solving linear inverse problems[END_REF] for presentations) is a generic Bayesian approach for solving linear inverse problem. It was first used for the MEG inverse problem in [START_REF] Lina | Waveletbased localization of oscillatory sources from magnetoencephalography data[END_REF], combined with wavelet decomposition, resulting in the so-called wMEM method from which this work is strongly inspired. Let d be a column of D and y the corresponding column of Y . MEM models the "source and noise" vector d as a random vector, and the observation equation d = Hy is replaced with E {Hy} = d, E {•} denoting expectation. MEM requires specifying a reference probability distribution on y. Noise and sources are assumed independent, and MEM seeks the probability distribution on y that minimizes the Kullback-Leibler divergence to the reference distribution, under the constraint E {Hy} = d. Sources are finally estimated as the expectation E {y} of y with respect to the so-obtained distribution. This constrained minimization problem turns out [START_REF] Gamboa | Bayesian methods and maximum entropy for ill-posed inverse problems[END_REF][START_REF] Le Besnerais | A new look at entropy for solving linear inverse problems[END_REF] to be equivalent to the maximization of a concave function D(λ) of an auxiliary variable λ, whose dimension equals the dimension of observations z,

D(λ) = λ • d -F * y (H T λ) , (4)
where λ • d is the inner product in the observation space, and the function F * y , called log-partition function, is fully specified by the reference probability distribution. Given the optimizer λ * = arg min λ D(λ), the estimate ŷ for the "sources and noise" vector finally reads ŷ = H∇D(λ * ) .

(5)

In [START_REF] Lina | Waveletbased localization of oscillatory sources from magnetoencephalography data[END_REF], this approach was proposed and tested for the MEG inverse problem. In this approach, the multi-sensor observations d are wavelet coefficients, leading to a J 0 -dimensional optimization problem (which is a major asset of the approach, given that J 0 ≪ K). The reference model is a Gaussian mixture model, based upon a parcellization of the cortex into independent regions. An important aspect is the parametrization of the reference model, for which the authors propose a generic strategy exploiting the cortex geometry and the source pre-localization MSP technique of [START_REF] Mattout | Multivariate source prelocalization (MSP): Use of functionally informed basis functions for better conditioning the MEG inverse problem[END_REF].

The extension to (vectorized) space-time data is straightforward, but increases significantly the dimension of the optimization problem, which in addition involves calculations in a very high dimensional space (cortex × time). The latter can be simplified by suitable prior choices (described in Section 3) and matrix formulation. However, dimension reduction in the (space-time) observation space is also necessary. Notice that the space-time extension also increases the number of parameters in the model.

Dimension reduction

Wavelet transform [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF] provides alternative representations for signals, that often have the property of concentrating the relevant information in a small number of coefficients, and/or enforcing decorrelation, yielding sparse or diagonal dominant covariance matrices. We perform a channel-wise orthonormal wavelet transform, followed by a selection of the relevant coefficients. Here the retained coefficients are selected a priori, and are the same for all trials (wavelet coefficients with largets trial averaged energy are selected, coefficients perturbed by boundary effects are not retained). The selection could also be done adaptively, we don't address this question here. Wavelet transform and coefficient selection result in an observation matrix W 1 ∈ R L×J0 , L being the number of retained coefficients, generally L < J 0 .

The sensor domain can also be reduced using a spatial filtering. Several approaches could be chosen (see e.g. [START_REF] Sekihara | Adaptive Spatial Filters for Electromagnetic Brain Imaging[END_REF] for a review), we limit ourselves here to simple PCA-based dimension reduction: projection onto the first J principal axes. Among possible extensions, discriminant filters such as introduced in [START_REF] Spinnato | Finding EEG spacetime-scale localized features using matrix-based penalized discriminant analysis[END_REF] are an interesting perspective.

After time and space dimension reduction, the observed data takes the form of a matrix D ∈ R L×J , and the inverse problem to be solved writes

D = W G T + N , (6)
where G and N are the projections of the lead-field matrix G 0 and the noise N 0 on the reduced sensor space. W ∈ R L×K represents the unknown source wavelet coefficients.

MATRIX WMEM FOR SPACE-TIME DATA

We now turn to the time-space model. The data to be processed is a time-space matrix D of size L × J. We denote corresponding vectorized matrices with boldface lowercase symbols: observations d = vec D T , noise n = vec N T , wavelet source coefficients w = vec W T and so on.

Matrix gaussian distribution, Kronecker covariance factorization

Our model is based upon the following assumptions and models on noise and cortical sources.

Noise model

The observation noise n is distributed following a matrix normal law [START_REF] Arjun | Matrix Variate Distributions[END_REF], i.e.

n = vec N T ∼ N (0, Σ t N ⊗ Σ s N ) (7)
with Kronecker covariance matrix (as proposed in [START_REF] Bijma | A mathematical approach to the temporal stationarity of background noise in MEG/EEG measurements[END_REF]) where Σ t N ∈ R L×L is the time covariance matrix, and Σ s N ∈ R J×J is the space (i.e. sensor) covariance matrix. Notice that this factorized model requires estimating and storing L 2 + J 2 numbers rather than L 2 J 2 . It is worth mentioning that Kronecker product covariance matrices can be estimated using dedicated algorithms (the so-called flip-flop algorithms), whose convergence has been studied and proven [START_REF] Dutilleul | The MLE algorithm for the matrix normal distribution[END_REF][START_REF] Muni | Models with a Kronecker product covariance structure: Estimation and testing[END_REF].

Source model

Recall that the source space is a mesh of the cortical surface, of size K. The source model is based upon a parcellization of the cortical surface, i.e. a segmentation into P connected regions, called parcels, of size K p , p = 1, . . . P .

The reference distribution on sources is then written as the product of P independent parcel laws. Following [START_REF] Lina | Waveletbased localization of oscillatory sources from magnetoencephalography data[END_REF], activity in each parcel w p ∈ R LKp , p = 1, . . . P is distributed according to a mixture of two matrix normal distributions, with respective means ω p and 0, and covariance matrices

Σ p = Σ t p ⊗ Σ s p (Σ t p ∈ R L×L , Σ s p ∈ R Kp×Kp ) and v p I LKp : w p ∼ α p N (ω p , Σ p ) + (1 -α p )N (0, v p I LKp ) . ( 8
)
State 1 is the "active state" (non-zero mean), and state 0 is the "silent state" (zero mean white noise), coefficient α p represents the probability for parcel p to be active. The parcel spatial covariance matrices Σ s p encode for correlations on the cortex. In the spirit of [START_REF] Lina | Waveletbased localization of oscillatory sources from magnetoencephalography data[END_REF], Σ s p is set to the restriction to parcel p of the covariance matrix exp(-ρ∆) of a diffusion process on the mesh graph, ∆ = D -A being the graph Laplacian (A and D are respectively the adjacency and degree matrices of the graph), and ρ some fixed parameter (unlike [START_REF] Lina | Waveletbased localization of oscillatory sources from magnetoencephalography data[END_REF] which used a truncated series expansion of the exponential, we prefer to stick here to the full matrix exponential, that yields better conditioned covariance matrices).

Kronecker-wMEM approach

We are now in position to describe our approach and solve the inverse problem associated with the model given in matrix form above. Using the parcellized Gaussian mixture model given above, the MEM principle leads to maximize the objective function λ ∈ R LJ → D(λ) defined in (4).

As in [START_REF] Lina | Waveletbased localization of oscillatory sources from magnetoencephalography data[END_REF], the independence of sources and noise, together with the independence of parcels, yield a splitting of the log partition function F * = P p=1 F * p + F * n . The assumptions on the noise yield closed form expressions for the objective function D(λ) as well as the parcel log-partition functions. However, these expressions are not easily amenable to numerical calculations as they involve algebraic manipulations in high dimensional space: each w p belongs to a high dimensional space of dimension or the order of 5000, corresponding matrix vector products are costly. Nevertheless, they can be conveniently re-formulated in matrix form. Defining Λ ∈ R L×J by λ = vec Λ T , the objective function writes

D(Λ) = Tr D T Λ - 1 2 Tr Λ T Σ t N ΛΣ s N - P p=1 F * p (ΛG p ) , (9)
where Tr denotes the matrix trace, and G p is the submatrix of G obtained by restricting to parcel p. The parcel log-partition functions are as follows:

F * p = ln α p exp F * p,1 +(1-α p ) exp F * p,0 , (10)
where the Gaussian log-partition functions are given by

F * p,0 (U ) = v p 2 Tr U T U (11)
F * p,1 (U ) = Tr U T Ω p + 1 2 Tr U T Σ t p U Σ s p , (12)
for all U ∈ R L×Kp , and Ω p is defined by ω p = vec Ω T p . Finally, denoting by Λ * the unique optimizer of D, the estimate W p ∈ R LKp for each parcel p is obtained in matrix form, and given as

W p = α * p Ω p + Σ t p Λ * G p Σ s p + (1 -α * p ) v p Λ * G p , (13)
where the updated activity probabilities read

αp = α p α p +(1-α p ) exp F * p,0 (Λ * G p )-F * p,1 (Λ * G p ) . (14)
From this the vector form is readily computed as

ŵp = vec W T p , (15)
and the time courses of the estimated sources are obtained by inverse wavelet transform.

NUMERICAL RESULTS

The algorithm was implemented in the Matlab R computing environment, using the WaveLab package [START_REF] Buckheit | About wavelab[END_REF] for wavelet transform. Numerical optimization of the objective function D(λ) was performed using the minFunc function [START_REF] Schmidt | minFunc: unconstrained differentiable multivariate optimization in Matlab[END_REF], that implements an adaptive step quasi-Newton (BFGS) algorithm. With the above data, computing time for an inversion is around 5 secs on a laptop (intel core i7-3687U CPU, 2.10GHz × 4, 16 Go RAM).

Data, and model specifications

The results presented here originate from a study of slow waves in deep sleep MEG data [START_REF] Marc | Electromagnetic neuroimaging of sleep slow waves generators from MEG-EEG data fusion[END_REF]. The original dataset (1 subject) consists in 180 trials, 4 seconds long, sampled at 50 Hz, recorded at 272 sensors. The trials were epoched and aligned by experts using EEG data. Besides, 109 "signal free" recordings were used for noise statistics estimation. Notice that noise contains both sensor noise, and non-interest background activity. Noise time and space covariance matrices were estimated using the flip-flop algorithm [START_REF] Dutilleul | The MLE algorithm for the matrix normal distribution[END_REF]. Time courses were wavelet transformed using Daubechies 6 orthonormal wavelet basis, after zero padding. Out of 256 wavelet coefficients, 62 most relevant coefficients were selected. Wavelet coefficients influenced by spurious boundary effects were not included. As for sensor domain dimension reduction, data were projected onto the 15 principal components with largest principal values, yielding a 62 × 15 data matrix for each trial. On original data 15 such virtual channels turn out to capture more than 98% of inertia, however since we use simulated data here the percentage is even larger (more than 99%).

The cortical surface (originating from MRI measurements) was discretized with 10002 mesh points, and parcellized into 156 connected parcels (of size ranging from 22 to 175 mesh points, the mean and median being around 60 mesh points), based on anatomical neighborhood criteria. The model parameters were set as follows, some being fixed, some being estimated from data using a simplified reference model (Gaussian instead of Gaussian mixture, equivalent to wMNE). Following [START_REF] Lina | Waveletbased localization of oscillatory sources from magnetoencephalography data[END_REF], parcel space covariance matrices Σ s p were constructed using the graph Laplacian of each parcel, with ρ = .3. Assuming time decorrelation (i.e. Σ t p = v p I L , with v p initialized to a constant value, used as signal to noise ratio estimate) and zero mean leads to a Gaussian reference model. Using that reference model, a preliminary estimate for sources was obtained. From the latter, active state parcel means ω p and parcel wavelet covariance matrices Σ t p were then estimated, to serve as parameters for the Gaussian mixture reference model. Parcel activity probabilities (which control the sparsity of the source estimates) were set to a constant value, α p = .25.

A simulation study

Running the proposed approach on the slow waves dataset turned out to produce a neat time course for the slow wave, which we used in the simulations described below. The time course ϕ(ℓ) is displayed in Fig. 1. Simulated datasets were created by 1) generating a connected cortical region Ω centered at a random seed, 2) generating a current distribution j 0 on the cortical surface, set to zero outside Ω and to j 0 (k, ℓ) = ϕ(ℓ) for k ∈ Ω, 3) propagating to the sensors using the leadfield matrix, and 4) adding a noise realization taken randomly in the above mentioned "signal free" trials, with a prescribed signal to noise ratio. We stress that these regions are independent from the parcels used in the inversion algorithm. From these simulations, performances were evaluated as follows. First a global space-time reconstruction index was computed as the normalized inner product of original cortical sources j 0 (k, ℓ) = ϕ(ℓ)1 Ω (k) and reconstructed ones j rec (k, ℓ):

ι = j 0 , j rec j 0 j rec ,
where inner product and norms are computed with respect to time and space. While the corresponding numerical values are hard to interpret quantitatively, they show a clear improvement when moving from the Gaussian reference to the Gaussian mixture one. The second criterion rests on decision theory ideas. The reconstructed sources being space-time data, detection of active/silent regions is difficult because activity is time dependent. We focus here on the ability of the method to detect activity with a given time profile, and report on detection performances based upon the following criterion: denoting by j rec (k, ℓ) the estimated source at time t and on the cortical mesh point k, and by ϕ(ℓ) the input time profile (the time course displayed in Fig. 1), we compute for each mesh point k = 1, . . . K an activity score, denoted by κ:

κ(k) = j rec (k, •), ϕ j rec (k, •) ϕ ,
where inner product and norms are computed with respect to time. This quantity, which ranges from 0 to 1, measures the ability of the method to recover sources with a given time course. For a fixed threshold τ , a cortical mesh point k is declared active if κ(k) ≥ τ . Varying τ within [0, 1] leads to Receiver Operating Characteristic (ROC) curves, and the area under curve (AUC) is used to assess the performances of the method in terms of detection. The closer AUC to 1 the better the performances, values near 80% being considered good.

We report here the quantitative evaluations of simulations, for three variants of the method. In the first one (denoted by G), the reference model is a Gaussian model with diagonal source covariance matrix, which provides a closed form solution (equivalent to a variant of the standard wMNE algorithm). The second one (GM) is based upon the Gaussian mixture model described above, whose parameters have been estimated from the Gaussian estimate. In the third one (uGM), the Gaussian mixture MEM algorithm has been run once again, using updated parameters estimated from the latter Gaussian mixture MEM estimate.

Results have been obtained on 100 different source trials (i.e. 100 different connected patches of active sources), averaged over 109 realizations of noise. The results of Table 1 provide the corresponding average ι-scores and κ-score based AUCs, which have also been averaged over source trials. The input SNR was set to SNR = 6dB (i.e. signal was 2 times larger than noise). The table also provides medians and standard deviations. As can be seen, moving from the Gaussian reference model (which does not account for time correlations in the source model) to Gaussian mixture priori significantly improves performances, both in terms of ιscore and κ-based AUCs, with an increase of the standard deviation. Refining the Gaussian mixture reference distribution third column) further sligthtly improves the average and median results, with a further increase in standard deviation. A closer look at results (not reproduced here) shows that this additional improvement is noticeable when the GM reference is already significantly better than the G rerefence. However, when the quality of reconstruction is not so good (which generally originates from poor quality parameters), updating parameters tends to degrade further the quality of reconstruction.

According to the discussion in [START_REF] Lina | Waveletbased localization of oscillatory sources from magnetoencephalography data[END_REF], AUC results can be biased by the imbalance of silent/active cortical mesh point. For that reason, we also provide corresponding results obtained by (randomly) selecting in each simulation a number of silent cortical mesh points equal to the number of active points. This leads to the so-called restricted AUC (AUC R in Table 1). Differences between AUC and AUC R are hardly noticeable.

Similar results have been obtained for higher values of input signal to noise ratio (namely, SNR = 14dB and SNR = 20dB). No striking difference with Table 1 is to be mentioned, except for the fact that higher SNR improves the ιindex, while AUC means and medians are not really affected.

Real data

While the present paper focuses on simulation results, the approach has also been tested on real sleep slow waves data. Results will be described and discussed in details in a forthcoming publication, we simply give a short example here. The complete dataset consists in 180 epochs (trials) recorded at J 0 = 272 sensors, which have been co-registered by an expert using additional EEG recordings, to be centered on a slow wave. Sampling rate was 50Hz, the duration of each epoch being 4secs. We report here on results obtained on the trial average.

CONCLUSION

We have presented in this paper first results on an extended wavelet MEM algorithm for time-space source localization from MEG measurements. The goal of this extension is to account explicitely for time correlations in the source space, which are not exploited (or exploited implicitely by a change of representation space, i.e. going to wavelet or time frequency space) in most classical approaches. The resulting curse of dimensionality is addressed using various dimension reduction tools. As a result, our numerical simulations (from realistic data) confirm that accounting for time correlations indeed improves precision in terms of time resolution, and that the sparsity properties induced by the gaussian mixture reference distribution also yields significant improvements in terms of detection performances.

A main difficulty of the approach lies in the choice of the model parameters. In this paper, some of these parameters were estimated from a first quick inversion, while some others were chosen by the user. Fully adaptive parameter choice would be desirable. Also, a natural follow up would be the study of spatio-temporal networks an the source level, for example exploiting space-time source covariance matrices. These questions will be addressed in a forthcoming publication, together with a more complete simulation study and extensive applications to real data.

Fig. 1 .

 1 Fig. 1. Time course of the simulated signals

  of 10002 grid points. Principal component analysis was performed on estimated sources (a matrix of size 201 × 10002), yielding a strong contribution of the first principal component (PC1, 56% of inertia, to be compared with 12% for the second PC). Corresponding time loadings are displayed in Fig.1, and form the time profile that was actually used in our simulations. Space (i.e. cortical) loadings are displayed in Fig.2, top left. The role of parcels appears clearly, as well as the spatial sparsity of the estimated sources (which is a by-product of the gaussian mixture model). Projections onto the cortical surface of PC1 show a strong localization in the frontal area, which is the expected localization for slow waves. PC1: source loadings (b) PC1 topography: front (c) PC1 topography: right (d) PC1 topography: left

Fig. 2 .

 2 Fig. 2. Spatial projection of the first principal component of estimated sources. a): loadings as a function of mesh points (ordered by increasing parcel number); b,c,d) projection of PC1 loadings onto the cortical surface (threshold set to 40% of maximal value).

Table 1 .

 1 The cortical surface was sampled, yielding a mesh Evaluation for the KMEM reconstruction method for three reference distributions: Gaussian, Gaussian mixture and updated Gaussian mixture. Input SNR was set to 6 dB.

	Criteria		G	GM uGM
		mean	.087 .207 .219
	ι-index median .077 .209 .233
		std-dev .037 .097 .121
		mean	.710 .833 .839
	AUC	median .692 .856 .891
		std-dev .098 .122 .159
		mean	.710 .833 .839
	AUC R	median .694 .857 .890
		std-dev .100 .123 .160

Acknowledgements

Part of this work was done while M.C. Roubaud and B. Torrésani were visiting the Centre de Recherches Mathématiques (CRM, UMI 3457) at Université de Montreal, both wish to thank CNRS for support and CRM for hospitality.


Jean-François Coulombel
email: jean-francois.coulombel@math.univ-toulouse.fr

Mark Williams
email: williams@email.unc.edu.

On the Mach stem configuration with shallow angle

The aim of this article is to explain why similar weak stability criteria appear in both the construction of steady Mach stem configurations bifurcating from a reference planar shock wave solution to the compressible Euler equations, as studied by Majda and Rosales [Stud. Appl. Math. 1984], and in the weakly nonlinear stability analysis of the same planar shock performed by the same authors [SIAM J. Appl. Math. 1983], when that shock is viewed as a solution to the evolutionary compressible Euler equations. By carefully studying the normal mode analysis of planar shocks in the evolutionary case, we show that for a uniquely defined tangential velocity with respect to the planar front, the temporal frequency which allows for the amplification of highly oscillating wave packets, when reflected on the shock front, vanishes. This specific tangential velocity is found to coincide with the expression given by Majda and Rosales [Stud. Appl. Math. 1984] when they determine the steady planar shocks that admit arbitrarily close steady Mach stem configurations. The links between the causality conditions for Mach stems and the so-called Lopatinskii determinant for shock waves are also clarified.

Introduction

The stability analysis of shock waves in gas dynamics now has a very long history, dating back to pioneering works for instance by D'yakov [D'y54], Erpenbeck [START_REF] Erpenbeck | Stability of step shocks[END_REF] and followers, see, e.g., [START_REF] Swan | Shock wave stability[END_REF][START_REF] Barmin | Stability of shock waves[END_REF] and references therein. At the linearized level, determining stability amounts to finding unstable and/or neutrally stable eigenmodes. The most favorable situation corresponds to nonexistence of unstable nor neutrally stable eigenmodes. Following Majda's memoirs [START_REF] Majda | The stability of multidimensional shock fronts[END_REF][START_REF] Majda | The existence of multidimensional shock fronts[END_REF], this regime will be referred to as that of uniform stability. When unstable eigenmodes (of positive real part) occur, violent instability is expected to take place. We shall mainly be concerned here with the intermediate situation where unstable eigenmodes do not occur but neutrally stable (that is, purely imaginary) eigenmodes do arise. This regime will be referred to as that of weak stability. From a more mathematical point of view, the regime we shall consider corresponds to the so-called WR (for Weak Real) class identified in [START_REF] Benzoni-Gavage | Generic types and transitions in hyperbolic initial-boundary-value problems[END_REF] but we shall only be concerned here with the particular problem of shock waves in gas dynamics.

In the first part [START_REF] Majda | A theory for spontaneous Mach stem formation in reacting shock fronts. I. The basic perturbation analysis[END_REF] of a series of papers, Majda and Rosales considered the regime of weak stability for reacting shocks and identified some weakly nonlinear waves that exhibited an amplification phenomenon. The weakly nonlinear waves considered in [START_REF] Majda | A theory for spontaneous Mach stem formation in reacting shock fronts. I. The basic perturbation analysis[END_REF] are approximate solutions to the evolutionary compressible Euler equations that are small, high frequency perturbations of a planar reference shock with zero tangential velocity. The analysis in [START_REF] Majda | A theory for spontaneous Mach stem formation in reacting shock fronts. I. The basic perturbation analysis[END_REF] has been recently generalized by the authors in [START_REF] Coulombel | The Mach stem equation and amplification in strongly nonlinear geometric optics[END_REF]. In the second part [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] of their series, Majda and Rosales considered the steady Euler equations in two space dimensions and proved that in the exact same regime of weak stability, and for some specific nonzero tangential velocity, a step shock could 'bifurcate' into a family of steady Mach stems with shallow angle (see Figures 1 and2 hereafter for an illustration), thus providing "a completely independent confirmation of that theory" (quote from [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF]). The links between the two problems in [START_REF] Majda | A theory for spontaneous Mach stem formation in reacting shock fronts. I. The basic perturbation analysis[END_REF][START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] and the appearance of the exact same weak stability condition are somehow hidden in lengthy calculations. It is the purpose of this article to clearly explain the role of the nonzero tangential velocity in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] and its connection with the so-called Lopatinskii determinant that has now been repeatedly computed for decades in the stability analysis of planar shock waves with zero tangential velocity. As should be clear from our analysis below, the problems studied in [START_REF] Majda | A theory for spontaneous Mach stem formation in reacting shock fronts. I. The basic perturbation analysis[END_REF] and [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] are not so independent as they might look at first glance.

Our main conclusions can be summarized as follows : the tangential velocity exhibited in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] for the analysis of step shocks bifurcating into steady Mach stems with shallow angle coincides with the root to the Lopatinskii determinant encoding the stability of the planar shock with same density/entropy/normal velocity but zero tangential velocity. This first result is proved in Proposition 2.2 below. Now a crucial property of the Euler equations is Galilean invariance which will reflect here into the result of Lemma 2.1 that connects the Lopatinskii determinant associated with two shock waves that share the same density/entropy/normal velocity but do not have the same tangential velocity. Combining the results of Proposition 2.2 and Lemma 2.1, we find that the tangential velocity exhibited in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] is equivalently determined by the requirement that the associated Lopatinskii determinant vanishes at the time frequency zero, hence the connection with the steady Euler equations. We review some of the arguments in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] and explain why the causality conditions used there as an admissibility criterion for Mach stems turn out to yield the same causality conditions used in [START_REF] Majda | A theory for spontaneous Mach stem formation in reacting shock fronts. I. The basic perturbation analysis[END_REF] to discriminate between incoming and outgoing wave packets.

The article is organized as follows. In Section 2, we briefly review the normal mode analysis which yields uniform/weak stability criteria for shock waves in evolutionary gas dynamics. In the weak stability regime, we verify that the expression of the tangential velocity given in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] coincides with the root to the Lopatinskii determinant. This first observation, based on 'brute force computations', is explained with further details in Section 3 where we prove that the bifurcation problem for the steady Euler equations considered in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] amounts to determining a tangential velocity for which the Lopatinskii determinant vanishes at the time frequency zero. We also give a complete construction of the family of Mach stems bifurcating from the reference planar shock, thus completing and clarifying some of the arguments in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF]. In particular, we make precise the assumptions on the pressure law under which the Mach stem construction can be achieved.

The normal mode analysis of shock waves

This Section could deal with any space dimension d ≥ 2, but since Section 3 will deal specifically with two-dimensional flows, we restrict from now on to d = 2 in order to keep the same notation throughout the whole article. We follow the presentation in [START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF] and consider a compressible inviscid fluid endowed with a complete equation of state e = e(τ, s). Here τ denotes the specific volume of the fluid, s denotes the specific entropy and e denotes the specific internal energy. The pressure p and temperature T are defined by the fundamental law of thermodynamics de = -p dτ + T ds .

The (evolutionary) compressible Euler equations in two space dimensions are written in the compact form (2.1)

∂ t f 0 (U ) + ∂ x 1 f 1 (U ) + ∂ x 2 f 2 (U ) = 0 ,
where U = (τ, u, s) is the four component vector1 of unknowns and u = (u, v) ∈ R2 is the fluid velocity.

The fluxes f α , α = 0, 1, 2, in (2.1) are given by

f 0 (U ) :=      ρ ρ u ρ v 1 2 ρ |u| 2 + ρ e      , f 1 (U ) :=       ρ u ρ u 2 + p ρ u v 1 2 ρ |u| 2 + ρ e + p u       , (2.2) f 2 (U ) :=       ρ v ρ u v ρ v 2 + p 1 2 ρ |u| 2 + ρ e + p v      
, where ρ := 1/τ denotes the density. Our assumptions on the equation of state are (part of) the classical Bethe-Weyl inequalities (we refer again to [START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF]):

p > 0 , T > 0 , ∂ 2 e ∂τ 2 > 0 , ∂ 2 e ∂s ∂τ < 0 , ∂ 3 e ∂ 3 τ < 0 .
We then define the sound speed c and the so-called Grüneisen coefficient Γ by setting

c 2 := τ 2 ∂ 2 e ∂τ 2 = -τ 2 ∂p ∂τ , Γ := - τ T ∂ 2 e ∂s ∂τ = τ T ∂p ∂s ,
both being positive quantities. A shock wave 2 is a piecewise constant solution to (2.1) satisfying Lax shock inequalities [START_REF] Lax | Hyperbolic systems of conservation laws[END_REF]. In the context of the Euler equations (2.1), one can always perform a Galilean change of frame and rotate the coordinate axes so that the shock wave is steady and reads as depicted in Figure 1. The Rankine-Hugoniot conditions, which ensure that (2.3) is a weak solution to (2.1), and Lax shock inequalities then read

(2.3) U = U 0 := (τ 0 , u 0 , s 0 ) if x 2 > 0 , U 1 := (τ 1 , u 1 , s 1 ) if x 2 < 0 , x 1 x 2 u 0 τ 0 , s 0 u 1 τ 1 , s 1
j := -ρ 0 v 0 = -ρ 1 v 1 > 0 , u 0 = u 1 =: u , (2.4a) j 2 (τ 1 -τ 0 ) = p 0 -p 1 , e 1 -e 0 + p 1 + p 0 2 (τ 1 -τ 0 ) = 0 , (2.4b) 0 < -v 1 c 1 < 1 < -v 0 c 0 . (2.4c)
The tangential velocity u could also be set to zero by a Galilean change of frame, but we do not do so here in order to highlight the links with the analysis of Section 3. Observe that up to changing x 1 into -x 1 and the tangential velocity accordingly, we can always assume u ≤ 0 without loss of generality. This is the convention in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] and we follow it here.

The (linear) stability properties of the particular solution (2.3) have been made precise after a long series of contributions which we have partly recalled in the introduction. The analysis is based on a normal mode decomposition that is briefly summarized below. We refer to the appendix of [START_REF] Zumbrun | Stability of large-amplitude shock waves of compressible Navier-Stokes equations[END_REF] and to [BGS07, chapter 15] for a detailed and complete analysis of this stability problem, and just introduce the notation that will be useful later on for our purpose.

We introduce the Fourier variable η ∈ R dual to x 1 , and the Laplace variable z = δ -i γ, γ ≥ 0, dual to t. We focus on the state '1' behind the shock since there are no stable modes ahead of the shock. The eigenmodes exp(i ω x 2 ) under consideration correspond to complex numbers ω of nonpositive imaginary part3 such that there exists a nontrivial solution ( τ , u, v, ṡ) ∈ C 4 to the linear system

                   (z + u η + v 1 ω) τ -τ 1 η u -τ 1 ω v = 0 , (z + u η + v 1 ω) u - c 2 1 τ 1 η τ + Γ 1 T 1 η ṡ = 0 , (z + u η + v 1 ω) v - c 2 1 τ 1 ω τ + Γ 1 T 1 ω ṡ = 0 , (z + u η + v 1 ω) ṡ = 0 .
All quantities with a '1' index are evaluated behind the shock, that is, in the region {x 2 < 0} for (2.3). The analysis of the latter linear system gives rise to two eigenmodes ω 0 and ω -, which are defined by

(2.5) z + u η + v 1 ω 0 = 0 , (z + u η + v 1 ω -) 2 = c 2 1 (η 2 + ω 2 -) ,
where the choice of ω -in (2.5) is such that ω -has negative imaginary part when z also has negative imaginary part, and ω -is extended continuously up to real values of z, see [START_REF] Benzoni-Gavage | Multidimensional hyperbolic partial differential equations[END_REF]chapters 14 & 15]. The corresponding eigenspaces associated with ω 0 and ω -are

(2.6) E 0 (z, η) := Span            0 ω 0 -η 0     ,     Γ 1 T 1 τ 1 0 0 c 2 1            , E -(z, η) := Span     τ 1 (z + u η + v 1 ω -) c 2 1 η c 2 1 ω - 0     .
Specifying to the case z = 0, η = 1, we obtain

(2.7)

E 0 (0, 1) = Span            0 u v 1 0     ,     Γ 1 T 1 τ 1 0 0 c 2 1            , E -(0, 1) = Span     τ 1 (u + v 1 ω -(0, 1)) c 2 1 c 2 1 ω -(0, 1) 0     , with (2.8) ω -(0, 1) = 1 c 2 1 -v 2 1 v 1 u + c 1 sgn (u) u 2 + v 2 1 -c 2 1 .
The expression (2.8) is valid as long as the tangential velocity u satisfies u 2 +v 2 1 > c 2 1 (here sgn denotes the sign function). In that case, (0, 1) is a hyperbolic frequency for the linearized Euler equations in {x 2 < 0} because all the roots to (2.5) are real and they depend smoothly on (z, η) near (0, 1). The reason for choosing sgn (u) in (2.8) rather than the opposite sign is forced by the fact that ω -should be continuous with respect to (z, η) on the set {Im z ≤ 0 , η ∈ R}. The determination of the appropriate sign in (2.8) can also be interpreted as a causality condition, meaning that oscillating wave packets associated with the phase 0

• t + 1 • x 1 + ω -(0, 1) x 2
should have a group velocity that points inside the half space {x 2 < 0} (in other words, the second component of the group velocity should be negative). This selection criterion is used when computing the weakly nonlinear expansions in [START_REF] Majda | A theory for spontaneous Mach stem formation in reacting shock fronts. I. The basic perturbation analysis[END_REF].

The subspace E s (z, η) of values ( τ , u, ṡ)| x 2 =0 -under consideration is then the direct sum of E 0 (z, η) and E -(z, η), except when the eigenmodes ω 0 and ω -coincide. We shall mainly be concerned here with the frequency (z, η) = (0, 1) for which ω 0 and ω -do not coincide, and therefore refer to the above mentioned works for the precise decomposition of E s (z, η) whenever ω 0 = ω -. We thus consider

E s (z, η) := E 0 (z, η) ⊕ E -(z, η) ,
and refer from now on to E s (z, η) as the stable subspace. (It is a three dimensional subspace of C4 .)

The stability analysis of the step shock (2.3) amounts to determining whether there exist some frequencies (z, η) for which one can find a nonzero pair (χ, U ) ∈ C × E s (z, η) that satisfies the linearized Rankine-Hugoniot conditions which, following [BGS07, pages 426-427] with our notation, read 4

(2.9)

                       - v 1 τ 2 1 τ + 1 τ 1 v = -(ρ 1 -ρ 0 ) (z + u η) χ , v 1 τ 1 u = -(p 1 -p 0 ) η χ , - v 2 1 + c 2 1 τ 2 1 τ + 2 v 1 τ 1 v + Γ 1 T 1 τ 1 ṡ = 0 , - 1 2 c 2 1 τ 2 1 (τ 1 -τ 0 ) + p 1 -p 0 τ + T 1 1 + Γ 1 (τ 1 -τ 0 ) 2 τ 1 ṡ = 0 .
For future use (we refer once again to [BGS07, chapter 15]), it should be kept in mind that (2.9) is an equivalent formulation of the relation

(2.10) df 2 (U 1 ) U = -χ z f 0 (U 1 ) -f 0 (U 0 ) + η f 1 (U 1 ) -f 1 (U 0 ) ,
after some elementary manipulations on the rows of (2.10). The Jacobian matrix df 2 (U 1 ) is invertible because of Lax shock inequalities (2.4c) (its eigenvalues are v 1 -c 1 , v 1 and v 1 +c 1 ), so any nonzero solution (χ, U ) ∈ C × E s (z, η) to (2.9) must satisfy χ = 0. In other words, either (2.9) has no nonzero solution in C × E s (z, η), or the set of solutions is a one-dimensional subspace of C × E s (z, η). Equivalently, the stability analysis of the step shock (2.3) amounts to determining whether there holds

(2.11) df 2 (U 1 ) -1 z f 0 (U 1 ) -f 0 (U 0 ) + η f 1 (U 1 ) -f 1 (U 0 ) ∈ E s (z, η) .
The shock wave linear stability problem may be encoded as the determination of the roots of an appropriate determinant, which is usually referred to as the Lopatinskii determinant. More precisely, one can define a complex number ∆(u, z, η) such that ∆(u, z, η) = 0 if and only if there exists a nonzero (χ, U ) ∈ C × E s (z, η) solution to (2.9). We have highlighted here the dependence of the determinant ∆ on the tangential velocity u, in order to state the following result, which is very simple though fundamental if one wants to understand the links between the calculations performed in [START_REF] Majda | A theory for spontaneous Mach stem formation in reacting shock fronts. I. The basic perturbation analysis[END_REF] and [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF].

Lemma 2.1. Let the Lopatinskii determinant ∆(u, z, η) be defined5 such that ∆(u, z, η) = 0 if and only if there exists a nonzero (χ, U ) ∈ C × E s (z, η) solution to (2.9). Then there holds ∆(u, z, η) = ∆(0, z + u η, η) , for all complex number z of nonnegative imaginary part and all real number η.

The proof of Lemma 2.1 is elementary and is based on inspection of (2.5), (2.6) and (2.9). Namely, it follows from the above definitions that the change of parameters z := z + u η in (2.5), (2.6) and (2.9) does not affect the sign of the imaginary part of z and it reduces the above analysis to the case u = 0. We leave this elementary verification to the reader.

Unsurprisingly, the existence/nonexistence of a zero to ∆, and accordingly the sign of the imaginary part of z to discriminate between violent instability and neutral stability, is independent of u (a consequence of Galilean invariance). However, the precise location of the zeroes to ∆, provided that they exist, does depend on u.

The stability analysis of step shocks with zero tangential velocity can be summarized as follows (see [START_REF] Majda | The stability of multidimensional shock fronts[END_REF] and [BGS07, chapter 15]), where from now on M 1 ∈ (0, 1) denotes the Mach number -v 1 /c 1 behind the shock:

• If M 2 1 (τ 0 /τ 1 -1) < 1/(1 + Γ 1 )
, then ∆(0, z, η) = 0 for any nonzero pair (z, η) ∈ C × R with Im z ≤ 0. This regime corresponds to uniform stability. In particular, the function ∆(0, •, 1) does not vanish on the real line.

• If M 2 1 (τ 0 /τ 1 -1) > (1 + M 1 )/Γ 1 , then ∆(0, z, η) = 0
for some pair (z, η) with Im z < 0 and η = 0. This regime corresponds to violent instability. Moreover, the function ∆(0, •, 1) does not vanish on the real line.

• If (2.12) 1 1 + Γ 1 < M 2 1 τ 0 τ 1 -1 < 1 + M 1 Γ 1 ,
then there exists a uniquely determined velocity6 V > 0 such that ∆(0, z, η) = 0 if and only if7 η = 0 and z = ±V η. This regime corresponds to neutral stability and falls into the so-called WR class of [START_REF] Benzoni-Gavage | Generic types and transitions in hyperbolic initial-boundary-value problems[END_REF].

In the latter weakly stable case, the velocity V > 0 can be characterized as follows, see [BGS07, Theorem 15.1]:

(i) V 2 > c 2 1 -v 2 1 , (ii) 1 + M 2 1 -M 2 1 Γ 1 τ 0 τ 1 -1 V 2 < v 2 1 1 -M 2 1 τ 0 τ 1 , (2.13) (iii) k -1 + M 2 1 V 2 -v 2 1 1 -M 2 1 τ 0 τ 1 2 = k 2 V 2 M 2 1 V 2 -v 2 1 1 -M 2 1 ,
where the parameter k in (2.13)-(iii) is defined by:

k := 2 -M 2 1 Γ 1 τ 0 τ 1 -1 .
In the limit case

(2.14) 1 1 + Γ 1 = M 2 1 τ 0 τ 1 -1 , V tends to the 'glancing' value c 2 1 -v 2 1 . The other limit case M 2 1 τ 0 τ 1 -1 = 1 + M 1 Γ 1 ,
corresponds to a transition from weak stability to violent instability for which the Lopatinskii determinant vanishes at the frequency (z, η) = (1, 0). This case is studied in [START_REF] Serre | La transition vers l'instabilité pour les ondes de choc multi-dimensionnelles[END_REF].

For what follows, the key point to keep in mind is that if the function ∆(0, •, 1) vanishes on the real line, then either the shock wave (2.3) satisfies (2.12) or it satisfies one of two the limit cases of (2.12), namely (2.14). In the other limit case of (2.12), the Lopatinskii determinant vanishes for some real number z but with a zero tangential frequency η.

Our main result in this Section asserts that the tangential velocity exhibited in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] for giving rise to a bifurcation from a rectilinear shock to a family of steady Mach stems coincides (up to the sign convention) with the root V to the Lopatinskii determinant in the regime (2.12). This is a preliminary connection between the shock wave stability analysis, that deals with an evolutionary problem, and the steady Mach stem bifurcation problem explored in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF].

Proposition 2.2. Assume that the inequalities (2.12) are satisfied. Then the tangential velocity c ⋆ defined in [MR84, page 124] coincides with the velocity V for which the Lopatinskii determinant ∆(0, •, •) vanishes at (±V, 1). Thanks to Lemma 2.1, the tangential velocity u := -c ⋆ in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] thus satisfies ∆(u, 0, 1) = 0.

Proof of Proposition 2.2. We check the items (i), (ii) and (iii) in (2.13) one by one. For the sake of completeness, let us first recall how the tangential velocity c ⋆ is defined in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF]page 124]. Assuming that the inequalities (2.12) are satisfied, one first computes the (unique real) root Φ ∈ (M 1 , 1) to the second order polynomial equation

(2.15) Φ 2 + M 1 Γ 1 Φ -1 -Γ 1 + 1 -M 2 1 (τ 0 /τ 1 -1) M 2 1 = 0 .
Then one defines β ∈ (-1, M 1 ) and c ⋆ > 0 by the relations:

(2.16)

β := 2 M 1 -(1 + M 2 1 ) Φ 1 + M 2 1 -2 M 1 Φ , c ⋆ := c 1 1 -M 1 β 1 -β 2 .
We are therefore going to show that the velocity c ⋆ given in (2.16) satisfies the characterization in (2.13).

For the sake of simplicity, we introduce the notation (2.17)

ν := τ 0 τ 1 -1 > 0 ,
which measures the compression ratio of the shock, the inequality ν > 0 following here from (2.12) (anyway under rather mild additional assumptions on the pressure law, shock waves are always compressive discontinuities, see [START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF]).

• Verifying (2.13)-(i): one uses the definition (2.16) of β and computes

c 2 ⋆ = c 2 1 (1 -M 1 β) 2 1 -β 2 = c 2 1 (1 -M 1 Φ) 2 1 -Φ 2 > c 2 1 1 -M 2 1 = c 2 1 -v 2 1 ,
where the final inequality follows from Φ ∈ (M 1 , 1).

• Verifying (2.13)-(ii): for future use, we introduce the rescaled parameter:

(2.18)

y := M 1 Φ ∈ (M 2 1 , M 1 ) .
Since Φ satisfies (2.15), y in (2.18) is a root (and it is even the largest one) to the polynomial Q that is defined by:

(2.19) Q(Y ) := Y 2 + M 2 1 Γ 1 Y -M 2 1 (1 + Γ 1 ) + 1 -M 2 1 ν .
Using the relation

(2.20) c 2 ⋆ = c 2 1 (1 -M 1 Φ) 2 1 -Φ 2 = v 2 1 (1 -y) 2 M 2 1 -y 2 ,
which we have found at the previous step, and using also the above definition (2.17) of the parameter ν, verifying (2.13)-(ii) amounts to proving the inequality

1 + M 2 1 -M 2 1 Γ 1 ν (1 -y) 2 -1 -M 2 1 (1 + ν) M 2 1 -y 2 < 0 .
We now use the fact that y is a root to the polynomial Q in (2.19), and write equivalently the latter inequality as (2.21)

y 2 -M 2 1 Γ 1 ν 1 + M 2 1 + M 2 1 Γ 1 > 1 ν -2 1 -M 2 1 + ν M 2 1 4 + 3 -M 2 1 Γ 1 -ν 2 M 4 1 Γ 1 (2 + Γ 1 ) = (2-M 2 1 Γ 1 ν) (M 2 1 -1+ν M 2 1 (2+Γ 1 ))
.

From the inequalities (2.12), we know that 2 -M 2 1 Γ 1 ν is positive, and therefore, proving (2.21) amounts to showing

(2.22) y > M 2 1 -1 + ν M 2 1 (2 + Γ 1 ) ν 1 + M 2 1 + M 2 1 Γ 1 ,
where, recalling the definition (2.18), we know that y is the positive and therefore largest root to the second degree polynomial Q in (2.19). The conclusion then follows from the relation

Q M 2 1 -1 + ν M 2 1 (2 + Γ 1 ) ν (1 + M 2 1 + M 2 1 Γ 1 ) = 1 -M 2 1 2 (1 + ν) ν 2 (1 + M 2 1 + M 2 1 Γ 1 ) 2 1 -ν M 2 1 (1 + Γ 1 ) < 0 ,
where the final inequality is a consequence of (2.12). This proves that the quantity on the right hand side in (2.22) lies between the two roots of Q, so (2.22) is satisfied. Rewinding the arguments above, we have verified (2.13)-(ii).

• Verifying (2.13)-(iii): factorizing 1 -M 2 1 , we first observe that the velocity V satisfies (2.13)-(iii) if and only if

(2.23) (k -1) 2 -M 2 1 V 4 + (k -1) 2 + 1 -2 M 2 1 -2 ν k -1 + M 2 1 v 2 1 V 2 + v 4 1 1 -M 2 1 (1 + ν) 2 = 0 ,
with the parameter k defined as in (2.13)-(iii). We are now going to verify that the velocity c ⋆ in (2.16) satisfies (2.23). We use the relation (2.20) and substitute this expression of c 2 ⋆ (recall that y is defined in (2.18) and is a root of the polynomial Q defined in (2.19)). Multiplying by the obviously nonzero quantity v -4 1 (M 2 1 -y 2 ) 2 , we are reduced to showing the relation

(2.24) 1 -ν M 2 1 Γ 1 2 -M 2 1 (1 -y) 4 + 1 -M 2 1 (1 + ν) 2 M 2 1 -y 2 2 + M 2 1 -y 2 (1 -y) 2 1 -ν M 2 1 Γ 1 2 + 1 -2 M 2 1 -2 ν 1 + M 2 1 -ν M 2 1 Γ 1 = 0 .
The left hand side of (2.24) can be factorized as

ν Q(y) a 2 y 2 + a 1 y + a 0 ,
where Q is defined in (2.19), and

a 2 := 4 + ν 1 -M 2 1 -2 M 2 1 Γ 1 , a 1 := -4 1 + M 2 1 + ν M 2 1 Γ 1 3 + M 2 1 , a 0 := 1 + M 2 1 2 -ν M 2 1 1 -M 2 1 -ν M 2 1 Γ 1 1 + M 2 1 . Since Q(y)
is zero, the left hand side of (2.24) is indeed zero, as expected. This shows that the velocity c ⋆ in (2.16) satisfies all three conditions in (2.13) and therefore coincides with V . Proposition 2.2 may look unsatisfactory since the links between the steady Mach stem bifurcation problem and the evolutionary stability analysis encoded in the function ∆ are still hidden. In the following Section, we explain in a more explicit way why determining the velocity c ⋆ in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] corresponds to finding a nonzero solution (χ, U ) ∈ C × E s (0, 1) to the linear system (2.9) in which the parameters z, η, u are respectively set equal to 0, 1, -c ⋆ . The existence of such a nontrivial solution is possible if and only if the Lopatinskii determinant vanishes, which gives a more satisfactory justification for Proposition 2.2 than our preliminary 'brute force' calculations. We also clarify below the links between the quantities in (2.16) and various expressions that arise in the shock wave stability problem.

Mach stems bifurcating from a steady planar shock

In all what follows, we keep the notation U for the four component vector (τ, u, s) and we keep the notation u = (u, v) for the velocity. We review the Mach stem bifurcation problem considered in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] and explain its link with the shock wave stability problem (and more precisely its formulation (2.11)).

The problem considered in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] is the following: under which conditions on the shock (2.3) is it possible to construct a family (U 0 , U 1 , U 2 , U 3 , Θ, Φ, Ψ)(ε), with ε ∈ [0, ε 0 ] for some ε 0 > 0, such that, for all ε ∈ (0, ε 0 ], the wave pattern depicted in Figure 2 defines a steady Mach stem, and

8 lim ε→0 U 0 (ε) = U 0 , lim ε→0 U 1 (ε) = lim ε→0 U 2 (ε) = lim ε→0 U 3 (ε) = U 1 , (3.1) lim ε→0 Θ(ε) = π , lim ε→0 Φ(ε) = Φ 0 ∈ (π, 3 π/2) , lim ε→0 Ψ(ε) = Ψ 0 ∈ (Φ 0 , 2 π) .
It is also required that the family depends smoothly, say at least in a C 1 way, on ε. By a Mach stem, we mean that Figure 2 represents a four wave interaction at the origin where, for all ε ∈ (0, ε 0 ]:

• S 1 (ε) is a steady shock front where U 0 (ε) is the state ahead of the shock and U 2 (ε) the state behind,

• S 2 (ε) is a steady shock front where U 0 (ε) is the state ahead of the shock and U 1 (ε) the state behind,

• S 3 (ε) is a steady shock front where U 1 (ε) is the state ahead of the shock and U 3 (ε) the state behind,

• CD(ε) is a steady contact discontinuity.

In that case, the wave pattern depicted in Figure 2 yields a weak solution to the steady Euler equations.

The following causality conditions are also imposed for all ε ∈ (0, ε 0 ], see [MR84, page 123]:

• The tangential velocity of the fluid along S 3 (ε) points away from the origin.

• The velocity on either side of CD(ε) points away from the origin.

x 1 S 2 (ε)

x 2 S 1 (ε) U 0 (ε) U 1 (ε) U 2 (ε) U 3 (ε) S 3 (ε) CD(ε) Θ(ε) Φ(ε) Ψ(ε)
Figure 2: The Mach stem configuration with a shallow angle.

By the so-called triple shock Theorem [START_REF] Henderson | Triple-shock entropy theorem and its consequences[END_REF], S 1 (ε) should be the shock with the largest amplitude while S 2 (ε) will have a slightly smaller amplitude. In other words, we should have:

p(U 0 (ε)) < p(U 1 (ε)) < p(U 2 (ε)) = p(U 3 (ε)) .
There is of course much freedom in the parametrization of the family of Mach stems (and we shall see below that even with the choice of one parametrization, there can be more than one family of Mach stems that satisfy the asymptotic behavior (3.1)). As will follow from the proof of Theorem 3.1 below, it turns out that a convenient way to parametrize the family of Mach stems amounts to choosing π -Θ as the 'bifurcation parameter'. In other words, we ask from now on the angle Θ to be given by

(3.2) Θ(ε) := π -ε .
Why ε should be positive, and therefore Θ in (3.2) less than π, will be made clear in the proof of Theorem 3.1 below (see also [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] for similar considerations).

The main conclusion in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] is that such a 'shock to Mach stem' bifurcation occurs if and only if the inequalities (2.12) are satisfied and the tangential velocity u in (2.4a) equals -c ⋆ , with c ⋆ defined by (2.15), (2.16). The proof of the 'if' part in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] is skipped. Let us observe that because of our result in Proposition 2.2, the requirement on u in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] can be rewritten as u = -V where V > 0 is characterized by (2.13). Our main conclusion is of course in agreement with [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] and is summarized as follows (we recall the convention u ≤ 0). Theorem 3.1. Assume that there exists a family of Mach stems depending smoothly on ε ∈ [0, ε 0 ] where ε is given by (3.2), and having the asymptotic behavior (3.1) where the steady shock (2.3) satisfies (2.4). Then there exists a solution U ∈ E s (0, 1) to the linear system (2.9) in which one specifies z = 0, η = 1 and χ = 1. In particular, there holds:

1 1 + Γ 1 ≤ M 2 1 τ 0 τ 1 -1 < 1 + M 1 Γ 1 ,
and u = -V , where the velocity V is characterized by (2.13). Equivalently, the step shock (2.3), (2.4) with tangential velocity u can bifurcate into a Mach stem with shallow angle only if there holds ∆(u, 0, 1) = 0.

Conversely, if the steady shock (2.3) satisfies (2.4) together with the strict inequalities (2.12), and if furthermore u = -V where the velocity V is characterized by (2.13), then there exists a one parameter family of Mach stems satisfying (3.1), (3.2), and the state U 0 (ε) ahead of the Mach stem can be chosen to have the particular form

(3.3) ∀ ε ∈ [0, ε 0 ] , U 0 (ε) = (τ 0 , u(ε), v 0 , s 0 ) .
Our proof slightly differs from [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] and intends to clarify the link between the Mach stem bifurcation analysis and the relation (2.11) which deals with the normal mode analysis for the shock wave stability problem. Let us recall for future use that the validity of (2.11) is equivalent to the existence of a nonzero pair (χ, U ) ∈ C × E s (z, η) satisfying (2.9), which is also equivalent to ∆(u, z, η) = 0. In the case where z is real, this corresponds to a weak stability property for the shock wave (2.3). In the proof of Theorem 3.1 below, we clarify some of the computations and arguments in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] and, above all, we provide a complete proof of the 'if' part in Theorem 3.1 without any additional assumption on the pressure law than those made in Section 2. The assumptions on the pressure law were not complete in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF].

Proof of Theorem 3.1. We first prove the 'necessity' part of Theorem 3.1 and assume that a smooth family of Mach stems with shallow angle is given. We first introduce some notation and write9 :

U 0 (ε) = U 0 + ε U ′ 0 + o(ε) , U j (ε) = U 1 + ε U ′ j + o(ε) , j = 1, 2, 3 ,
where U 0 , U 1 are the two states of the step shock in (2.3). Following [MR84, page 130], we write the stationary Rankine-Hugoniot conditions for each of the four discontinuities S 1,2,3 (ε), CD(ε) and expand at the first order in ε. This yields

df 2 (U 0 ) U ′ 0 -df 2 (U 1 ) U ′ 1 = 0 , f 1 (U 0 ) -f 1 (U 1 ) + df 2 (U 0 ) U ′ 0 -df 2 (U 1 ) U ′ 2 = 0 , (3.4) -sin Ψ 0 df 1 (U 1 ) + cos Ψ 0 df 2 (U 1 ) (U ′ 1 -U ′ 3 ) = 0 , (3.5) -sin Φ 0 df 1 (U 1 ) + cos Φ 0 df 2 (U 1 ) (U ′ 3 -U ′ 2 ) = 0 . (3.6)
We can also pass to the limit in the zero normal velocity constraint for the contact discontinuity CD(ε) and in Lax shock inequalities for S 3 (ε), which yields (recall u 1 = (u, v 1 )):

(3.7) -u sin Φ 0 + v 1 cos Φ 0 = 0 , -u sin Ψ 0 + v 1 cos Ψ 0 = -c 1 .
Our goal now is to analyze the relations (3.4), (3.5), (3.6) and (3.7). Observe in particular that we shall never use the precise expression of U ′ 0 , meaning that the assumption made in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] that U ′ 0 only has one nonzero component is useless.

We can first equivalently rewrite (3.4) as

df 2 (U 0 ) U ′ 0 -df 2 (U 1 ) U ′ 1 = 0 , (3.8a) U ′ 1 -U ′ 2 = df 2 (U 1 ) -1 f 1 (U 1 ) -f 1 (U 0 ) . (3.8b)
Equation (3.8a) will ultimately determine U ′ 0 once we know U ′ 1 . Observe now that U ′ 0 does not appear in (3.8b), (3.5) and (3.6) so we may forget about U ′ 0 from now on and focus on the vectors U ′ j , j = 1, 2, 3. We are going to show below that the relations (3.5), (3.6) imply U ′ 1 -U ′ 2 ∈ E s (0, 1) where the stable subspace E s (0, 1) has been defined in Section 2, see (2.7). Therefore (3.8b) will give df 2 (U 1 ) -1 f 1 (U 1 ) -f 1 (U 0 ) ∈ E s (0, 1) , which is exactly (2.11) with z = 0 and η = 1. In other words, we shall have ∆(u, 0, 1) = 0, which equivalently means that (2.12) holds and u = -V with our convention for the sign of u. Let us make all these arguments precise.

We are first going to make the relations (3.5), (3.6) more explicit. Recalling that U denotes the vector (τ, u, s), we compute

df 1 (U ) = P (U ) -1       u -τ 0 0 -c 2 τ u 0 Γ T 0 0 u 0 0 0 0 u       , df 2 (U ) = P (U ) -1       v 0 -τ 0 0 v 0 0 -c 2 τ 0 v Γ T 0 0 0 v       , with (3.9) P (U ) :=       -τ 2 0 0 0 -τ u τ 0 0 -τ v 0 τ 0 τ 2 T ρ |u| 2 2 -e -p -τ u T -τ v T τ T       .
Multiplying (3.6) by P (U 1 ) on the left and using the first equation in (3.7), (3.6) is seen to be equivalent to

U ′ 3 -U ′ 2 ∈ Span            0 cos Φ 0 sin Φ 0 0     ,     Γ 1 T 1 τ 1 0 0 c 2 1            , with cos Φ 0 sin Φ 0 = u v 1 .
Recalling the decomposition (2.7) of the vector space E 0 (0, 1), it thus turns out that (3.6) equivalently reads U ′ 3 -U ′ 2 ∈ E 0 (0, 1). Let us now turn to (3.5). Multiplying (3.5) by P (U 1 ) on the left and using the second equation in (3.7), (3.5) is seen to be equivalent to

(3.10) U ′ 1 -U ′ 3 ∈ Span     τ 1 c 1 sin Ψ 0 -c 1 cos Ψ 0 0     ,
where the angle Ψ 0 satisfies the second equation in (3.7). We wish to show that (3.10) equivalently reads

U ′ 1 -U ′ 3 ∈ E -(0, 1)
where the vector space E -(0, 1) is given in (2.7). We thus need to show that the two vectors

    τ 1 (u + v 1 ω -(0, 1)) c 2 1 c 2 1 ω -(0, 1) 0     ,     τ 1 c 1 sin Ψ 0 -c 1 cos Ψ 0 0     ,
are collinear, or in other words, we need to show the relations

(3.11) sin Ψ 0 = c 1 u + v 1 ω -(0, 1) , cos Ψ 0 = - c 1 ω -(0, 1) u + v 1 ω -(0, 1)
.

This requires further knowledge on the angle Ψ 0 than just (3.7). We are therefore going to determine the angle Ψ 0 completely. By the causality conditions for the family of Mach stems, the velocity along the contact discontinuity CD(ε) must point away from the origin. Passing to the limit ε → 0, this means that the velocity u 1 is a positive multiple of the unit vector (cos Φ 0 , sin Φ 0 ). In other words, we have

cos Φ 0 = u u 2 + v 2 1 , sin Φ 0 = v 1 u 2 + v 2 1 .
Let us now recall that the angle Ψ 0 should satisfy Ψ 0 ∈ (Φ 0 , 2 π) so we must have cos Ψ 0 ∈ (cos Φ 0 , 1) , and sin Ψ 0 < 0 .

Let us now go back to (3.7). Since sin Ψ 0 is negative, cos Ψ 0 must satisfy (3.12)

-(c 1 + v 1 cos Ψ 0 ) = u 1 -cos 2 Ψ 0 < 0 .
Moreover, the existence of Ψ 0 implies, by the Cauchy-Schwarz inequality, the condition10

(3.13) u 2 ≥ c 2 1 -v 2 1 .
Recalling now the expression (2.8) of ω -(0, 1) and the expressions (3.14) for the sine and cosine of Ψ 0 , it is a simple exercise to verify that the relations (3.11) hold. Going then back to (3.10), we have obtained U ′ 1 -U ′ 3 ∈ E -(0, 1), and combining with U ′ 3 -U ′ 2 ∈ E 0 (0, 1), we have shown

U ′ 1 -U ′ 2 ∈ E s (0, 1)
. Summarizing what we have done so far, we have shown that the existence of a family of Mach stems satisfying (3.1), (3.2) and the causality conditions implies that the tangential velocity u satisfies (3.13) and that there exists a vector U ′ 1 -U ′ 2 ∈ E s (0, 1) such that (3.8b) holds. This means that the relation (2.11) holds for the frequencies (z, η) = (0, 1) hence the step shock (2.3) satisfies ∆(u, 0, 1) = 0. Using Lemma 2.1, we thus have ∆(0, u, 1) = 0 and we shall now use the conclusions of the stability analysis for shock waves with zero tangential velocity.

There are two cases. Either (3.13) holds with an equality sign, and in that case, one gets11 :

u = -c 2 1 -v 2 1 , ∆(0, u, 1) = 0 .
This situation occurs only if (we refer again to [BGS07, chapter 15])

1 1 + Γ 1 = M 2 1 τ 0 τ 1 -1 ,
which is the limit case (2.14) for (2.12). Or (3.13) holds with a strict inequality, which means that the Lopatinskii determinant ∆(0, •, 1) vanishes in the hyperbolic region. In that case, the strict inequalities (2.12) hold and we necessarily have u = -V , see again [BGS07, chapter 15] for further details. This completes the proof of the 'necessity' part of Theorem 3.1.

We are now going to give a complete construction of a family of Mach stems satisfying (3.1), (3.2) and the causality conditions. We are even going to show that the state U 0 (ε) ahead of the Mach stem can be chosen to have the particular form (3.3) for some smooth function u verifying u(0) = u. This justifies the calculations made in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF]. We assume from now on that the inequalities (2.12) hold and that the tangential velocity u equals -V with V characterized by (2.13). Where all these assumptions come into play will be made precise below. The construction of the family of Mach stems splits in several steps. Let us recall some important facts. The shocks S 1 (ε) and S 2 (ε) will have a 'large' amplitude, meaning that their amplitude will not tend to zero with ε. At the opposite, both discontinuities S 3 (ε) and CD(ε) will have small amplitude, which explains why we deal with those two sets of discontinuities by using separate arguments.

We now wish to construct a one parameter family of Mach stems satisfying (3.1), (3.2). We shall see later on which degrees of freedom on the state U 0 can be spared. The family of Mach stems will be constructed by repeatedly applying the implicit function Theorem (and at the very end of the argument, a somehow degenerate version of that Theorem). Some arguments in the construction are inspired from [Ser00, chapter 4]. We split the whole construction in several steps.

• Step 1. The large amplitude shocks. Let us first set some notation. By applying the implicit function Theorem 12 , we know that for all U close to U 0 and for all ε close to 0, there exist some uniquely determined states U 1 (U ) and U 2 (ε, U ) that are close to U 1 and satisfy the Rankine-Hugoniot relations

f 2 (U 1 (U )) -f 2 (U ) = 0 , f 1 (U 2 (ε, U )) -f 1 (U ) sin ε + f 2 (U 2 (ε, U )) -f 2 (U ) cos ε = 0 , (3.15) U 1 (U 0 ) = U 2 (0, U 0 ) = U 1 .
By uniqueness in the implicit function Theorem, we even have the relation

U 2 (0, U ) = U 1 (U ) ,
for all U sufficiently close to U 0 . The relations (3.15) correspond to admissible discontinuities for the stationary Euler equations across the lines {x 2 = 0} and { sin ε x 1 + cos ε x 2 = 0}, as depicted in Figure 2. At the very end of the analysis, we shall specify the state U close to U 0 by choosing U = U 0 (ε) for some appropriate smooth function U 0 . Then the discontinuity:

U 0 (ε) , x 2 > 0 , U 1 (U 0 (ε)) , x 2 < 0 ,
will satisfy the Rankine-Hugoniot jump conditions for (2.1), and since (U 0 (ε), U 1 (U 0 (ε))) will be close to (U 0 , U 1 ), this discontinuity will necessarily be a shock wave since it will have a nonzero mass flux across {x 2 = 0} and it will satisfy Lax shock inequalities. In the same way, the discontinuity:

U 0 (ε) , sin ε x 1 + cos ε x 2 > 0 , U 2 (ε, U 0 (ε)) , sin ε x 1 + cos ε x 2 < 0 ,
will be a shock wave for (2.1). For both discontinuities, U 0 (ε) will be the upstream state. Up to the choice of the function (ε → U 0 (ε)), which will be specified later on, this means that we have already constructed the two large amplitude shocks of the Mach stem.

• Step 2. The small amplitude shock.

Let us now move on to the small amplitude shock S 3 (ε) in the Mach stem configuration. The arguments follow [Ser00, chapter 4] here. Given the state U 1 with u = -V (and therefore u 2 + v 2 1 > c 2 1 ), we define the angle Ψ 0 ∈ (π, 2 π) by (3.14). In particular, Ψ 0 satisfies the second equation in (3.7) and there holds

u cos Ψ 0 + v 1 sin Ψ 0 = u 2 + v 2 1 -c 2 1 = 0 , use (3.14).
For all pair (U 1 , U 3 ) sufficiently close to (U 1 , U 1 ), we define the matrices

(3.16) A j (U 1 , U 3 ) := 1 0 df j t U 1 + (1 -t) U 3 dt , j = 1, 2 .
Given an angle Ψ, the Rankine-Hugoniot conditions

-f 1 (U 3 ) -f 1 (U 1 ) sin Ψ + f 2 (U 3 ) -f 2 (U 1 ) cos Ψ = 0 , with U 1 = U 3 , equivalently read det -sin Ψ A 1 (U 1 , U 3 ) + cos Ψ A 2 (U 1 , U 3 ) = 0 , (3.17) U 3 -U 1 ∈ Ker -sin Ψ A 1 (U 1 , U 3 ) + cos Ψ A 2 (U 1 , U 3 ) .
Let us recall the expression:

det -sin Ψ A 1 (U, U ) + cos Ψ A 2 (U, U ) = det P (U ) -1 × -u sin Ψ + v cos Ψ -c × -u sin Ψ + v cos Ψ 2 × -u sin Ψ + v cos Ψ + c ,
where the invertible matrix P (U ) is given in (3.9). We can therefore apply the implicit function Theorem and determine a unique angle Ψ(U 1 , U 3 ) that satisfies the eikonal equation (3.17) for any pair (U 1 , U 3 ) sufficiently close to (U 1 , U 1 ), together with Ψ(U 1 , U 1 ) = Ψ 0 . In the case

U 1 = U 3 = U , the angle Ψ(U, U ) is obtained by solving -u sin Ψ + v cos Ψ + c = 0 , with Ψ(U, U ) close to Ψ 0 (u is close to u, v is close to v 1 and c is close to c 1 )
. This means that the angle Ψ(U, U ) is determined by the relations:

(3.18) cos Ψ(U, U ) = -c v + u √ u 2 + v 2 -c 2 u 2 + v 2 , sin Ψ(U, U ) = c u + v √ u 2 + v 2 -c 2 u 2 + v 2 .
For the angle Ψ = Ψ(U 1 , U 3 ), the matrix -sin Ψ A 1 (U 1 , U 3 ) + cos Ψ A 2 (U 1 , U 3 ) has a one-dimensional kernel that is spanned by a vector R(U 1 , U 3 ) which we can choose to satisfy the normalization condition

(3.19) R(U, U ) = R(U ) :=     τ c sin Ψ(U, U ) -c cos Ψ(U, U ) 0     ,
where τ denotes the specific volume for the state U and c the sound speed associated with U . Then for any real number λ close to 0 and any state U 1 close to U 1 , the implicit function Theorem shows that there exists a uniquely determined U 3 (λ, U 1 ) satisfying the Rankine-Hugoniot jump conditions

(3.20) U 3 (λ, U 1 ) -U 1 = λ R U 1 , U 3 (λ, U 1 ) .
In particular, the function U 3 in (3.20) satisfies:

U 3 (0, U 1 ) = U 1 , d U 1 U 3 (0, U 1 ) = I , ∂ λ U 3 (0, U 1 ) = R(U 1 ) ,
for all U 1 close to U 1 (and the vector R(U 1 ) is given in (3.19)). At this point, given any U 1 sufficiently close to U 1 and any λ close to 0, we have constructed a weak solution to (2.1) of the form

U 1 , -sin Ψ(U 1 , U 3 (λ, U 1 )) x 1 + cos Ψ(U 1 , U 3 (λ, U 1 )) x 2 > 0 , U 3 (λ, U 1 ) , -sin Ψ(U 1 , U 3 (λ, U 1 )) x 1 + cos Ψ(U 1 , U 3 (λ, U 1 )) x 2 < 0 .
It is not clear yet whether this discontinuity is a contact discontinuity or a shock wave (at this stage, it could even be a nonadmissible shock, that is it could be a noncharacteristic discontinuity that violates Lax shock inequalities). Ultimately, the state U 1 will be chosen as U 1 (U 0 (ε)) and the state U 3 will therefore be determined as U 3 (λ, U 1 (U 0 (ε))) for some appropriate amplitude λ that will depend on the small parameter ε. As a matter of fact, we are now going to choose the amplitude λ appropriately in order to make the states U 2 and U 3 have equal pressures. 4: Towards the construction of the Mach stem configuration with a shallow angle. Here Ψ is a short notation for Ψ(U 1 (U ), U 3 (λ, U 1 (U ))). At this stage, there is no separation between the states U 2 (ε, U ) and U 3 (λ, U 1 (U )).

x 1 S 2 (ε) x 2 S 1 (ε) U U 1 (U ) U 2 (ε, U ) U 3 (λ, U 1 (U )) S 3 Θ(ε) Ψ Figure
• Step 3. Adapting the small amplitude to get equal pressures for the states 2 and 3. Let us still keep the upstream state U free, close to U 0 . Then we specify the state U 1 as U 1 (U ) and the state U 2 as U 2 (ε, U ). We also specify the state U 3 as U 3 (λ, U 1 (U )), as shown in Figure 4 below.

In order to construct a Mach stem as depicted in Figure 2, we need the states U 2 (ε, U ) and U 3 (λ, U 1 (U )) to have equal pressures and collinear velocities, which will determine the angle Φ as their common argument13 . Let us first adapt the (small) amplitude λ in order to make the pressures of U 2 (ε, U ) and U 3 (λ, U 1 (U )) equal. This relies again on the implicit function Theorem. Indeed, we consider the function

(λ, ε, U ) -→ p U 2 (ε, U ) -p U 3 (λ, U 1 (U )) ,
where p(U ) denotes the pressure associated with a state U (which, of course, only depends on the first and fourth coordinates of U , namely the specific volume and specific entropy). We wish to compute the partial derivative of the above function with respect to λ at (0, 0, U 0 ). We use ∂ λ U 3 (0, U 1 ) = R(U 1 ), with R(U 1 ) as in (3.19), and we find that the partial derivative with respect to λ at (0, 0, U 0 ) equals c 2 1 /τ 1 = 0. By applying the implicit function Theorem, we thus find that for all ε close to zero and for all U close to U 0 , there exists a uniquely determined amplitude λ(ε, U ) close to zero such that

p U 2 (ε, U ) = p U 3 (λ(ε, U ), U 1 (U )) .
Because U 2 (ε, U ) coincides with U 1 (U ) for ε = 0, and because there holds U 3 (0, U 1 ) = U 1 , we find λ(0, U ) = 0 , for all U close to U 0 . It is therefore convenient to rewrite λ(ε, U ) as ε λ(ε, U ) for a new function λ, so that there holds

(3.21) p U 2 (ε, U ) = p U 3 (ε λ(ε, U ), U 1 (U )) ,
for all (ε, U ) close to (0, U 0 ).

• Step 4. Making the velocities collinear. At this point, for all 'bifurcation parameter' ε small enough (ε has no prescribed sign so far), and for all state U close to U 0 , we have constructed:

• the discontinuity S 1 (ε) connecting U to U 2 (ε, U ) (the corresponding angle is Θ(ε) = π -ε as required in (3.2)),
• the discontinuity S 2 (ε) connecting U to U 1 (U ) (the corresponding angle equals zero, which is no loss of generality up to rotating the axes),

• the discontinuity

S 3 (ε) connecting U 1 (U ) to U 3 (ε λ(ε, U ), U 1 (U )) (the corresponding angle equals Ψ(U 1 (U ), U 3 (ε λ(ε, U ), U 1 (U )))).
Whether the discontinuity S 3 (ε) is a shock wave is still undetermined at this point (it is not even known so far whether this discontinuity is admissible in the entropy sense). Let us also recall that the amplitude ε λ(ε, U ) of the discontinuity S 3 (ε) has been tuned so as to make the pressures of the states U 2 (ε, U ) and

U 3 (ε λ(ε, U ), U 1 (U )) equal, see (3.21).
In what follows, we write u(U ) to denote the velocity associated with any state U , that is the vector formed by the second and third coordinates of U . For any sufficiently small ε and any state U close to U 0 , we define the quantity

(3.22) δ(ε, U ) := det u(U 2 (ε, U )) , u(U 3 (ε λ(ε, U ), U 1 (U ))) .
Our goal is to construct a state U 0 (ε), possibly of the form (3.3), such that δ(ε, U 0 (ε)) vanishes for all sufficiently small ε. However, we are facing here a rather degenerate situation because the function δ satisfies (3.23) δ(0, U ) = 0 , for any state U close to U 0 . Therefore there is no chance of proving, for instance, ∂ u δ(0, U 0 ) = 0 and conclude straightforwardly by the implicit function Theorem. As a matter of fact, the cancellation property (3.23) shows that we can write

δ(ε, U ) = ε δ(ε, U ) ,
and our only hope is to apply the implicit function Theorem to the rescaled function δ. There is however a price to pay, which is reminiscent of the analysis in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] and which confirms why u has to be fixed in some very specific way. Let us assume indeed that we are able to construct some smooth U 0 (ε) such that U 0 (0) = U 0 and δ(ε, U 0 (ε)) = 0 for all ε close to 0. Then because the partial derivatives of δ with respect to U at (0, U 0 ) vanish (use (3.23)), we must necessarily have ∂ ε δ(0, U 0 ) = 0. It turns out, see below, that this 'compatibility' condition on the function δ is an equivalent formulation of the constraint u = -V for the tangential velocity of U 0 and U 1 . Let us therefore assume for now that we can prove the relation ∂ ε δ(0, U 0 ) = 0. In view of the factorization δ = ε δ, the relation ∂ ε δ(0, U 0 ) = 0 also reads δ(0, U 0 ) = 0 .

In order to apply the implicit function Theorem to δ, we need to verify that some partial derivative of δ with respect to one of the components of U is nonzero. In what follows, we compute the partial derivative with respect to the second coordinate of U , which corresponds to the tangential velocity. We are going to show ∂ u δ(0, U 0 ) = 0, or equivalently:

∂ 2 ε,u δ(0, U 0 ) = 0 .
This will yield the existence of a state U 0 (ε) of the form (3.3) such that δ(ε, U 0 (ε)) = 0 for all sufficiently small ε. In particular, this will ultimately explain why in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] it was legitimate to expand only the tangential velocity of U 0 with respect to the small bifurcation parameter 14 . Summarizing the above arguments, there are two main points that remain to be clarified. We first need to prove the compatibility condition ∂ ε δ(0, U 0 ) = 0, and we also need to prove the invertibility condition ∂ 2 ε,u δ(0, U 0 ) = 0.

The compatibility condition. Let us compute the partial derivative of δ with respect to ε. Specifying U = U 0 in (3.22), we compute

(3.24) ∂ ε δ(0, U 0 ) = det u 1 , u( U3 -U2 ) ,
where we use the notation

U2 := ∂ ε U 2 (0, U 0 ) = -df 2 (U 1 ) -1 f 1 (U 1 ) -f 1 (U 0 ) , (3.25a) U3 := ∂ ∂ε U 3 (ε λ(ε, U 0 ), U 1 ) ε=0 = λ(0, U 0 ) R(U 1 ) . (3.25b)
By differentiating the relation of equal pressures (3.21), we have

- c 2 1 τ 1 τ U2 -U3 + Γ 1 T 1 s U2 -U3 = 0 ,
where τ U2 -U3 , resp. s U2 -U3 , denote the first, resp. fourth, coordinate of U2 -U3 . This means that the vector U2 -U3 can be decomposed as:

U2 -U3 = μ     Γ 1 T 1 τ 1 0 0 c 2 1     +     0 u U2 -U3 v U2 -U3 0     ,
for some scalar μ, where U3 satisfies (3.25b) and therefore belongs to E -(0, 1) (because the vector R(U 1 ) belongs to E -(0, 1) as seen in the first part of the proof of Theorem 3.1).

14 To our knowledge, it is open, though likely, that the implicit function Theorem could also be applied with respect to the v variable, which would yield another family of Mach stems bifurcating from the same reference step shock (2.3). In any case, there are necessarily infinitely many families of Mach stems that bifurcate from the reference shock (2.3) since ∂ 2 ε,u δ(0, U 0 ) = 0 implies that infinitely many linear combinations of partial derivatives such as µ ∂ 2 ε,τ δ(0, U 0 ) + ∂ 2 ε,u δ(0, U 0 ) (µ small enough) are nonzero.

It appears from (3.24) that the partial derivative ∂ ε δ(0, U 0 ) vanishes if and only if the velocity u( U3 -U2 ) is parallel to u 1 . Hence the partial derivative ∂ ε δ(0, U 0 ) vanishes if and only if the vector U2 -U3 belongs to the vector space

Span            0 u v 1 0     ,     Γ 1 T 1 τ 1 0 0 c 2 1            = E 0 (0, 1) .
Using (3.25a) and (3.25b), this means that the partial derivative ∂ ε δ(0, U 0 ) vanishes if and only if the solution U2 to the linear system

df 2 (U 1 ) U2 = -f 1 (U 1 ) -f 1 (U 0 ) ,
belongs to E s (0, 1), which, as we have already seen, occurs if and only if u = -V . We therefore impose u = -V , which yields ∂ ε δ(0, U 0 ) = 0. It now only remains to verify the invertibility condition ∂ 2 ε,u δ(0, U 0 ) = 0.

The invertibility condition. We shall need in the calculations below the precise expression of several quantities that have arisen earlier. Consistently with the notation in [MR84, page 124], we define

β := cos Ψ 0 ,
where the angle Ψ 0 is defined by (3.14) (we recall that u = -V so we have u 2 + v 2 1 > c 2 1 ). The expression found for β in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF] is exactly the one which we have recalled in (2.16). One can check that it can be equivalently defined by (3.14) provided that u coincides with -V . (The significance of the parameter Φ in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF], which is defined as a root to (2.15), is less clear.) In order to factorize many expressions below, we also introduce the notation (3.26) Υ := u c 1 .

Recalling (3.14), we have

β = M 1 + Υ Υ 2 + M 2 1 -1 Υ 2 + M 2 1 , 1 -β 2 = -Υ + M 1 Υ 2 + M 2 1 -1 Υ 2 + M 2 1 , which, in particular, implies (compare with the definition of c ⋆ in (2.16)) (3.27) Υ = - 1 -M 1 β 1 -β 2 , Υ 2 + M 2 1 -1 = M 1 -β 1 -β 2 .
Since we have u = -V , we know that the Lopatinskii determinant ∆(u, 0, 1) vanishes. As recalled earlier, this means that the linear system (2.9) (with z = 0, η = 1 and χ = 1) has a solution in E s (0, 1). For later purposes, we shall need the expression of that solution, which is nothing but the vector U2 defined in (3.25a). The solution to (2.9) with z = 0, η = 1, and χ = 1, is given by: (3.28)

                     τ = 2 + Γ 1 1 - τ 0 τ 1 v 2 1 u v 0 (c 2 1 -v 2 1 ) (τ 1 -τ 0 ) , u = v 1 -v 0 , v = u c 2 1 -v 2 1 τ 1 τ 0 -1 c 2 1 + v 2 1 + v 2 1 Γ 1 1 - τ 0 τ 1 , ṡ = j 2 T 1 u v 0 (τ 1 -τ 0 ) 2 .
The mass flux j is defined in (2.4). The vector U2 is decomposed on E s (0, 1) as follows:

(3.29) U2 = α 0     0 u v 1 0     + µ 0     Γ 1 T 1 τ 1 0 0 c 2 1     + α -     τ 1 c 1 sin Ψ 0 -c 1 cos Ψ 0 0     , with (3.30) α 0 := - u v 0 (v 1 -v 0 ) 2 u 2 + v 2 1 , α -:= v 1 -v 0 c 1 sin Ψ 0 v 1 v 0 u 2 + v 0 v 1 u 2 + v 2 1 .
The coefficient µ 0 equals ṡ/c 2 1 , with ṡ given in (3.28), but its expression will not be relevant in the subsequent analysis. Recalling the definition (2.17) of the positive parameter ν and using the expression (3.27) of Υ, we can rewrite the coefficients α 0 and α -in (3.30) as

(3.31) α 0 = - M 1 ν 2 1 + ν (1 -M 1 β) 1 -β 2 1 + M 2 1 -2 M 1 β , α -= - M 1 ν 1 + ν 1 1 -β 2 + M 2 1 ν 1 -β 2 1 + M 2 1 -2 M 1 β .
For future use, we note that the vectors U2 , U3 in (3.25) satisfy

(3.32) u( U2 -U3 ) = α 0 u 1 , λ(0, U 0 ) = α -,
where the coefficients α 0 , α -are given in (3.31).

Let us now compute ∂ 2 ε,u δ(0, U 0 ) by differentiating (3.22). Observing that ∂ u U 1 (U 0 ) equals (0, 1, 0, 0), we obtain

∂ 2 ε,u δ(0, U 0 ) = det 1 0 , u( U3 -U2 ) -det u 1 , u(∂ 2 ε,u U 2 (0, U 0 )) + det u 1 , ∂ u λ(0, U 0 ) u(R(U 1 )) + det u 1 , λ(0, U 0 ) u(∂ u R(U 1 )) = -α 0 v 1 -det u 1 , u(∂ 2 ε,u U 2 (0, U 0 )) -c 1 (u cos Ψ 0 + v 1 sin Ψ 0 ) ∂ u λ(0, U 0 ) + α -det u 1 , u(∂ u R(U 1 )) ,
where we have used (3.32) and the expression (3.19) of R(U 1 ). Using now (3.14), we have already simplified the expression of ∂ 2 ε,u δ(0, U 0 ) into:

1 c 1 ∂ 2 ε,u δ(0, U 0 ) = M 1 α 0 - 1 c 1 det u 1 , u(∂ 2 ε,u U 2 (0, U 0 )) -c 1 Υ 2 + M 2 1 -1 ∂ u λ(0, U 0 ) + α - c 1 det u 1 , u(∂ u R(U 1 )) ,
that is, using (3.27),

(3.33) 1 c 1 ∂ 2 ε,u δ(0, U 0 ) = M 1 α 0 - 1 c 1 det u 1 , u(∂ 2 ε,u U 2 (0, U 0 )) - M 1 -β 1 -β 2 c 1 ∂ u λ(0, U 0 ) + α - c 1 det u 1 , u(∂ u R(U 1 )) .
There are three still undetermined quantities on the right hand side of (3.33): the scalar ∂ u λ(0, U 0 ), and the vectors ∂ u R(U 1 ), ∂ 2 ε,u U 2 (0, U 0 ).

The vector ∂ 2 ε,u U 2 (0, U 0 ) is obtained by differentiating the second equation in (3.15) with respect to ε and u, and by therefore solving the linear system:

df 2 (U 1 ) ∂ 2 ε,u U 2 (0, U 0 ) = -∂ u df 2 (U 1 ) U2 + ∂ u f 1 (U 0 ) -∂ u f 1 (U 1 ) .
This means that ∂ 2 ε,u U 2 (0, U 0 ) is a solution to the linear system

df 2 (U 1 ) ∂ 2 ε,u U 2 (0, U 0 ) =       ρ 0 -ρ 1 (ρ 0 -ρ 1 ) u 0 ρ 0 1 2 |u 0 | 2 + e 0 -ρ 1 1 2 |u 1 | 2 + e 1      
, and we eventually obtain the very simple expression

∂ 2 ε,u U 2 (0, U 0 ) = 1 u     τ 0 v ṡ     , v c 1 = -M 1 α 0 + β α -,
where τ , v, ṡ are the first, third and fourth coordinates of the vector U2 and are given in (3.28) or (3.29). This already simplifies (3.33) into (3.34)

1 c 1 ∂ 2 ε,u δ(0, U 0 ) = 2 M 1 α 0 + β α -- M 1 -β 1 -β 2 c 1 ∂ u λ(0, U 0 ) + α - c 1 det u 1 , u(∂ u R(U 1 )) .
Differentiating now the pressure equality (3.21) with respect to ε and u, we get

- c 2 1 τ 2 1 τ (∂ 2 ε,u U 2 (0, U 0 )) + c 2 1 τ 1 ∂ u λ(0, U 0 ) + Γ 1 T 1 τ 1 s(∂ 2 ε,u U 2 (0, U 0 )) = 0 ,
or equivalently (here we use the expression of ∂ 2 ε,u U 2 (0, U 0 ) and (2.9) with z = 0, η = 1 and χ = 1)

c 1 ∂ u λ(0, U 0 ) = 1 u c 1 c 2 1 τ 1 τ -Γ 1 T 1 ṡ = v 1 u c 1 v - v 1 c 1 1 - τ 1 τ 0 = 1 Υ M 2 1 α 0 + M 1 β α -+ M 1 ν 1 + ν .
This simplifies (3.34) into

1 c 1 ∂ 2 ε,u δ(0, U 0 ) = 2 M 1 α 0 + β α -+ M 1 -β 1 -M 1 β M 2 1 α 0 + M 1 β α -- M 1 ν (M 1 -β) (1 + ν) 1 -β 2 + α - c 1 det u 1 , u(∂ u R(U 1 )) ,
Recalling that the vector U 0 (ε) is of the form (3.3), we compute U ′ 0 (0) = 0, u ′ (0), 0, 0 , U ′ 1 (0) = 0, u ′ (0), 0, 0 , U ′ 3 (0) = α -τ 1 , α -c 1 sin Ψ 0 + u ′ (0), -α -c 1 cos Ψ 0 , 0 , where we have used the relation ∂ u U 1 (U 0 ) = (0, 1, 0, 0), (3.32), (3.19) and the properties of U 3 . We know that the states U 0 (ε) and U 1 (ε) satisfy the Rankine-Hugoniot jump relations across S 2 (ε) = {x 2 = 0}. Similarly U 0 (ε) and U 2 (ε) satisfy the Rankine-Hugoniot jump relations across S 1 (ε) = {sin ε x 1 + cos ε x 2 = 0}. Since the reference state (2.3) satisfies (2.4), the discontinuities S 1 (ε) and S 2 (ε) are necessarily shock waves with U 0 (ε) being the state ahead of the shock. We also know from our previous construction that the states U 2 (ε) and U 3 (ε) are connected by a contact discontinuity, which we denote CD(ε), whose angle Φ(ε) is close to Φ 0 . The causality conditions are easily verified since we have

u 1 • (cos Φ 0 , sin Φ 0 ) = |u 1 | > 0 , u 1 • (cos Ψ 0 , sin Ψ 0 ) = u 2 + v 2 1 -c 2 1 > 0 ,
by the definition of the angles Φ 0 and Ψ 0 . By continuity of the states U 1 (ε), U 2 (ε), U 3 (ε) with respect to ε, we thus have u(U 2 (ε)) • (cos Φ(ε), sin Φ(ε)) > 0 , u(U 1 (ε)) • (cos Ψ(ε), sin Ψ(ε)) > 0 , for all ε (up to restricting the interval [-ε 0 , ε 0 ]). The only remaining fact to clarify is to determine whether the discontinuity S 3 (ε) that connects the states U 1 (ε) and U 3 (ε) is a shock wave. This is where the sign condition on ε will arise, as in [START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF]. Let us introduce the function F(Ψ, U ) := det -sin Ψ A 1 (U, U ) + cos Ψ A 2 (U, U ) = det P (U ) -1 -sin Ψ u + cos Ψ v 2 -sin Ψ u + cos Ψ v -c -sin Ψ u + cos Ψ v + c .

This function vanishes at (Ψ 0 , U 1 ) because of the factor -sin Ψ u + cos Ψ v + c. All other factors are nonzero at (Ψ 0 , U 1 ). Differentiating with respect to ε the relation det -sin Ψ(ε) A 1 (U 1 (ε), U 3 (ε)) + cos Ψ(ε) A 2 (U 1 (ε), U 3 (ε)) = 0 , and using the symmetry of A 1 , A 2 with respect to U 1 , U 3 , see (3.16), we end up with the relation

∂ Ψ F(Ψ 0 , U 1 ) Ψ ′ (0) + 1 2 d U F(Ψ 0 , U 1 ) • U ′ 1 (0) + U ′ 3 (0) = 0 .
Using the decomposition of F as a product and the expressions of U ′ 1 (0) and U ′ 3 (0), we obtain the relation (3.37)

u 2 + v 2 1 -c 2 1 Ψ ′ (0) + sin Ψ 0 u ′ (0) = - 1 2 α -c 1 G 1 , G 1 := - τ 1 2
∂ 3 τ τ τ e(τ 1 , s 1 ) ∂ 2 τ τ e(τ 1 , s 1 ) , which connects the derivative Ψ ′ (0) of the angle Ψ(ε) with the first order variation u ′ (0) of the state U 0 (ε).

Let us observe that because of our assumptions on the equation of state, G 1 is a positive quantity (this quantity is a measure of the genuine nonlinearity of the characteristic fields associated with the acoustic waves, see [START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF]).

Figure 1 :

 1 Figure 1: The steady planar shock.

Vectors are written either as rows or columns when no confusion is possible.

The problem studied in[START_REF] Majda | A theory for spontaneous Mach stem formation in reacting shock fronts. I. The basic perturbation analysis[END_REF] deals with reacting flows, for which the pressure laws ahead and behind the shock do not necessarily coincide. The normal mode analysis in[START_REF] Majda | A theory for spontaneous Mach stem formation in reacting shock fronts. I. The basic perturbation analysis[END_REF] though is independent of the presence of a chemical reaction as long as Lax shock inequalities (2.4c) below are satisfied. We thus restrict to the more standard framework of the compressible Euler equations without reaction and with same pressure law on either side of the shock for simplicity.

Here we have extracted the i factor from all frequency parameters, and consider the half space x2 ≤ 0, which is the reason why we consider Im z ≤ 0 and Im ω ≤ 0. The conventions in[START_REF] Zumbrun | Stability of large-amplitude shock waves of compressible Navier-Stokes equations[END_REF][START_REF] Benzoni-Gavage | Multidimensional hyperbolic partial differential equations[END_REF] are different but passing from one to the other is harmless.

Here we have kept track of the possibly nonzero tangential velocity u while linearizing the Rankine-Hugoniot conditions and made some elementary linear combinations between several equations to get (2.9).

As explained in[START_REF] Benzoni-Gavage | Multidimensional hyperbolic partial differential equations[END_REF], there are many possible ways to define ∆ but all possible definitions give rise to the same zeroes, if there are any.

We refer to V as a velocity since it has the physical homogeneity of a velocity, see equation (2.13)-(iii).

The final conclusion of the shock wave stability analysis is written here in two space dimensions. In higher space dimensions, one should rather write z = ±V |η| instead of z = ±V η (η becomes a real vector in dimension d ≥ 3.)

Up to a rotation, the shock S2(ε) is always chosen as the (0x1) axis, which yields the convention for the angles Φ0 and Ψ0. The location of the contact discontinuity with respect to S3(ε) follows from the convention u ≤ 0, the symmetric situation being obtained by changing u into its opposite and x1 into -x1.

In[START_REF] Majda | A theory for spontaneous Mach-stem formation in reacting shock fronts. II. Steady-wave bifurcations and the evidence for breakdown[END_REF], the state U ′ 0 is assumed to have only one nonzero component, which corresponds to the tangential velocity. This assumption will be justified when we construct the family of Mach stems but it is actually not needed at this point of the argument.

This means that the velocity field u1 behind the shock in (2.3) is supersonic though the normal velocity v1 is subsonic by Lax shock inequalities.

This limit case corresponds to a Lopatinskii determinant that vanishes at a glancing frequency. At such frequencies, the eigenmode ω-is not locally a smooth function of (z, η).

For simplicity, we assume that the specific internal energy e is a C ∞ function of (τ, s) so all functions involved in the calculations below have C ∞ regularity.

Both velocities will be small perturbations of u1 = (u, v1) so collinearity here means proportionality with a positive scaling factor, which will uniquely determine the angle Φ close to Φ0 ∈ (π, 3 π/2).

and the author gratefully thanks this institution for the stimulating working conditions it has provided. Research of J.-F. C. was supported by ANR project BoND, ANR-13-BS01-0009, and ANR project Nabuco, ANR-17-CE40-0025.

Solving (3.12) for cos Ψ 0 , we obtain the two possible values

both values belonging to the interval (cos Φ 0 , 1). Since both possible values correspond to a positive value of c 1 + v 1 cos Ψ 0 , we obtain

where the choice between ± is the same in both expressions of cos Ψ 0 and sin Ψ 0 . At this stage, it is unfortunately still not possible to select the appropriate determination of Ψ 0 between the two possible choices. However, the causality conditions for the family of Mach stems imply

so passing to the limit, we must have

Therefore the appropriate determination of Ψ 0 corresponds to

which amounts to choosing the smallest possible value of cos Ψ 0 among the two possible ones, see the geometric interpretation in Figure 3. that is, after a little bit of factorizations

Using (3.18), we now differentiate the definition (3.19) with respect to u and get

Incorporating this last expression in (3.35), we obtain

We now use the expressions of α 0 and α -in (3.31) and derive the final expression

.

Recalling that M 1 belongs to (0, 1) and that β belongs to (-1, M 1 ), we see that up to a nonzero multiplicative factor, ∂ 2 ε,u δ(0, U 0 ) can we written as Ω 0 + ν M 1 Ω 1 where both Ω 0 and Ω 1 are positive quantities. Consequently, we have ∂ 2 ε,u δ(0, U 0 ) = 0 and we can apply the implicit function Theorem to the function δ. We can therefore construct a smooth function U 0 (ε) of the form (3.3) such that δ(ε, U 0 (ε)) = 0 for all ε close to 0.

Since the velocities of the two states U 2 (ε, U 0 (ε)) and U 3 (ε λ(ε, U 1 (U 0 (ε))), U 1 (U 0 (ε))) are collinear and close to u 1 , we can determine the angle Φ(ε), close to Φ 0 , as their common argument and therefore determine the location of the contact discontinuity CD(ε).

• Step 5. Conclusion.

It remains to verify that the corresponding four wave pattern which we have just constructed defines a Mach stem for any positive ε. For the sake of clarity, let us denote

We also define the angle

which represents the location of the discontinuity between U 1 (ε) and U 3 (ε). All these functions are defined on a common interval [-ε 0 , ε 0 ] and are smooth functions of ε on that interval.

We can now verify that the discontinuity S 3 (ε) is a shock wave for any sufficiently small ε > 0. Indeed, the only thing to verify is that Lax shock inequalities

are satisfied for any sufficiently small ε > 0. Using the expressions of U ′ 1 (0) and U ′ 3 (0), we compute

where we have used (3.37), and we also compute 15

where we have used (3.37) again. We observe from the expression (3.31) that α -is negative (recall that the parameter ν in (2.17