Two Dimensional Riemann Problems for the Nonlinear Wave System: Rarefaction Wave Interactions

We analyze rarefaction wave interactions of self-similar transonic irrotational flow in gas dynamics for the two dimensional Riemann problems. We establish the existence result of the supersonic solution to the prototype nonlinear wave system for the sectorial Riemann data, and study the formation of the sonic boundary and the transonic shock. The transition from the sonic boundary to the shock boundary inherits at least two types of degeneracies (1) the system is sonic, and in addition (2) the angular derivative of the solution becomes zero where the sonic and shock boundaries meet.


Introduction
It is well known that the wave patterns created by initial value problems for the multidimensional compressible Euler system are complicated and challenging to study. As such, numerous works are focused on special initial data that allow us to reduce the space dimensions. In particular two-dimensional Riemann problems can be reduced to self-similar problems, which are interesting to study on their own. For instance, the flow in self-similar coordinates changes its type; it is supersonic (hyperbolic) in the far-field and becomes mixed near the origin. A seminal work by Zhang and Zheng [34] illustrated various complicated wave patterns on four sectorial Riemann data for two-dimensional self-similar Euler systems. Their conjectures are later validated numerically, see for example [18,19,24]. In particular Lax and Liu [19] attested the complicated wave patterns including Mach reflections, rolling up and instability of slip lines, and much more. They also noted that a general theory of Glimm type for one-dimensional problems would not be likely established for the multidimensional problems.
While establishing a comprehensive analysis to understand the complicated wave patterns in two dimensional flows with Riemann data still has a long way to go, there is recent progress describing the solution structures of transonic regular shock reflection problems for various systems, [2,3,4,5,6,7,12,14,15,16,17,22,26,27,35]. This is an incomplete list of related work specifically on transonic shocks by Riemann problems or by a wedge, Date: November 6, 2018. The work of Kim was supported by the National Science Foundation under the Grants DMS-1109202, 1615266.
The work of Tsikkou was supported by the National Science Foundation under the Grant DMS-1400168. 1 and interested readers can refer to the references therein. On the other hand, relatively little is known about transonic rarefaction waves. The computational result by Glimm et al [13] showed the shock formations from the rarefaction wave interactions. This paper addresses the understanding of rarefaction wave interactions and their shock formations in a specific self-similar problem, the nonlinear wave system in two space dimensions. More precisely, the simple wave created by the planar rarefaction wave becomes a different wave with two families of non-trivial characteristics. Both families of characteristics merge into becoming sonic near the origin and the type of the system then changes to subsonic, that is the characteristics no longer transport the data any further. We call this wave a transient wave. On the other hand, there exists a family emanating from this transient wave region, becoming compressive and forming a shock downstream. Thus the type changes and forms a sonic boundary in some part and a transonic shock in the other.
We note that a similar configuration was studied by [30] for the pressure gradient system. They [30,36] called the region -where a family of characteristics starts on sonic curves and ends on transonic shock curves -the "semi-hyperbolic" region. [30] established the existence of a local solution in a given semi-hyperbolic region, provided a smooth convex boundary and small Riemann data, and [31] established the local regularity result for this solution. Our work is motivated by the work of [1,11,30,31]. The nonlinear wave system, which can be considered as wave motions of shallow water and multidimensional p-systems, is a reduced system of the compressible Euler system for isentropic, irrotational flow in two space dimensions [4,5]. The nonlinear wave system can be also considered as a part of an operator splitting scheme in numerics, where the compressible Euler system can be split into the nonlinear wave system (the pressure system) and the pressure-less system (the gradient flow). In fact [36] noted that the Euler system can be split into the pressure-gradient system and the pressure-less system, see [29,35] and the references therein. The pressure-gradient system is a special case of the nonlinear wave system. The pressure-less system is well understood by [28]. Hence if one understands the solution structure of the nonlinear wave system, then one can construct the solution of the Euler system successively by using the splitting method. Furthermore, there are many similarities on the structures of both the nonlinear wave system and the Euler system, see [4,25]. As such, it is crucial to understand the nonlinear wave system in order to study the Euler system.
We focus on the wave patterns created by planar rarefaction waves. For the configuration, we impose four sectorial Riemann data U i , i = 1, 2, 3, 4 for each ith quadrant, see the left figure in Figure 1. We consider planar waves with a horizontal rarefaction wave R 14 created by the constant data U 1 and U 4 , and another vertical rarefaction wave R 34 created by U 3 and U 4 . The contact discontinuities reside along the positive y−axis created by U 1 and U 2 , and along the negative x−axis by U 2 and U 3 , which have no effect on the system (the reduction of the system makes the contact discontinuities trivial), see the right figure in Figure 1. The data are symmetric with respect to x = −y and thus it suffices to focus only on R 14 . Near the locus of the sonic circle (it is the origin for the nonlinear wave system) the type of the flow changes and becomes subsonic, creating a transonic PSfrag replacements shock downstream. We formulate the boundary value problem and establish the supersonic solution for the entire hyperbolic region of this configuration.
Since the change of the type is not known a priori, the problem has two different types of free boundaries: the sonic boundary and the transonic shock boundary. We show that the sonic boundary in this configuration inherits at least two types of degeneracies.
The first obvious one is that the wave across the sonic boundary becomes degenerate meaning it is neither hyperbolic nor elliptic. The degeneracy of this type is well known as a Tricomi type problem, where the characteristics enter the sonic boundary perpendicularly. The Tricomi type degeneracy appears, for example, in the airflow over a wing where the steady subsonic flow creates a supersonic region over the convex surface of the wing creating a shock, see Courant and Friedrichs [9]. This is a long standing open problem. Numerical studies [32,33] suggest that similar wave patterns are observed for the Mach reflection problem near the Mach stem. We establish the existence of the supersonic solution, which becomes sonic, and that the solution and the sonic boundary are C 1 .
The second one is when a family of characteristics creates a compression wave downstream, the sonic boundary and the transonic shock merge into a point, and at that point, the angular derivative of the solution also disappears. To our knowledge, this is a new type of free boundary problem. We note that the Tricomi type degeneracy, since the characteristics enter the sonic boundary perpendicularly, implies that the sonic boundary is never the characteristics. Thus we utilize both the directional derivatives along the sonic boundary and the data transported along the characteristics to have the solution and the sonic boundary to be C 1 . The solution may not be C 1 at the point where the sonic and shock boundaries meet, and at that point, the sonic boundary no longer has a Tricomi type degeneracy but the angular derivative of the solution becomes zero. Our results provide an insight to understand how the compressive wave is created by the expansion wave.
For multidimensional conservation laws, the entropy conditions are insufficient to answer whether our solution is the physically relevant one. However, our solution captures the numerics, see Section 8. This paper provides a framework to establish existence of the supersonic solution suggested by the numerics, and an analysis to understand how the type of the rarefaction wave changes. The complete analysis to construct the transonic wave in the entire region including the subsonic region will be discussed in our forthcoming paper.
The main contributions of this paper are the following. We first formulate the boundary value problem for the self-similar nonlinear wave system. We next discuss the wave patterns, monotonicity properties, regularity and existence results. The solution will be constructed locally and then assembling the pieces together along the characteristics and the sonic boundary. We show that the sonic boundary created by the transient waves is C 1 and is strictly increasing radially. This sonic boundary is terminated and radially tangential when it merges to the transonic shock downstream. Numerical results by using CLAWPACK for certain pressures are presented as well.
We believe our results will serve as a vehicle for understanding transonic flows in particular the long standing open problem of the flow over the convex wing, and lead to further developments of systematic theories for multi-dimensional conservation laws. Interested readers can refer to the survey paper [8] for the comprehensive references and recent progress in transonic problems.
2. Description of the problem 2.1. Nonlinear wave system: Configuration. From the compressible Euler system for isentropic flow in two space dimensions, ignoring the nonlinear velocity terms (assuming low velocities) and assuming irrotational flow, we can deduce a simpler system, the nonlinear wave system [4]; (1) Here ρ(t, x, y) is the density, u(t, x, y) and v(t, x, y) are the x and y components of velocity, respectively, and p(ρ) is the pressure satisfying a polytropic gas law with constants k (we let k = 1 for simplicity), 1 < γ < ∞, (typically 1 < γ < 2, γ = 5/3 is air), and a local sound speed c 2 (ρ). This system can be considered as wave motions of shallow water and multidimensional p-systems. We let the momentum (ρu, ρv) = (m, n) and use U to denote (ρ, m, n): U = (ρ, m, n) = (ρ, ρu, ρv).
Specifically we have the following Riemann data U i = (ρ i , m i , n i ) where i = 1, 2, 3, 4 for each quadrant satisfying; That is, Hence we impose four sectorial data; We let c 1 = c(ρ 1 ) and c 4 = c(ρ 4 ) and denote two sonic circles We replace c 2 (ρ) = γp κ , κ = (γ − 1)/γ, and write the system (1) in self-similar coordinates ξ = x/t, η = y/t; The system can be written in a second order equation (by applying ∂ ξ to the second equation and ∂ η to the third equation in (2) and replacing the momentum terms with their derivatives from the first equation in (2)) which can be written in polar coordinates r 2 = ξ 2 + η 2 , θ = tan η/ξ: The system is hyperbolic (supersonic) when c 2 < r 2 , sonic when c 2 = r 2 , and elliptic (subsonic) when c 2 > r 2 .
In the following section we discuss the characteristic equations for the system in the supersonic region.

2.2.
Characteristic equations in the supersonic region. When the state is hyperbolic, that is c 2 < r 2 , we let λ ± = ±r r 2 − c 2 c 2 , or simply λ = r r 2 − c 2 c 2 , and the positive and negative characteristics derived by integrating dr dθ = λ ± = ±λ, respectively, so that (4) reads In addition, by letting ∂ ± = ∂ θ ± λ∂ r , we have, as in [30], We denote and which can be written as  In the following section, we discuss and formulate the boundary value problem for the system.

Derivations of the boundary value problem and the main result
When the rarefaction wave R 14 : c 2 (p) = η 2 , where p 4 ≤ η ≤ p 1 , enters the sonic circle C 1 , a new wave pattern is created. Specifically, the wave pattern changes from the simple wave to two wave families of non-constant solutions, and these two families of characteristics merge and become sonic near the positive η-axis. We denote the sonic boundary created by these two families of characteristics by σ. In addition, there exists a wave family connecting the sonic boundary and the transonic shock boundary. The transient wave region (consists of two wave families of non-constant solutions) is enclosed PSfrag replacements Figure 2. Configuration with details. by two limiting characteristics and the sonic boundary. More precisely, let Γ 12 be the positive characteristic emanating from Ξ 1 = (r 1 , π/2) to Ξ 2 = (r 2 , θ 2 ), which separates the simple wave R 0 = {c 2 (p) = η 2 } and the transient wave. Γ 12 is completely determined by the simple wave R 0 . Ξ 2 = (r 2 , θ 2 ) is the point at which this characteristic crosses η = c 4 , where r 2 = c 4 / sin θ 2 and θ 2 = sin −1 c 4 /c 1 . Let Γ 23 be the negative characteristic emanating from Ξ 2 to Ξ 3 where Ξ 3 = (r 3 , θ 3 ) is the point at which the characteristic speed becomes zero, that is, a sonic point. The transient wave region is the hyperbolic region enclosed by characteristics Γ 12 and Γ 23 , and the sonic boundary σ.
On the other hand, there exists yet another new simple wave created downstream, which is adjacent to the constant state p = p 4 . More precisely, when the governing hyperbolic system is reducible to a first order homogeneous system, the new state adjacent to the constant state forms a simple wave, see Courant and Friedrichs [9]. We refer to Dafermos [10] for details of the general framework of the simple wave. This simple wave becomes compressive and creates a transonic shock downstream.
In summary, the transonic boundary is created by the transient wave (the sonic boundary) in part, and the simple wave (the shock boundary) in the other. Hence it is crucial to identify the corresponding characteristics first to separate these different wave regions, in order to formulate the correct boundary problems. Below we list the notation used for the characteristics and the corner points that we just discussed. Let The supersonic region is divided into three regions based on the characteristics where R = ∂ + p and S = ∂ − p are either strictly positive or zero, and denoted by; • Rarefaction wave region R 0 : R > 0 and S = 0, see Section 4.
• Transient wave region R 1 : Both families R and S are non-trivial, enclosed by the positive characteristic Γ 12 emanating from Ξ 1 to Ξ 2 ; the negative characteristic Γ 23 from Ξ 2 to Ξ 3 ; and the sonic boundary σ on which R = S. See Sections 5 and 6. • Simple wave region R 2 : R = 0 and S > 0, see Section 7. The boundary of this region consists of the negative characteristic Γ 23 emanating from Ξ 2 to Ξ 3 ; the positive characteristic, denoted by ; and a part of the shock Σ ′ ⊂ Σ from Ξ 3 to Ξ 4 . In region R 1 , we have the following Goursat boundary value problems: where We discuss the compatibility conditions in Section 5.
In R 2 , the positive characteristics, emanating from Γ 23 , form an envelope, where the positive characteristics satisfy for each (r 23 (θ i ), θ i ) ∈ Γ 23 , θ 2 ≤ θ i ≤ θ 3 , and r = r 23 (θ) satisfies (11) dr Remark 3.1. We note that the simple wave does not create a sonic curve. More precisely if it is sonic somewhere, that is λ = 0, in the simple wave region. Then S = ∂ − p = p θ −λp r and R = ∂ + p = p θ + λp r imply p θ = S = R, which yields a contradiction: simple wave.
The family S > 0 remains to be non-trivial and this is the one that carries the data from σ to Γ 23 between the two regions R 1 and R 2 . The other family R becomes zero, and therefore carries no information in this simple wave region R 2 . Thus we have We can write the solution to this Riccati type equation in the following form: Let (θ,r) a point on the negative characteristic Γ 23 and integrate (12) along the positive characteristic lines in region R 2 , then We finally state the main theorem.
. For a convex Γ 23 ∈ C 2 and the data (p, R, S) ∈ C 2 (Γ 23 ) satisfying the compatibility conditions at Ξ 2 and Ξ 3 and on Γ 23 , there exists a supersonic solution (p, R, S) the Goursat boundary problems (7)- (9). The solution (p, R, S) creates the sonic boundary The sonic boundary σ and the transonic shock Σ ′ merge into a point Ξ 3 at which the solution (p, R, S) holds Furthermore there exists a simple wave creating a transonic shock Σ ′ in region R 2 .
Our existence result of the supersonic flow is established with the given convex negative characteristics Γ 23 and the data on Γ 23 holding the compatibility conditions. In our forthcoming paper, we establish the global transonic solution and provide the scheme to select the correct data and Γ 23 .
In what follows, we use the condition 2c(p 4 ) > c(p 1 ) and discuss the existence results for each region.

Transient wave Region R 1
We first discuss many useful properties of the characteristics in the transient wave region R 1 , in the same spirit as in [1]. More precisely we discuss the monotonicity and convexity properties of the characteristics, and the monotonicity of p along the characteristics in polar coordinates and cartesian coordinates (different coordinates provide different aspects of the characteristics) in Lemmas 5.1 -5.4. We also state a priori bounds at the end of the section, in Lemmas 5.5, 5.6.
Lemma 5.1. The hyperbolic solution p ∈ C 1 to the Goursat problem satisfies Proof. Let (θ,r) be a point on the positive characteristic Γ 12 , and (θ,r) a point on the negative characteristic Γ 23 . Integrate (8) and (9) along the negative and positive characteristic lines, respectively, to obtain The data (16) ensures R > 0 only. Thus we check the positiveness of S to complete the proof. From (19), if S = 0 at (θ,r) then S = 0 along the positive characteristic passing through (θ,r) in region R 1 . Thus R = S = 0 at a point different from Ξ 1 on σ which is a contradiction. Therefore S = 0 along Γ 23 , possibly excluding the endpoints. If S < 0 at (θ,r) then R = S < 0 somewhere on the sonic curve which is again a contradiction to R > 0. So we conclude that S > 0 along Γ 23 for θ = θ 2 , θ 3 . Thus R, S > 0 and consequently p θ > 0 in the interior of region R 1 .
By Lemma 5.1, we deduce R = 0 in the simple wave region R 2 . We next discuss monotonicity properties of the characteristics. To ease the analysis, we write the characteristics in self-similar coordinates. From dr dθ = ±λ, the characteristics Let the corresponding directional derivatives in the self-similar coordinates be; We observe that in region and (c 2 − ξ 2 )Λ 2 ± + 2ξηΛ ± + c 2 − η 2 = 0, which can be solved for c 2 to get or We next discuss properties of the characteristics.
Lemma 5.2. The hyperbolic solution p ∈ C 2 (R 1 ) to the Goursat problem has the following properties: (1) Along the negative characteristics dη dξ = Λ − starting from any point on Γ 12 \ Ξ 1 : (2) Along the positive characteristics dη dξ = Λ + starting from any point on Γ 23 : and by Lemma 5.1, we obtain the strict inequality thus using (21), we conclude that , which means that the negative characteristics are convex.
(iii)-(iv) Evaluate the first equation of (25) on This immediately implies dp and Λ − = − ξ η = 0, respectively, at the points on the η-axis. On the other hand, If not then we have unbounded dp d − η when Λ − = 0 at some point (r * , θ * ) on the negative characteristics. However c 2 = (η * ) 2 when Λ − = 0, and at the same time dp d − ξ = p ξ < 0. At the points on Γ 12 we know that c 2 = η 2 so by (20) the negative characteristics satisfy Λ − = 0. Hence a convex characteristic has at least two points with Λ − = 0 which leads to a contradiction. Thus, which is again a contradiction to the convexity and the behavior of the negative characteristics along Γ 12 . We therefore conclude that the change of type occurs in the first quadrant, R 1 is located in the first quadrant and Λ − < 0, dp (v) By Lemma 5.1, and we conclude that dp d + η > 0 everywhere in the interior of R 1 .
(vi) In addition differentiating (24) gives By (22), we conclude that which means that the positive characteristics are concave.
Thus by a similar argument as before we conclude the claim. More precisely, from 0 < dp we show that Λ −1 + = 0 to obtain Λ −1 + < 0 and dp d + ξ < 0 in the interior of R 1 by a contradiction argument. Suppose not. Then c 2 = ξ 2 and dp d + η = p η > 0 at the contradiction point (r * , θ * ) on the positive characteristics. On the other hand for the negative characteristics passing from (r * , θ * ) the following hold: However we have shown that dp d − η > 0.
In the following two lemmas, we discuss the properties of the sonic boundary σ; in particular the monotonicity and the corner point Ξ 3 with the level curve where {p = p(Ξ 3 )}. In addition R = S on the sonic boundary. In particular R = S > 0 on σ \ {Ξ 1 , Ξ 3 }.
If p η = 0 at Ξ 3 then p ξ = 0 and thus by (20) we have If p η = 0 at Ξ 3 then We thus conclude that a level curve meets tangentially the sonic boundary at Ξ 3 .
Proof. By Lemma 5.1 we first note that in region R 2 we have R = 0 and S > 0 along the positive characteristics and therefore the sonic boundary along which R = S cannot extend below Ξ 3 .
Assume now that there is a point (ξ * , η * ) on σ such that c 2 = (c * ) 2 and dη dξ = 0 at this contradiction point on the sonic boundary. Then from the unboundedness of dξ dη , we have d 2 ξ dη 2 < 0 when ξ > ξ * and In the latter case, let Closely related to the above, in spirit as well as in technique, we end up with a contradiction.
In (32), notice that if (c 2 ) ′ p η = 2η then (c 2 ) ′ p ξ = 2ξ and thus R + S = 0 which does not hold for points along the sonic boundary different from Ξ 1 and Ξ 3 . We therefore conclude that (c 2 ) ′ p η = 2η and does not change sign along σ excluding the endpoints, thus Let us assume that (c 2 ) ′ p η > 2η everywhere along σ, then p η is positive, and we also know from Lemma 5.3 that the sonic boundary satisfies dξ dη < 0, in the neighborhood of Ξ 3 . On the other hand, in the neighborhood of Ξ 1 , since c 2 = η 2 in region R 0 and c 2 = r 2 on σ, we expect that the level curves have positive slopes and thus p ξ < 0 by (33). This combined with (34) gives dξ dη > 0 for the sonic boundary, in the neighborhood of Ξ 1 . This implies that there must be a point along σ, different from the endpoints, such that dξ dη = 0. which is a contradiction. We therefore conclude by (30) that (c 2 ) ′ p η < 2η and (c 2 ) ′ p ξ < 2ξ along σ.  Lemma 5.5. The solution p ∈ C 1 (R 1 ) satisfies Proof. By Lemma 5.1, we have Thus p is strictly increasing in the θ direction, and p > p(Ξ 2 ) = p 4 in R 1 . Note that the change of type occurs in the first quadrant as the characteristics enter the sonic circle, where c 2 (p) = r 2 , in the first quadrant, and becomes subsonic holding c 2 (p) > r 2 . Thus when θ ≥ π/2 we are now in the subsonic region. Therefore p < p(Ξ 1 ) = p 1 in R 1 .
Similarly, these bounds hold in the simple wave region R 2 . While it is not straightforward to see whether these bounds remain valid in the subsonic region, we note that this lemma holds in the entire domain and refer to the forthcoming paper on the transonic problem.
We next cite estimates from [30].
Lemma 5.6. The maximum values of ∂ ± p of the solution p ∈ C 2 are attained on the characteristic boundaries Γ 12 ∪ Γ 23 . Furthermore the solutions (p, R, S) ∈ C 1 satisfy where t = r 2 − c 2 (p).

The existence result in the transient wave region R 1
We formulate the Goursat boundary problem to construct the transient wave in R 1 . We first discuss the boundary data on the negative characteristic Γ 23 . Let f ∈ C 2 (θ 2 , π/2) be convex, and and Additionally, we require that f and g satisfy the following compatibility condition: [G1.] There exists θ 2 < θ 3 < π/2 such that The compatibility condition [G1.] implies We note that due to the conditions that df /dθ < 0, c 2 (g) ≤ f 2 while g θ > 0 (since ∂ − g > 0 and R[f, g] > 0), we have immediate bounds for f and g We construct a solution that creates a sonic boundary where R = S becomes zero as the solution approaches Ξ 3 on the sonic boundary.
The second condition may be rewritten in the form Hence the Goursat problem has a local solution near Ξ 2 . Next we establish a local solution to an initial boundary value problem. Let the initial position be I = {(r(θ), θ)} which is to be determined, and W I = W (r(θ), θ). Let X = I ∩Γ 12 and Y = I ∩ Γ 23 be the points where the initial and boundary positions meet. Now the compatibility conditions are as follows: for i = 0, −, we have and for i = 0, +, we have The initial position I ∈ C 1 is chosen to be the constant level set of r 2 − c 2 (p) so that r ′ = λ ± , which allows us to match the compatibility conditions. Thus we have the local existence result.
In order to establish the global solution to the entire region R 1 , which is enclosed by Γ 12 , Γ 23 and σ, we first discuss the regularity results near the sonic boundary due to [31]. We note however that the regularity result is limited to strictly positive R and S. Hence the estimates depend on δ 0 where R, S ≥ δ 0 > 0. Equipped with these regularity results near the sonic boundary, noting that p θ > 0, and the characteristics entering the sonic boundary in the radial direction and never parallel to the tangential direction of the sonic boundary, the sonic boundary σ is estimated by the gradient of the pressure which then leads to the existence of a C 1 solution and σ ∈ C 1 . 6.2. Regularity results. We discuss the Hölder gradient estimates near the sonic boundary. The estimates are established by [31] for the pressure gradient system provided that ∂ ± p are strictly positive, and under certain smoothness of the solution and the sonic boundary. While our result relies on [31], we provide insights which signify the estimates and their consequences.
We first change the coordinate system to flatten the sonic boundary as it was done in [31]. This new coordinate system brings a couple of technical advantages. The first obvious one is that it simplifies the geometry of the sonic boundary. The next one is not immediate: the new coordinates enable us to derive the corresponding system that provides estimates on t β |R r |, t β |S r |, where t = r 2 − c 2 (p), for sufficiently small 0 ≤ t ≤ t 0 and uniformly bounded R, S such that R, S ≥ δ > 0, where 1 < β = β(t 0 , δ) < 2. This is essential to establish the sonic boundary to be in C 1 .
Since the estimates depend on the strict positiveness of R, S, we consider the region R 1 excluding small neighborhoods of Ξ 1 and Ξ 3 , where δ is the distance from these corner points, see Figure 4. Let R 1 [δ] be neighborhoods of Ξ 1 and Ξ 3 with δ > 0 distance away from these points where Γ δ − is the negative characteristic with δ = dist(Ξ 1 , Γ δ − ), and Γ δ + is the positive characteristic with δ = dist(Ξ 3 , Γ δ + ). The Hölder gradient estimates are established in the region R 1 \R 1 [δ], in particular near the sonic boundary. Let p be smooth enough to derive the characteristic equations for the derivatives (for simplicity we may assume p ∈ C 3 (R 1 ) for now, which can be weaken). Let t = √ r 2 − c 2 , and consider the new coordinate system (r, t). From simple calculations, PSfrag replacements , the region R 1 excluding the small neighborhoods of Ξ 1 and Ξ 3 .
the system in (r, θ) coordinates becomes .
The corresponding characteristics read and the last two equations can be rewritten as As discussed in [31], the following equations will be useful to establish the estimates. Let and obtain and Evaluate and write dG where wherẽ We now establish, in the same spirit as in [31], the following regularity result. Theorem 6.2. For given Riemann data 2c(p 4 ) > c(p 1 ), t 0 > 0 sufficiently small, and δ > 0, there exist a positive constant C and β ∈ (1, 2) depending only on t 0 , δ, the Riemann data and max Γ 12 ∪Γ 23 (R, S), such that the solutions R, Moreover, R, S and t −1 (R − S) are uniformly continuous in Proof. Recall that Hence in order to have the sonic line to be in C 1 , we show that (R − S)/t is uniformly bounded, and R, S and (R − S)/t are uniformly continuous in R 1 \ R 1 [δ], For each fixed t = t b > 0, we consider the level curve {r 2 − c 2 (p(r, θ)) = t 2 b } where θ = θ(r). We then have Since R and S are positive and bounded in region R 1 , we have p θ = (R + S)/2 positive and bounded. Thus θ ′ (r) is well-defined on each level curve.
Recall that V = 1/S − 1/R and V satisfies (57). We have m i , i = 0, 1 positive constants such that |(µ 0 − 1)/t| ≤ m 0 δ −1 , and |µ j | ≤ m 1 δ −3 where j = 1, 2, in region R 1 \ R 1 [δ]. Hence we can write equation (57) in the form Integrating the last equation from t to t b , we have We then deduce Observe that for t 0 > 0, there exist positive constants L 0 and F 0 such that |l(r, t)| ≤ L 0 δ −1 , Hence for t ≤ t b ≤ t 0 , we get from the last inequality We then deduce Thus for 1 < β < 2, by choosing t b ≤ t 0 sufficiently small if necessary, we then have Hence we now have established the bound of M ; This uniform bound (similarly to H) holds for any 1 < β < 2, and immediately gives Now with this uniform bound, from inequality (64) for V /t, noting |e| ≤ E 0 δ −1 , we have Hence we get We next show that R, S and V /t are uniformly continuous. We now have for some constant C > 0 uniformly in t. Hence integrate the last inequalities along the negative and positive characteristics respectively to obtain Since R = S along the sonic line and both R and S are continuous inside region R 1 , we have as |r 1 − r 2 | → 0. Thus R and S are continuous on the sonic line and uniformly continuous in Next, integrating (57) from 0 to t b where t b is arbitrary chosen, we have which then becomes Hence with Therefore we have established the claim.
6.3. The supersonic solution in region R 1 . The existence result in the entire region R 1 is established in two steps. We first show Lemma 6.3 to establish the sonic boundary σ where R = S. We next show that σ ∩ Γ 12 = Ξ 1 and σ ∩ Γ 23 = Ξ 3 . The proof of lemma 6.3 is inspired by the work in [1,11]. [11] and later [1], established the global existence result for the degenerate hyperbolic system of the pressure gradient equation, where the pressure becomes zero at the origin which makes the system degenerate only at the origin and hyperbolic elsewhere. Furthermore since the pressure gradient system is quasilinear, the type of the system (whether supersonic or subsonic) must be also identified. While our system may appear to have similar technical difficulties as in [1,11], we note that it is not straightforward to apply their result to our case. We also note that the sonic boundary will be determined by the choice of the data on Γ 23 and thus the compatibility conditions at Ξ 1 and Ξ 3 will play crucial roles in selecting the correct data on Γ 23 to find the supersonic solution in the entire region R 1 .
Let u = r 2 − c 2 (p) so that u = 0 on the sonic boundary σ. The characteristic equations in u become Thus we have , there exists a level curve l τ = {(τ (θ), θ) : u(τ (θ), θ) = d} ⊂ R 1 and τ = τ (θ) ∈ C 1 such that the following hold: (1) there exists U = (u, R, S) ∈ C 1 (D τ ) satisfying the Goursat boundary problem Proof. Let L be the set where d ∈ L satisfying the assertions (1)- (2) in this lemma. Since the system for W = (p, R, S) is equivalent to the new system for U = (u, R, S), the local existence result from the system of W = (p, R, S) ensures the existence of d 0 ∈ L (that is L = ∅) and [d 0 , d m ] ∈ L (since u θ < 0). Hence we only need to show that inf L = 0 to establish the claim. We show the claim by contradiction. Suppose that inf L = d * > 0. The proof consists of two main steps. By extracting a limit to first show that d * ∈ L, and then show that this d * violates the infimum assertion.
Since d * is assumed to be the infimum of L, there exists a monotone decreasing sequence {d n } ⊂ L satisfying lim n→∞ d n = d * . Then for each d n , there exists τ n = τ n (θ) which satisfies the assertions (1)-(2) where l n is enclosed in the domain bounded by Γ 12 and Γ 23 , that is l n = {(τ n (θ), θ), θ n 23 ≤ θ ≤ θ n 12 } where l n ∩ Γ 12 = (r n 12 , θ n 12 ) and l n ∩ Γ 23 = (r n 23 , θ n 23 ). We let D n be the domain enclosed by l n , Γ 12 and Γ 23 . The uniqueness of the local existence results ensure the monotonicity of D n such that D n ⊂ D n+1 . Thus we let l * be the graph of τ * (θ) where lim n→∞ τ n (θ) = τ * (θ), for all θ * 23 ≤ θ ≤ θ * 12 . Let D * be the closed region enclosed by l * , Γ 12 and Γ 23 .
Hence in D * \ l * , there exists the solution U = (u, R, S) satisfying (due to Lemmas 5.1-5.6 and Theorems 6.1, 6.2) where c 0 , c 1 > 0 depends only on p 1 , p 4 and d * . Since the tangential derivative of u along l τ is we now obtain τ C 1,β ≤ C. By using the standard compactness argument we have τ C 1,α ≤ C 1 , for any 0 < α < β. Repeating a similar compactness argument we have U = (u, R, S) ∈ C 1,α with the uniform bound established from the regularity results. Thus we extend U in D * to satisfy the governing system. Now observe that l * is the constant level set of u, and thus On the other hand (75) u θ + λ ± u r = 0 on l * .
Next we show that there exists ε > 0 small such that d * − ε ∈ L to contradict inf L = d * . Since we have shown that d * ∈ L, by the local existence result, we extend the solution in C 1 by solving the initial boundary value problems where the initial value on l * satisfies the compatibility conditions on Γ 23 ∩ l * and Γ 12 ∩ l * . Furthermore the solution is unique (since the data is prescribed uniquely) and u θ < 0, the corresponding solution determines the level curve where u = d * − ε > 0 for some small ε > 0 satisfying the assertions (1)- (2) in this lemma. Hence repeating the same argument as to d * we can show d * − ε ∈ L, which is a contradiction.
Therefore inf L = 0 and this completes the proof. Lemma 6.3 implies the existence of the sonic boundary σ where R = S. Due to the monotonicity properties of the characteristics and the solution, the sonic boundary σ is connected, bounded and will not form a closed loop. We now check whether this sonic boundary satisfies the compatibility conditions at Ξ 1 and Ξ 3 .
We first consider Ξ 3 . If σ meets at X 0 ∈ Γ 23 where X 0 = Ξ 3 . This then violates the data at X 0 since S = R on Γ 23 . So if σ meets somewhere on Γ 23 it must be only at Ξ 3 . Now if σ does not intersect with Γ 23 , that is the sonic boundary is terminated before it reaches to Γ 23 and dist(σ, Γ 23 ) = δ > 0. In that case, since R and S must be positive in the region R 1 , we have R = S > 0 on σ. In particular we can find the positive characteristic Γ δ + connecting the boundary point, denoted by X 1 , of σ, and Γ 23 . Hence we can treat this as a new boundary value problem where the boundaries consist of Γ δ + and a segment, denoted by Γ δ 23 , of Γ 23 from Ξ 3 to Y 1 = Γ δ + ∩ Γ 23 . We note that the data on Γ δ + is now prescribed with the given data on Γ 23 \ Γ δ 23 . Note that R = S > 0 on σ. Thus we may write R(X 1 ) = S(X 1 ) = q 0 > 0 and τ ′ ≥ ε 0 > 0 for some constants q 0 , ε 0 (however note that these constants q 0 , ε 0 depend on the choice of the data g on Γ 23 ). Hence we find 0 < δ 1 < δ 0 and ε 1 such that we can extend τ ′ (θ) ≥ ε 1 > 0, where ε 1 < ε 0 , in a small neighborhood of X 1 , denote B δ 1 (X 1 ) to be the circle centered at X 1 with radius δ 1 . We now find the positive characteristics Γ δ 2 + connecting Y 2 ∈ Γ 23 and ∂B δ 1 /2 (X 1 ), see Figure 5. If there exists δ 1 such that B δ 1 (X 1 ) covers the entire neighborhood of Ξ 3 then our choice of g on the data to prescribe Γ 23 needs to be adjusted (scale down the strength on S).
Repeat this process by finding ε n → 0 to construct σ[δ n ] for each ε n > 0 and find sequences {X n }, {Y n }, and q n , where R = S ≥ q n on σ[δ n ]. By the construction (these sequences are monotone and bounded) we know that there exist limits for these sequences. Let X * , Y * , q * be the limits of these sequences respectively. Note that the tangential derivative on r 2 − c 2 (p) along σ is (c 2 ) ′ (p θ + τ ′ p r ) = 2rτ ′ which implies R = S = p θ on σ and p θ = q n → q * = 0, as ε n → 0.
Hence we are left to check whether X * = Y * . Suppose X * = Y * . We have two possibilities either (1) X * surpasses below Γ 23 , or (2) the sequence {X n } is terminated before it reaches to Γ 23 . (We may treat this as a shooting method, that is, case (1) is overshoot, and case (2) is undershoot). We note however by the construction (we have adjusted the data g so that δ i can be chosen appropriately), X * should not be located below Γ 23 . Hence the sonic boundary constructed by this sequence is then terminated before it reaches Γ 23 . In other words we need to readjust the choice of the data g. Therefore there exists data g such that the limits of X n and Y n match, and denote the limit by Ξ 3 .
For the case Ξ 1 . As before, since R > S = 0 on Γ 12 , if σ meets somewhere on Γ 12 , it must be at Ξ 1 . Now if σ does not intersect with Γ 12 , as we did for Ξ 3 , we formulate a boundary value problem where the boundaries are now with the corresponding characteristics Γ δ − and Γ δ 12 . Repeating a similar argument as we did for Ξ 3 and noting that the data on Γ 12 is given (independent of the data on Γ 23 ) while the data on Γ δ − depends on the data on Γ 23 we find the correct data so that it matches the compatibility condition at Ξ 1 .
Therefore we finally establish the existence result in the entire transient wave region R 1 . Lemma 6.4 provides the schematics, as illustrated in Figure 5, that is, how to construct PSfrag replacements Figure 5. Schematics of constructing the solution near Ξ 3 .
the solution near Ξ 1 and Ξ 3 by selecting the correct data on Γ 23 to hold the compatibility conditions at Ξ 1 and Ξ 3 . Furthermore Lemmas 6.3, 6.4 utilize the properties of the tricomi type degeneracy, that is, the sonic boundary cannot be a characteristic curve, and as a consequence R = S = p θ remain strictly positive. Hence the natural compatibility condition where the different types of boundaries meet is R = S = 0 which gives rise to an additional degeneracy.
We further find the positive characteristic emanating from Ξ 2 , denoted by Γ 24 , by integrating  PSfrag replacements Γ 12 Figure 6. Envelope formation by the simple wave.
These results are produced by using CLAWPACK [21]. We implement Roe average methods [23] and finite volume methods on quadrilateral grids [21]. More precisely, we implement Roe average methods in a uniform grid in polar coordinates as our computational domain, together with a coordinate mapping and appropriate scaling of the flux differences. The scaling is done by using the area ratio "capacity" of the computational cell which is determined by the size of the corresponding physical cell [20].
In Figure 8: the right figure is the enlargement from the left figure near which the shock appears. In the right figure, the density flattens out near the shock, while there exists a compression (a dip in the cross section) which merges to the shock. The numeric suggests that the angle of the location of Ξ 3 where the sonic boundary and the shock boundary meet is between 50 and 60 degrees in this configuration.