Stability of transonic jets with strong rarefaction waves for two-dimensional steady compressible Euler system

We study supersonic flow past a convex corner which is surrounded by quiescent gas. When the pressure of the upstream supersonic flow is larger than that of the quiescent gas, there appears a strong rarefaction wave to rarefy the supersonic gas. Meanwhile, a transonic characteristic discontinuity appears to separate the supersonic flow behind the rarefaction wave from the static gas. In this paper, we employ a wave front tracking method to establish structural stability of such a flow pattern under non-smooth perturbations of the upcoming supersonic flow. It is an initial-value/free-boundary problem for the two-dimensional steady non-isentropic compressible Euler system. The main ingredients are careful analysis of wave interactions and construction of suitable Glimm functional, to overcome the difficulty that the strong rarefaction wave has a large total variation.

We are concerned with two-dimensional steady non-isentropic compressible supersonic Euler flows passing a convex corner surrounded by static gas. When the pressure of the upcoming supersonic flow is larger than that of the quiescent gas, there may appear a strong rarefaction wave to rarefy the upcoming flow to the quiescent gas. Meanwhile, a characteristic discontinuity, which is a combination of vortex sheet and entropy wave, is generated to separate the supersonic flow behind the rarefaction wave from the still gas, see Figure 1. Under suitable assumptions on the upstream supersonic flow and the surrounding quiescent gas, we wish to establish structural stability of such a flow pattern in the class of functions with bounded variations by considering an initial-value/free-boundary problem for the two-dimensional steady full compressible Euler system.
Recall that the Euler system, consisting of conservation of mass, momentum and energy, reads as where ρ, p and (u, v) represent respectively the density of mass, scalar pressure and velocity of the flow, and E = 1 2 (u 2 + v 2 ) + e is the total energy per unit mass, with e the internal energy. Let S be the entropy. For polytropic gas, the constitutive relations are given by with κ, c v , γ > 1 being positive constants. The corner and the Cartesian coordinates (x, y) of the plane are illustrated in Figure 1.
System (1) can be written in the general form of conservation laws: where W (U ) = ρu, ρu 2 + p, ρuv, ρu(E + p ρ ) ⊤ , H(U ) = ρv, ρuv, ρv 2 + p, ρv(E + p ρ ) ⊤ . The eigenvalues of system (2), namely, the roots of the polynomial det(λ∇ U W (U )− ∇ U H(U )), are provided that u > c > 0, where c = γp/ρ is the local sonic speed. The corresponding right eigenvectors (null vectors of the matrix λ∇ U W (U ) − ∇ U H(U )) are So for supersonic flow with u > c > 0, all the characteristics are real and the four eigenvectors are linearly independent. Hence by definition, the system (1) is hyperbolic in the positive x-direction if u > c > 0, with constant characteristic multiplicities. We may introduce the constants k j > 0 to normalize r j so that ∇ U λ j (U ) · r j (U ) ≡ 1 for j = 1, 3 (see the Appendix for detailed computations), which means the first and the third characteristic families are genuinely nonlinear. The second characteristic family is linearly degenerate, namely ∇ U λ 2 (U )·r 2k (U ) ≡ 0 for k = 1, 2.
Next we formulate the domain and boundary conditions to describe the transonic jet problem. Since we consider only time-independent flow, it is reasonable to assume that the state below the characteristic discontinuity is always static, and unchanged even if the upcoming supersonic flow experienced perturbations. The static state is denoted byŪ = (0, 0,p,ρ) ⊤ . The characteristic discontinuity C itself is a free-boundary, which is the graph of a Lipschitz continuous function y = g(x) on x ≥ 0, with g(0) = 0. So the domain we consider is Ω g = {(x, y) ∈ R 2 : x > 0, y > g(x)}.
From the Rankine-Hugoniot jump conditions across a characteristic discontinuity, we have the following two boundary conditions wherep is the constant pressure of the quiescent gas. The second condition implies that C is a characteristic boundary for the Euler equations [4, p.3]. Regarding x as time, suppose that the supersonic flow on I = {(x, y) ∈ R 2 : x = 0, y > 0} is given by U (0, y) = U 0 (y), with u 0 (y) > c 0 (y), which is the initial data. Hence the main problem we will study in this paper is the Euler system (1) in the unknown domain Ω g , subjected to the following initialvalue/free-boundary conditions (u, v, p, ρ)(0, y) = (u 0 (y), v 0 (y), p 0 (y), ρ 0 (y)), for y ≥ 0; p =p, g ′ (x) = v u (x, g(x)), on y = g(x).
We will seek global entropy solutions of this problem.
(ii) The following entropy inequality holds in the sense of distributions: In the next section, we will show that for suitably chosen constant supersonic state U 0 (y) = U + and static gasŪ , there exists an entropy solution (U b , g b ) to problem (2)(4), which consists of a strong rarefaction wave of the third characteristic family to decrease the pressure, and a transonic characteristic discontinuity y = g b (x) of the second characteristic family to separate the constant supersonic state U − behind the rarefaction wave from the static gasŪ . Such a special solution is called a background solution in the sequel. The purpose of this work is to show structural stability of these background solutions, which is described in the following theorem.
Theorem 1.2. For given constant supersonic state U + = (u + , 0, p + , ρ + ) ⊤ , there is a positive number p * . Letp be the pressure of the static gasŪ , with p * <p < p + . Then there exists a background solution where TV.U 0 (·) denotes the total variation of the vector-valued function U 0 (y), then problem (2)(4) admits an entropy solution (U (x, y), g(x)), containing a 3-strong rarefaction wave, which is a small perturbation of the background solution, in the sense that for almost all x ≥ 0, it holds that and for a.e. (x, y) ∈ Ω g , (20) in §2, which is a neighborhood of the background 3-rarefaction wave curve in the state space {U ∈ R 4 }, and is an invariant region for this problem.
We now review some known results on the stability of large solutions for the stationary compressible Euler equations. The transonic characteristic discontinuity separating supersonic flow from static gas was firstly studied in [4,5]. The case that the jet contains a strong shock was solved in [16]. This paper is a continuation of these works. The new feature here is the appearance of strong rarefaction waves. To our knowledge, Zhang firstly studied an initial-boundary value problem for the isentropic irrotational Euler equation with strong rarefaction waves in [22]. Later, in [3], the authors studied strong rarefaction wave in steady exothermically reacting Euler flows, and in [12], the authors considered a piston problem for the case that the piston was drawn away from the gas which filled a straight thin tube ahead of the piston. There appears a strong rarefaction wave in the tube. See also [10,11,13] for some generalizations. There are also many works on flow fields containing strong supersonic shocks. Except [16], one may consult, for instance, the work of Wang-Zhang [21] on steady supersonic flow past a curved cone, and Chen et.al. [7] on supersonic flow past a wedge. Considerable progress has been made on the existence and stability of transonic shocks in steady full Euler flows (see, for example, [17] and references therein), which, unlike our case, are treating non-characteristic freeboundaries. Structural stability of supersonic characteristic discontinuities over Lipschitz walls were thoroughly investigated by Chen et.al. in [6]. There are also some results on transonic characteristic discontinuities for three-dimensional steady compressible Euler system, see [20,18].
In this paper we will employ a wave front tracking algorithm to construct a family of approximate solutions to problem (2)(4). Then by modifying appropriately the Glimm functional introduced in [14,2] and showing its monotonicity, we obtain compactness in the space of functions with bounded variations so that a subsequence of the approximate solutions converges to an entropy solution. One of the main difficulties here is to obtain rather accurate estimates when weak waves interact with strong rarefaction waves. One needs to introduce appropriate interaction potentials to take into account of the fact that the total variation of the strong rarefaction wave is not small, and choose carefully weights for various terms in the Glimm functional. Some of the ideas were inspired by [1,22]. However, different from previous works with strong rarefaction waves, in our model a characteristic free-boundary, namely the strong transonic characteristic discontinuity occurs. Therefore, the analysis is somewhat subtle and different. We recommend [2,15] for a general introduction of the front tracking method. This paper is organized as follows. In §2, we review some basic properties of elementary waves of the steady Euler system. Then we establish existence of background solutions and show solvability of free-boundary Riemann problems for system (1). In §3, we outline the construction of approximate solutions by wave front tracking method. In §4, we analyze local interaction estimates of perturbed waves and reflections of waves on the transonic characteristic discontinuity in the front tracking process. Then we construct a Glimm functional and prove its monotonicity. This manifests that the approximate solutions can really be computed for all x > 0. In §5, we show consistency, namely the limit of the approximate solution is actually an accurate entropy solution. §6 is a short appendix, contains some quantities used in the main text.
2. Riemann Problems and Background Solutions. In this section, we find a background solution to problem (2)(4) with constant initial data, and study reflections of waves from the characteristic free-boundary.
2.1. Wave curves in state space. Firstly we consider Riemann problem of (2) with initial data: where and right (upper) states, respectively. Solvability of the Riemann problem for general strictly hyperbolic conservation laws with genuinely nonlinear or linearly degenerate characteristics can be found in [2,15,19,9] when |U L − U R | is sufficiently small. The basic idea is to introduce several wave curves in the state space, which can also be used to construct large solutions, as shown below. For any given left (resp. right) state U l (resp. U r ), the set of all possible states U which can be connected to U l (resp. U r ) on the right (resp. left) by 1-or 3-shock wave, is denoted by S 1 (U l ) or S 3 (U l ) (resp. S −1 1 (U r ) or S −1 3 (U r )). Similarly, we denote R 1 (U l ) or R 3 (U l ) (resp. R −1 1 (U r ) or R −1 3 (U r )) the (inverse) wave curves of 1-or 3-rarefaction wave. We can parameterize R j (U l ) (j = 1, 3) by For our use, the 1-inverse rarefaction wave curve R −1 1 (U r ) is where q .
= arctan v u , and Similarly, the inverse wave curve of 3-rarefaction wave R −1 3 (U r ) is given by We remark that the above representation of the rarefaction wave curve is not complete. It only contains the part where u > 0. The second characteristic field is linearly degenerate. So it only supports characteristic discontinuities. The corresponding wave curve passing through U l for vortex sheet (i.e. integral curve of vector field r 21 in state space R 4 ) is 1 and the entropy wave curve (integral curve of r 22 ) through U l is given by The two waves coincide in the physical (x, y)-plane, with the same line {(x, y) : The Rankine-Hugoniot jump conditions across the j-shock give the wave curves S −1 j (U r ): where [h] = h r − h stands for the jump of a quantity h across a shock front; s j is the speed of the shock front: and σ(1) = 1, σ(3) = 0. The entropy inequality means the pressure increases when particles pass a shock front, namely [p] > 0 for 1-shock and [p] < 0 for 3-shock (cf. [15, p.276]). A solution to the Riemann problem (2)(9) is given by at most four constant states connected by shocks, characteristic discontinuities, and/or rarefaction waves. Exactly speaking, there exist piecewise C 2 curves α j → Φ j (α j ; U ), j = 1, 3, and a C 2 surface α 2 . = (α 21 , α 22 ) → Φ 2 (α 2 ; U ) (the latter is the point with parameter α 22 on C 22 (U m ), while U m is the point with parameter α 21 on curve C 21 (U )), such that whenever |U L − U R | ≪ 1.

Remark 1.
Hereinafter, we denote by α i , β i , γ i etc. the parameters for the iwave curve/surface, i = 1, 2, 3, while by their absolute values the corresponding strengths of the waves. It should be noted that since the Euler system is not strictly hyperbolic, we introduce the convention that strength of the characteristic discontinuity with parameter α 2 is |α 2 | . = |α 21 | + |α 22 |. We also use the parameters to represent the i-waves.
To study reflection of waves on free-boundary, we introduce the notation U l = Ψ(α 1 , α 2 , α 3 ; U r ) to represent the inverse wave curves; that means, the left state U l 8 MIN DING AND HAIRONG YUAN and the right state U r can be connected by 1-wave α 1 , 2-wave α 2 and 3-wave α 3 . From the above constructions, we have where α 2 = (α 21 , α 22 ), and r i , r 2k are the right eigenvectors of the system (1), given by (3).
. Particularly, we can parameterize the 3-inverse rarefaction wave curve R −1 3 (U ) by solving 2.2. Background solutions. We now consider problem (2)(4) with uniform upcoming supersonic flow whose pressure is quite large.
Here k b , k 1 , k 2 are constants, while U ba ( y x ) is a 3-rarefaction wave connecting the constant states U + and U − , with p − =p. Furthermore, there exists a positive constant C such that Proof. From the second and third equations in (11), the entropy S and Bernoulli constant B are invariant across 3-rarefaction waves, hence they are the same as those of U + . Set M 1 = u/c. It is a function of q, namely, By Bernoulli law, we also have . It is obvious that χ(p + ) = u + /c + > 1. By continuity, there is a number p * > 0 so that M 1 = χ(p) > 1 for p ∈ (p * , p + ]. Now by requiring p − =p > p * , and p − < p + , we determine q − = χ 2 (p − ) > q + . Using the first equation in (11): we solve θ − , which is negative, hence k b = tan θ − and we find the downstream state U − in (18). Then k 1 = λ 3 (U + ) and k 2 = λ 3 (U − ). For given θ ∈ (θ − , 0), we solve the speed q ♯ and the velocity (u ♯ , v ♯ ) = (q ♯ cos θ, q ♯ sin θ) from and hence determine the state U ♯ (θ). Then U ba (y/x) = U ♯ (θ), where we could solve θ from λ 3 (U ♯ (θ)) = y/x, by virtue of the genuine nonlinearity of λ 3 . The estimate (19) follows from monotonicity of p, ρ, q and constancy of S, B (see (10)) along the 3-rarefaction wave curve (cf. Lemma 2.2 below). Since the pressure is decreasing along rarefaction waves, one checks easily that u − > c − by our choice of p * and p − .

Remark 2.
We observe that no vacuum could appear. This is different from the piston problem.
For small δ 0 > 0 (which is then fixed so that all the Riemann problems are solvable in the sequel), we introduce the perturbation domain D(U + , δ 0 ) in state space as follows: which may be considered as a neighborhood in R 4 of the 3-rarefaction wave curve passing through U + , with pressure lying in [p, p + ]. The following lemma shows we could use the variation of pressure to measure the strength of large rarefaction waves in D(U + , δ 0 ), just as in (19). We use V (k) to denote the k-th argument of a vector V .
Moreover, there exist two positive constants C 1 and C 2 such that Proof. For α 3 ≥ 0, by the properties of the rarefaction waves, we have for any U l and Φ 3 (α 3 ; U l ) ∈ D(U + , δ 0 ). So by continuous differentiability of the curve Φ 3 (α 3 ; U l ), there is a δ ′ * > 0 so that for −δ ′ * ≤ α 3 < 0, we still know the derivative calculated above is positive. The estimate (21) follows from the mean value theorem and inverse function theorem of differentiable monotonic functions.
The following lemma shows that strength of the rarefaction wave in the background solution is Lipschitz continuous with respect to U + . Figure 2. Reflection on the free-boundary.
Then there is a constant C depending only on δ 0 so that Proof.
Thus, by the inverse function theorem of differentiable monotonic functions, there exists a C 1 function h such that and mean value theorem implies (22).

2.3.
Reflection of waves on the characteristic free-boundary. Assume that a 1-wave front α 1 hits the transonic free-boundary at (x 0 , y 0 ) (cf. §3.1 below that a rarefaction wave is replaced by several fronts in the (x, y)-plane, and we need only consider one front meets the free-boundary here, as shown in Figure 2.2). Let U l and U r be the left and right states of α 1 , respectively. The reflected wave is a 3-wave, denoted by β 3 . Suppose the left and right states of β 3 are U m and U r , respectively. Consider the free-boundary Riemann problem where k is a constant to be solved.
where K b1 is positive and uniformly bounded, with a bound C b depending only on the background solution and δ 0 .
We recall that δ 0 appeared in (20), which is small and fixed.
Proof. By the boundary conditions on the free-boundary, we have Differentiate (25) with respect to β 3 , and let β 3 = 0 (hence α 1 = 0), we see By the implicit function theorem, close to α 1 = 0, there exists a function f 1 ∈ C 1 such that Hence U m is solved by Now consider the left-hand side of (25) as a function of α 1 , and differentiate it with respect to α 1 , we have Let α 1 = 0. Then From Taylor's formula, we have The explicit expressions of k 1 and k 3 in (83) were used here to determine the sign. Since K b1 is continuous with respect to α 1 and U r , for ε 0 > 0 small, if |α 1 | < ε 0 and |U r − U − | < ε 0 , we still have Thus, K b1 is positive and bounded by a constant C b depending only on the background solution and δ 0 .
Remark 3. This lemma says that 1-wave front is changed to 3-wave after reflection, and rarefaction front becomes shock front, and vice versa. It shows how to solve the free-boundary as "times" evolves. We remark that U m ∈ B(U − , M ε) implies that |k − k b | ≤ M 0 ε, for a constant M 0 depending only on the background solution, by noting that the velocity u has a positive lower bound in D(U + , δ 0 ), and using the mean value theorem for continuously differentiable functions. This fact finally leads to (6) claimed in Theorem 1.2.
3. Construction of Approximate Solutions. In this section, following the general ideas presented in [2,15], we define accurate and simplified Riemann solvers, which are both appropriate modifications of the exact solutions of the Riemann problems indicated in §2. They are building blocks to construct approximate solutions of problem (2)(4)  Case 1. Accurate Riemann solver Let δ > 0 be given, which measures accuracy of each approximate solution constructed by front tracking later. The accurate Riemann solver is as mentioned in §2, except that every rarefaction wave R i , i = 1, 3, with strength α i > 0, has been divided into ν equal parts and replaced by ν rarefaction fronts. Here ν is the smallest positive integer larger than or equal to α i /δ.
More precisely, suppose that the left state U L and the right state U R are connected by a 1-rarefaction wave α 1 . Let U 0,0 = U L , U 0,ν = U R , and for any 1 ≤ k ≤ ν, set Then we replace the 1-rarefaction wave by piecewise constant states for x > x 0 : where λ * 1 ∈ (max U∈D(U+,δ0) λ 1 (U ), min U∈D(U+,δ0) λ 2 (U )) is a fixed number. Similarly, we can approximate 3-rarefaction wave by ν 3-rarefaction fronts in the domain

Case 2. Simplified Riemann solver
Letλ (strictly larger than all the eigenvalues of system (2) in D(U + , δ 0 )) be the speed of artificial discontinuities called non-physical fronts, which are introduced so that the total number of wave fronts (discontinuities) is finite for all x ≥ 0 for a given approximate solution of problem (2)(4). The strength of a non-physical wave is the error due to the following simplified Riemann solver. It occurs in three cases: Case a. A j-wave β j and an i-wave α i interact at (x 0 , y 0 ), 1 ≤ i ≤ j ≤ 3. Suppose that U L , U M and U R are three constant states, satisfying We define an auxiliary right state The simplified Riemann solver U S (U L , U R ) at (x 0 , y 0 ) of problem (2)(9) is is constructed by the accurate Riemann solver as in Case 1 (the rarefaction waves are split, while all shocks and characteristic discontinuities remain unchanged). The non-physical front is defined by and the strength of the non-physical front is ǫ = |U R − U ′ R |. Case b. A non-physical front interacts with an i-wave front α i (i = 1, 2, 3) coming from the above/right at (x 0 , y 0 ). Suppose that the three states U L , U M and U R satisfy Then the simplified Riemann solver U S (U L , U R ) of problem (2)(9) is where x > x 0 . This means that strength of the physical front is unchanged after interaction, namely still to be α i , while strength of the non-physical front becomes to be ǫ ′ = |Φ i (α i ; U L ) − U R |. In particular, if α i is a rarefaction front, there is no splitting of rarefaction waves in the solver U δ A (U L , Φ i (α i ; U L )), since |α i | ≤ δ. Case c. A 1-wave front α 1 hits the free-boundary. In Lemma 2.4 we have shown the accurate Riemann solver, where β 3 is the reflected wave, separating U m and the new free-boundary from the states U r behind α 1 . For the simplified Riemann solver, we just replace β 3 by a non-physical front ǫ travelling with speedλ. So the state ahead and behind of ǫ is still U r and U m respectively, and the free-boundary is unchanged. By Lemma 2.4, we have the estimate The non-physical fronts will be considered as fronts of the fourth family in the rest of the paper.

Approximate solutions.
For any sufficiently small δ > 0, we construct a δ-approximate solution (U δ (x, y), g δ (x)) to (2)(4) by induction in the region {x > 0, y ∈ R} as follows: Step 1. For x = 0, we approximate the initial data U 0 by U δ (0, y) (y ≥ 0), a piecewise constant function with finite jumps. It is required that Note that our assumptions in Theorem 1.2 imply that lim y→+∞ U 0 (y) = U + . So U δ (0, y) = U + for large y. The number ε 0 (see Theorem 1.2) is chosen small so that all the standard Riemann problems at the discontinuous points of U δ (0, y) are solvable. At the corner (0, 0), as in Lemma 2.1, we solve a free-boundary Riemann problem with a strong rarefaction wave.
Then we approximate all the rarefaction waves appeared in these Riemann problems by rarefaction fronts as described by Case 1 in §3.1.
Step 2. By induction, we assume that (U δ , g δ ) has been constructed for x < τ , for some τ > 0, and assume that U δ | x<τ consists of a finite number of wave fronts and for the first time, some of them interact, at x = τ . As shown in §3.1, we solve a Riemann problem when two wave fronts interact, or a free-boundary Riemann problem when a wave front hits the free-boundary. Thus we extend the approximate solution (U δ , g δ ) beyond x = τ .

Remark 4.
We may adjust speeds of wave fronts by a quantity arbitrarily small, so that there are no more than two discontinuities intersect at a point (τ, y 0 ), and only one wave front hits the free-boundary for each "time" x > 0. Also, at any "time" x = τ , only one interaction occurs [2, p.132, Remark 7.1]. Otherwise, we need technically more complicate estimates of wave interactions as in [15, p.290]. Also, by our choice of speedλ, a non-physical front will never hit the free boundary.
To distinguish those fronts obtained by splitting the strong rarefaction wave that issued from the corner and the other weak rarefaction fronts coming from perturbations of the initial date, we assign to each front of an approximate solution an integer called generation order in the following way, cf. [15, p.300]. This is also used as a tool to estimate the total strength of non-physical fronts, at any given "time" x, of an approximate solution. (A) All wave fronts issued from x = 0, y > 0 have order 1; the free-boundary, as well as all the rarefaction fronts obtained from splitting the rarefaction waves issued from the corner (0, 0), have order 0. (B) A wave front of order k hits the free-boundary, then the generation order of the reflected fronts are still k. (C) An i-wave front α i of order k 1 interacts with a j-wave front β j of order k 2 at a point (τ, y 0 ). We then assign generation order to the produced l-wave fronts (l = 1, 2, 3, 4) to be Definition 3.1. (Strong rarefaction front) A front s is called a strong rarefaction front provided that s is a 3-rarefaction wave front with generation order 0. Otherwise, it is called a weak front.
There may appear weak front of generation order zero, but all non-physical fronts have order at least one.
We now indicate what solver is used to solve Riemann problems encountered in Step 2 above. To be more specific, there is at most one of the following six cases occurring at the interaction "time" x = τ : Case 1. Two weak fronts α i and β j (1 ≤ i, j ≤ 3) interact; Case 2. A 1-weak front α 1 hits the characteristic free-boundary; Case 3. A strong rarefaction front s interacts with an i-weak wave front α i (i = 1, 2) from the above at "time" x = τ ; Case 4. A strong rarefaction front s interacts with a 3-weak shock front α 3 from the below (above) at "time" x = τ ; Case 5. A strong rarefaction front s interacts with a non-physical front ǫ; Case 6. A non-physical front ǫ interacts with an i-weak wave α i from the above (i = 1, 2, 3).
For the given approximate solution (U δ , g δ ), set Here µ δ is a small parameter depending on δ, which is to be specified in (78).
Rule 2. For Case 5-Case 6, we always use the simplified Riemann solver.

Remark 5.
To make sure we could construct the approximate solution (U δ , g δ ) in {0 ≤ x < T } for any given large T , by the above front tracking algorithm, we need to guarantee the following: (A). For any τ > 0, U δ (τ ) ∈ D(U + , δ 0 ), and the perturbation of total variation is small, so each Riemann problem could be solved; (B). At any "time" x = τ , the total number of physical and non-physical front is finite, with a bound independent of τ . So only a finite number of interactions occur in the approximate solution.
We notice that (A) is demonstrated in the next section by showing a Glimm functional F (τ ) is decreasing, if the total variation of the initial perturbation is small, by specifying ε 0 in Theorem 1.2. Then (B) follows easily. Furthermore, to prove consistency of the approximate solutions, the key point is to show that the total strength of the non-physical front at any non-interaction "time" x = τ is of the order O(δ). This is achieved later by choosing carefully µ δ , see (78).
To end this section, we state the following lemma, which says how to estimate the distance from a state U 3 by the 3-rarefaction wave curve to a state U 1 . It can be used to prove (8) claimed in Theorem 1.2.
Proof. From the expression of rarefaction wave curves, for any α 3 ≥ 0, we have For α 3 < 0 and other waves, we use the standard results, namely, employ wave curve/surface parameters α i , i = 1, 2, 3 to measure variations of physical states (see, for example, [15, p.256, (5.144)]). This completes proof of the lemma.

Monotonicity of Glimm Functional.
In this section, we construct a Glimm functional and prove its monotonicity based on local interaction estimates. Then by induction we demonstrate that for any δ > 0 small, we could construct a global approximate solution (U δ , g δ ) as outlined in §3.2.

A Glimm functional.
For convenience of writing, for any weak front α, denote its position and magnitude by y α (x) and α, respectively. Similarly, for a front s of the strong rarefaction wave, denote its location and magnitude by y s (x) and s. We also define approaching waves as follows.
The concept here of approaching waves is the same as that introduced by Glimm in [14]. We next introduce the most important functionals to study strong rarefaction waves. For any weak front α and any non-physical front ǫ, denote R(x, α, b) = {s|s is a front of the strong rarefaction wave with y s (x) < y α (x)}, R(x, ǫ, a) = {s|s is a front of the strong rarefaction wave with y s (x) > y ǫ (x)}, We observe that R(x, α, b) is the set of those strong rarefaction fronts that lie below the weak front α at "time" x, and R(x, ǫ, a) is the set of strong rarefaction fronts that lie above the non-physical front ǫ at "time" x. The idea is, for example, ǫ is approaching all the strong 3-rarefaction fronts in R(x, ǫ, a). The crucial functionals W (α, x) and W (ǫ, x), utilize fast growth of the exponential functions, could drastically magnify the decrease of strength of the strong rarefaction wave, when a weak front penetrating into it, by further choosing the constants K ω and K np large, thus overtaken the difficulty that the 3-strong rarefaction wave has a large total variation, for which the original Glimm interaction potential (see Q 0 below) fails.
We further set The two functionals are used to control the total variation of all the weak waves introduced by the perturbations of the initial data (excluding the corner) in our problem. We also need the following functionals: which represent respectively the interaction potentials between weak physical wave fronts and/or non-physical fronts, a weak wave front and the 3-strong rarefaction waves, a non-physical front and the 3-strong rarefaction waves. The functional Q 0 was originally introduced by Glimm, while Q i (i = 1, 2, 4) resemble those appeared in [2, p.139, (7.65)], and has been used in many previous works [3,10,12,22]. We do not need the interaction potential between 3-weak shocks and the strong 3rarefaction waves, since for this case there are cancellations and it is not necessary to consider second-order terms in interactions of waves. Define S > 0 to be the strength of the background strong rarefaction wave, namely Φ 3 (S; U − ) = U + . Then we set {s(x) : s is a strong 3-rarefaction wave front}, The latter is used to measure the perturbation of the 3-strong rarefaction wave at each "time" x.

Finally we introduce
and the Glimm functional is defined as where K, K 0 , K 1 , K 2 , K 3 , K 4 , K ω , K np and K * are positive constants called weights that need to be chosen later (cf. (59)-(67)). The weight K is used to handle the reflections of 1-wave fronts on free-boundary, while K 3 is used to magnify the cancellation between 3-weak shock fronts and 3-strong rarefaction waves. Although it turns out that we may take K 1 = K 2 = K 4 = K * = 1 later, we retain them for easy to track estimates of each term in the Glimm functional in the following computations.

4.2.
Decreasing of Glimm functional. Note that the Glimm functional experiences changes only if two fronts interact, or a physical front hits the free-boundary, at some interaction "time" x = τ . We have the following crucial result.
To prove (37), we now check the six cases listed before in §3.2.

Case 1. Interaction between weak physical fronts.
Let the two weak fronts α i and β j interact at a point on the line x = τ , and γ l be the generated waves, l = 1, 2, 3, and ǫ be the outgoing non-physical front (if we use the simplified Riemann solver). By a standard procedure (see [2, p.133] or [15, p.290]), we have the following estimates, even if we use the convention made in Remark 1 on strengths of 2-waves.

Lemma 4.3. It holds that
and • when i = j, • when i = j, All the quantities O(1) here are bounded in D(U + , δ 0 ) with a uniform bound C 1 .
Without loss of generality, we may assume in the following that C 1 ≥ 1. Based on the estimates (38)-(41), we have and no matter the interaction happens above/below/in the middle of the strong rarefaction waves, it always holds that Using now standard arguments as in [15, pp.294-295], and assumption (36), we have By bounds of S(τ −), we also get It follows that where we choose K 0 large enough so that Case 2. Reflection of a front on the free-boundary.
Assume that a weak 1-wave front α 1 hits the characteristic free-boundary at a point (τ, g δ (τ )). For the accurate Riemann solver, denote the reflected wave by β 3 . From Lemma 2.4, we have It follows that If the simplified Riemann solver is used, then β 3 is replaced by a non-physical front ǫ, and from (32), we have as well as Without loss of generality, we may assume C ′ 1 > 1 in (32), so this condition implies (44). Case 3. Interaction between a strong 3-rarefaction front and 1-or 2weak fronts from above.
In this case, without loss of generality, we consider a 3-strong rarefaction front s interacts with a 1-weak front α 1 . It is similar to analyze the other situation. Let the below and above states of the 3-strong rarefaction front s be U l and U m , respectively. The incoming 1-wave front α 1 connects the states U m and U r . Denote the outgoing waves by γ j (1 ≤ j ≤ 2), s ′ , and the non-physical front by ǫ, respectively.
The following lemma is standard [2, p.133, Lemma 7.2]. Based on this lemma, we get which implies that by triangle inequality. It is easy to see that The estimates of Q k (k = 1, 2, 4) are more complicated. Let ǫ ′ be any nonphysical front lying below the interaction point at "time" τ , and S ǫ ′ (τ ) the total strength of strong 3-rarefaction fronts lying above ǫ ′ . Then Here we assumed that which implies that |O(1)|K np |α 1 ||s| < 1, and used the fact that |e x − 1| ≤ 3|x| for |x| ≤ 1. It follows that Similarly, by considering all those weak fronts β 2 lying above α 1 at x = τ , we could obtain provided that We notice that the term W (α 1 , τ −) shall be retained to balance a term appeared in Q 1 below.

MIN DING AND HAIRONG YUAN
Now we turn to Q 1 . We firstly observe that by definition, Then Here in the second last inequality we used e −x − 1 < − 1 2 x for 1 > x > 0, and O(1)|s|e −Kωs ≤ C 1 |s| since s > 0. The weights will be chosen independent of δ, so we may choose δ small such that K ω s < 1. We remark that the computation carried here is one of the key point in dealing with large rarefaction waves.
Then considering any 1-wave front β 1 lying above α 1 at x = τ , like what we obtained before, one has provided that (47) holds. Hence we have Here we also used the assumption (57). Finally, from the above estimates of Q k , we get here for the last inequality, we used W (α 1 , τ −) ≥ 1, and the term with underline will be negative. Therefore We wish to choose K ω large so that which implies that Here we used the fact that C 1 ≥ 1.
If we consider a weak 2-wave front interacts with a 3-strong rarefaction front, we may get similar estimates as above, provided that Case 4. Interaction between a 3-strong rarefaction front and a 3-weak shock front from the below (above).
Suppose that a 3-strong rarefaction front s > 0 and a 3-weak shock front α 3 < 0 interact at a point on x = τ . Let γ k and ǫ be the outgoing k-waves (k = 1, 2, 3) and non-physical front, respectively. Depending on whether the produced 3-wave is a rarefaction wave or a shock, we consider the two subcases.
We do not know exactly the interaction potential of α 3 at x = τ − , but which is anyway nonnegative. Hence we have Now for k = 1, 2, let β k be any weak front lying above α 3 at x = τ −, and S ′ the total strength of strong 3-rarefaction fronts lying below β k at x = τ −. Then

MIN DING AND HAIRONG YUAN
Here we used the fact that e −x − 1 ≤ − 1 2 x for 0 < x < 1, and the assumption that It follows that Similarly, we have Summing up, one has if K 3 is sufficiently large so that Case 4.2. γ 3 < 0. For this case, and Direct calculation yields Here the summation is over all those non-physical fronts ǫ ′ lying below ǫ at x = τ . Suppose S ′ is the total strength of strong 3-rarefaction fronts lying above ǫ ′ . We used the fact that Similarly, for k = 1, 2, Then we conclude So if there follows Case 5. Interaction between a 3-strong rarefaction front and a nonphysical front from the below.
Suppose that a front s of the 3-strong rarefaction waves and a non-physical front ǫ interact when x = τ . Let s 0 and ǫ 0 be respectively the outgoing rarefaction wave front and non-physical front. Then we have the following standard lemma.

It is clear that
To calculate Q k for k = 1, 2, note that since s 0 = s, we always have W (β k , τ +) = W (β k , τ −), for any weak k-front β k lying above s at x = τ . Since no weak front from k-th family is involved in this case, we have Similarly, for any non-physical front ǫ ′ lying below ǫ, we have W (ǫ ′ , τ +) = W (ǫ ′ , τ −). This implies that Here we used W (ǫ 0 , τ +) ≥ 1, 1 − e x < −x for x > 0, and positiveness of ǫ, s. It follows that Therefore, if K np is large enough so that then Case 6. Interaction between a non-physical front and an i-weak front from the above (i = 1, 2, 3).
Suppose that a non-physical front ǫ interacts with an i-weak front α i (i = 1, 2, 3) from the above at some point on {x = τ }. Let the outgoing waves be γ i and ǫ ′ respectively. From the definition of the simplified Riemann solver, we have the following lemma (cf. [2, p.133 where O(1) is bounded, with a bound C 1 depending only on the background solution.
This completes the proof of Theorem 4.2.

Finiteness of fronts and interactions.
For the δ * determined in Theorem 4.2, and any δ < δ * , suppose (U δ , g δ ) is an approximate solution constructed by the front tracking algorithm. Let 0 < τ 1 < τ 2 < · · · < τ k < · · · be the interaction "time". Note that for x = 0, we solve all the Riemann problems by accurate Riemann solver. Particularly, by Lemma 2.3, we have Also recall that ε 0 by property of Riemann problems and the assumption in Theorem 1.2; and by definition, . All the constants C 3 appeared here depend only on the Euler system and the background solution, through the constants C 1 and weights chosen in previous sections. Hence we could choose ε 0 ≤ 1 claimed in Theorem 1.2 so that C 3 (ε 0 + ε 2 0 ) < δ * . By Theorem 4.2, we infer that for any approximate solution (U δ , g δ ) defined on 0 ≤ x ≤ τ .
Corollary 1. For given δ ∈ (0, δ * ), the number of fronts in the approximate solution (U δ , g δ ) at "time" x = τ is bounded by a constant independent of τ , and there are only a finite number of interactions of fronts.
Proof. Let N (τ ) be the number of fronts at x = τ . Then N 0 = N (0+) is finite and depends on the number of initial jumps in U δ (0, y) and ε 0 , δ.
The changes of N (τ ) can only occur in two situations.
2) The simplified Riemann solver is used. For Cases 2,4,5,6, the number of fronts does not change. For Cases 1, 3, only one new non-physical front is born. Therefore, we know the number of physical wave fronts is still N p = N p (δ * , δ, µ δ , N 0 ). Now, two physical fronts can only interact at most twice (possibly one before and one after reflection from the free-boundary). Hence the total number of new non-physical fronts is at most N p + 2N 2 p . Therefore we have Recall that an approximate solution U δ (x, y) is piecewise constant, with discontinuities across finite (say N ) fronts y = y k (x), k = 1, . . . , N , which are labeled so that y 1 (x) > y 2 (x) > · · · > y N (x) = g δ (x). Each y k connects the state U k−1 above it to the state U k below it. Let H(s) be the Heavside step function, whose value is zero for negative argument and one for positive argument. Then we could write Now set Θ k = 1, y k is a weak wave front or non-physical front, 0, y k is a 3-strong rarefaction front. We consists of only weak fronts, and consists of only 3-strong rarefaction fronts. Let p δ a (respectively p δ b ) be the pressure ahead of the upmost (respectively behind the lowermost) 3-strong rarefaction fronts in U δ . Then by decreasing of pressure across rarefaction waves (from above to below), and note that pressure increases only passing a 3-shock front in the middle of the rarefaction waves fans (cf. Lemma 2.2), we have . By triangle inequality and (70), Particularly this implies that TV.p δ (x, ·) is uniformly bounded. By Lemma 2.2, the total variation of U δ introduced by all strong rarefaction fronts is bounded by the total variation of p δ , while those introduced by 1, 2, 4-wave fronts could be controlled by L 0 (x), so TV.U δ (x, ·) is uniformly bounded (independent of x and δ).
Remark 6. As shown in [2, Section 7.5, p.146] or [9, p.530, Section 14.5], we can similarly prove that the total variations of U δ on space-like curves are uniformly bounded (independent of δ).

5.
Existence of Entropy Solutions. The results in the previous section guarantee that we could construct (U δ (x, y), g δ (x)) for all 0 ≤ x < +∞. Then by standard compactness arguments as shown in [5, Section IV, A], there is a subsequence {δ j } ∞ j=1 that converges to zero, and functions U, g so that g δj (x) converges uniformly to g(x) on any bounded interval, with the estimate (6) following from Remark 3, and U δj converges to U in C([0, T ]; L 1 ) after suitable shifts in y-variable. The estimates (7) and (8) then follow directly from (71)(73).
Therefore, to complete proof of Theorem 1.2, we need only to show that the limit (U, g) found above is actually an entropy solution to problem (2)(4). We follow the idea presented in [15, pp.299-305]. In fact, once (70) is established, there is little difference from our situation to the standard theory.
Let N P(τ ) be the set of all non-physical fronts lying on x = τ in the approximate solution (U δ , g δ ). The key point is to show the total strength of non-physical fronts (called ghost waves in [15]) is bounded by O(1)δ, rather than δ * proved before.
To this end, denote G m (τ ) to be the set of front that lies on x = τ , has generation order m, in the approximate solution (U δ , g δ ). Since for m ≥ 1, all fronts of order m are weak ones, we could use (6.32) in [15, p.301], (with n = 3 and T there replaced by δ * ,) to obtain the inequality here ♯(A) is the cardinal number of a set A. For T m (τ ) being the total strength of fronts in G m (τ ), Lemma 6.6 in [15, p.301] claims that (with T (t) there replaced by δ * ) and C 1 , C 2 are constants depending only on the background solution.

STABILITY OF TRANSONIC JETS 31
For a non-physical wave ǫ, denote its generation order to beǫ. Then T m (τ ) We require δ * small so that C 2 δ * < 1 as we done before. Now we choose k 0 large so that C 1 (C2δ * ) k 0 1−C2δ * ≤ δ/2. Then as δ and k 0 are fixed, we could choose µ δ small so that C 2 C 1 µ δ k0−1 m=1 3δ * δ 2m−1 ≤ δ/2. So by using the simplified Riemann solver judiciously, one could obtain that for the approximate solution (U δ , g δ ). Since by our construction, the free-boundary is accurate in each approximate solution, we could proceed in the same way as in [2,Section 7.4] or [15, pp. 304-305] to show consistency, namely the limit of (U δ , g δ ) must be an entropy solution to problem (2)(4).