Entropy-based moment closure for kinetic equations: Riemann problem and invariant regions

Jean-François Coulombel and Thierry Goudon

CNRS & Université Lille 1, Laboratoire Paul Painlevé, UMR CNRS 8524 Cité scientifique, 59655 VILLENEUVE D’ASCQ Cedex, France and Team SIMPAF, INRIA Futurs

        E-mails: jfcoulom@math.univ-lille1.fr, thierry.goudon@math.univ-lille1.fr

March 1, 2006

Abstract

We study a nonlinear hyperbolic system of balance laws that arises from an entropybased moment closure of a kinetic equation. We show that the corresponding homogeneous Riemann problem can be solved without smallness assumption, and we exhibit invariant regions.

AMS subject classification: 35L60, 35L67 82C40,

1 Introduction

This paper is devoted to the analysis of the following PDEs system:

\begin{matrix}∂_{t}ρ + ∂_{x}J = 0 , \\ ε^{2} ∂_{t}J + ∂_{x} (ρ ψ(\frac{εJ}{ρ} )) = −J , \end{matrix}  (1.1)

where the unknown quantities are the density ρ and the current J, while ε is a positive scaling parameter, and ψ is a smooth function whose definition will be given later on. Our motivation for studying (1.1) comes from kinetic theory. The system (1.1) is indeed derived as an intermediate model between a fully microscopic description and its diffusion asymptotics ε → 0. Let f(t, x, v) ≥ 0 stand for a (microscopic) density of particles having at time t a position x ∈ ℝ and a velocity v ∈ (−1, +1). In all the paper, dμ stands for the normalized Lebesgue measure on (−1, +1), that is dμ = dv/2. Let us assume that the evolution of f is governed by the linear kinetic equation:

ε ∂_{t}f + v ∂_{x}f = \frac{1}{ε} Q(f) , Q(f)(v) = ∫_{−1}^{+1} f(v′) dμ(v′) − f(v) .  (1.2)

The parameter ε is related to the mean free path: 0 < ε ≪ 1 means that particles suffer more and more collision events per time unit. As a consequence of the penalization of the collision term in (1.2), we guess that for small values of ε, the microscopic density f tends to an element of Ker(Q), which is reduced to constants with respect to the velocity variable:

f(t, x, v) ≃ ϱ(t, x) .  (1.3)

Let us introduce the moments of the density f:

(\begin{matrix}ρ \\ J \\ p \end{matrix}) (t, x) := ∫_{−1}^{+1} (\begin{matrix}1 \\ v/ε \\ v^{2} \end{matrix}) f(t, x, v) dμ(v) .

Integration of (1.2) yields the following system:

\begin{matrix}∂_{t}ρ + ∂_{x}J = 0 , \\ ε^{2}∂_{t}J + ∂_{x}p = −J . \end{matrix}  (1.4)

The formal ansatz (1.3) leads to p(t, x) ≃ ϱ(t, x) ∫ v^{2} dμ(v) = ϱ(t, x)/3, so that formally, as ε tends to 0, (1.4) becomes:

∂_{t}ϱ + ∂_{x}J = 0 , \frac{1}{3} ∂_{x}ϱ = −J .

Hence, the limit density ϱ should satisfy the heat equation:

∂_{t}ϱ − \frac{1}{3} ∂_{xx}^{2}ϱ = 0 .  (1.5)

The convergence of the solutions of (1.2) to the solution ϱ of the heat equation (1.5) can be rigorously established, for this simple model as well as for much more complicated (linear and nonlinear) kinetic models, see e.g. [2].

An important issue, for instance in nuclear engineering or in radiative transfer where this question naturally arises, consists in describing, with enough accuracy, intermediate regimes for small but nonzero ε. Clearly the solution of the heat equation does not provide such a nice approximation, first of all because it gets rid of the velocity variable. Note also that (1.5) propagates information with infinite speed while in (1.2), the speed of propagation does not exceed 1/ε (that is large but finite!). To examine intermediate regimes, a possible strategy consists in closing the moment system (1.4), by expressing in a suitable way the second moment p as a function of the zero and first moments ρ, and J. If the resulting system is found to be hyperbolic, we shall keep a finite speed of propagation. There exists a huge variety of such closure methods. The system (1.1) is derived by an entropy minimization principle, as described by Levermore in [14]. This is the reason why we shall always refer to (1.1) as Levermore’s system. Let us set:

H(f) = ∫_{−1}^{+1} f(v) ln f(v) dμ(v) ,

and observe that H is dissipated by (1.2). Then, for (ρ, J) given in ℝ^{+} × ℝ, let us define the function v ↦ F (v) such that:

H(F) = min \{H(g) , g : (−1, +1) → ℝ^{+} , ∫_{−1}^{+1} (\begin{matrix}1 \\ v/ε \end{matrix}) g(v) dμ(v) = (\begin{matrix}ρ \\ J \end{matrix})\} .

Defining:

∀ β ∈ ℝ , 𝔽(β) := \frac{sinh(β)}{β} ,  𝔾(β) := coth(β) − \frac{1}{β} = \frac{𝔽′(β)}{𝔽(β)} ,  (1.6)

the above minimizer F is well-defined provided that ε|J| < ρ, and is given by the formula:

F (v) = \frac{ρ}{𝔽 ∘ 𝔾^{−1}(εJ/ρ)} exp (v 𝔾^{−1}(\frac{εJ}{ρ} )) .

Observe that 𝔾 is an increasing C^{∞} diffeomorphism from ℝ onto (−1, 1), so the use of the inverse 𝔾^{−1} is legitimate. We also define a function ψ by the formula:

\begin{aligned}ψ : & (−1, +1) & ⟶ & ]0, +∞[ \\ & u & ⟼ & u^{2} + 𝔾′(𝔾^{−1}(u)) = \frac{𝔽′′}{𝔽} (𝔾^{−1}(u)) . \end{aligned}  (1.7)

Eventually, we close (1.4) by requiring p to be the second moment of the minimizer F:

p = ∫_{−1}^{+1} v^{2} F (v) dμ(v) = ρ ψ(\frac{εJ}{ρ} ) ,

and we thus get the system (1.1), with ψ defined by (1.7). Of course, the closure of the moments system (1.4) highly depends on the set in which the velocity variable lies. For instance, if one replaces the set (−1, +1), equipped with the normalized Lebesgue measure, by the whole line ℝ, equipped with the Gaussian measure exp(−v^{2}/2) dv/\sqrt{2π}, then the entropy minimization procedure leads to the isothermal Euler system with relaxation:

\begin{matrix}∂_{t}ρ + ∂_{x}J = 0 , \\ ε^{2}(∂_{t}J + ∂_{x} \frac{J^{2}}{ρ} ) + ∂_{x}ρ = −J . \end{matrix}

In [7], we have shown that (1.1) admits global smooth solutions that are bounded away from vacuum, and that the density converges towards the solution to the heat equation (1.5) as ε tends to zero. The intermediate system (1.1) is thus consistent with the diffusion limit (1.5). However, this first result is not fully satisfactory because both the kinetic equation (1.2), and the heat equation (1.5) are perfectly defined up to ρ = 0, and it is therefore natural to try to justify a similar asymptotic result for solutions to (1.1) that may contain vacuum regions. More precisely, we would like to show the existence of global weak solutions of (1.1) for bounded initial data that may contain vacuum, and derive uniform L^{∞} bounds with respect to ε. As we shall see later on, vacuum corresponds to the boundary of the domain of hyperbolicity of (1.1). In view of earlier results, see e.g. [9, 16, 15, 17, 4, 5, 13], it seems natural to solve our problem by showing that the vanishing viscosity method for (1.1) is convergent. A key ingredient in proving this convergence is to obtain a priori L^{∞} bounds that are uniform with respect to the vanishing viscosity coefficient. Such bounds can be derived by showing the positive invariance of some domains, using the result of [6]. In this paper, we compute such positive invariant domains, and we postpone the complete justification of the vanishing viscosity method to a future work. There is a similar question for the homogeneous Riemann problem, that is for the system (1.1) without relaxation term. We shall show that this homogeneous problem can be solved without smallness assumption, and we shall compute invariant regions. The main result of this paper is the following:

Theorem 1.1. Let ρ_{l}, ρ_{r} > 0, and J_{l}, J_{r} ∈ ℝ be such that ε|J_{l}| < ρ_{l} and ε|J_{r}| < ρ_{r}. Then the homogeneous Riemann problem:

\begin{matrix}∂_{t}ρ + ∂_{x}J = 0 , \\ ε^{2} ∂_{t}J + ∂_{x} (ρ ψ(\frac{εJ}{ρ} )) = 0 , \end{matrix}

with initial data:

(ρ, J)_{|t=0} = \begin{cases}(ρ_{l}, J_{l}) & if x < 0, \\ (ρ_{r}, J_{r}) & if x > 0, \end{cases}

is uniquely solvable, and its solution does not contain vacuum. Moreover, for all ρ_{0} > 0, and for all J_{0} that satisfies ε|J_{0}| < ρ_{0}, there exists a region I_{ρ_{0},J_{0}} ⊂ ℝ^{2} that contains (ρ_{0}, J_{0}) and that is invariant for the Riemann problem.

For all ρ_{0} > 0, and all ν > 0, there exists a region that contains (ρ_{0}, 0) and that is positively invariant for the viscous approximation:

\begin{matrix}∂_{t}ρ + ∂_{x}J = ν ∂_{xx}ρ , \\ ε^{2} ∂_{t}J + ∂_{x} (ρ ψ(\frac{εJ}{ρ})) = −J + ν ∂_{xx}J . \end{matrix}  (1.8)

The precise definition of the invariant regions requires some notations, so we postpone it to the next sections. We can also deal with the case where the Riemann data for (1.1) contains vacuum:

Proposition 1.1. Let ρ_{r} > 0, and J_{r} ∈ ℝ be such that ε|J_{r}| < ρ_{r}. Then the homogeneous Riemann problem:

\begin{matrix}∂_{t}ρ + ∂_{x}J = 0 , \\ ε^{2} ∂_{t}J + ∂_{x} (ρ ψ(\frac{εJ}{ρ} )) = 0 , \end{matrix}

with initial data:

(ρ, J)_{|t=0} = \begin{cases}(0, 0) & if x < 0, \\ (ρ_{r}, J_{r}) & if x > 0, \end{cases}

is uniquely solvable, and its solution is a rarefaction wave that connects the states (0, 0) and (ρ_{r}, J_{r}). An analogous result holds when the right state is the vacuum, and the left state has a positive density.

To contruct global weak solutions to the system (1.1) by the vanishing viscosity method, one introduces its regularized version (1.8), and one tries to pass to the limit ν → 0^{+}. For the isentropic gas dynamics system (with no relaxation term), this procedure was first justified by DiPerna [9]. DiPerna’s work was later extended in [16, 15, 4, 5] to more general pressure laws. We also refer to [8] for the case with a source term. In the isothermal case (where the pressure is a linear function of the density), the analysis of [16, 15, 4, 5] does not apply, due to the singular nature of the equation for the entropies (see Section 4). Nevertheless, the existence of global weak solutions with vacuum for the isothermal Euler equations was proved recently in [13] (see also [11]). Since the system (1.1) is ”close” to the isothermal Euler system with a relaxation, it is really tempting to believe that the vanishing viscosity method (1.8) can also be justified for Levermore’s system (1.1).

The rest of this paper is organized as follows. In Section 2, we prove some general properties of the relaxation free system corresponding to (1.1). We show that this system is strictly hyperbolic, and that its characteristic fields are genuinely nonlinear. We also exhibit some Riemann invariants. In Section 3, we investigate the Riemann problem, and show that it is uniquely solvable with invariant regions. In Section 4, we prove the positive invariance of some regions for the parabolic regularization (1.8).

2 Some general facts on Levermore’s system

In this section, we focus on the system (1.1) with no relaxation term:

\begin{matrix}∂_{t}ρ + ∂_{x}J = 0 , \\ ε^{2} ∂_{t}J + ∂_{x} (ρ ψ(\frac{εJ}{ρ})) = 0 . \end{matrix}  (2.1)

The functions ψ, 𝔽, and 𝔾 are defined by (1.7), (1.6). In particular, 𝔽 and ψ are even functions, while 𝔾 is an odd function. The following relations will be often used throughout the paper:

𝔽(0) = 1 ,  𝔾(0) = 0 ,  ψ(0) = 𝔾′(0) = \frac{1}{3} ,  ψ′(0) = 0 .

In all what follows, we use the notation u := εJ/ρ ∈ (−1, 1) to denote the (rescaled) velocity. System (2.1) meets the classical hyperbolicity properties, as shown in the following:

Proposition 2.1. The system (2.1) is stricly hyperbolic in the convex open set {(ρ, J)/ρ > 0 , ε |J| < ρ}. Its characteristic speeds λ_{1}^{ε}, λ_{2}^{ε} are given by:

λ_{i}^{ε}(ρ, J) := \frac{1}{ε} λ_{i}(u) , i = 1, 2 ,

λ_{1}(u) := \frac{ψ′(u) − \sqrt{ψ′(u)^{2} − 4uψ′(u) + 4ψ(u)}}{2} , λ_{2}(u) := \frac{ψ′(u) + \sqrt{ψ′(u)^{2} − 4uψ′(u) + 4ψ(u)}}{2} , (2.2)

and the corresponding eigenvectors can be defined as follows:

r_{i}^{ε}(ρ, J) = (\begin{matrix}1 \\ λ_{i}^{ε}(ρ, J) \end{matrix}) , i = 1, 2 .  (2.3)

Moreover, the function:

ℍ(ρ, J) := ρ ln ρ − ρ ln [𝔽 ∘ 𝔾^{−1}(u)] + εJ 𝔾^{−1}(u) ,  (2.4)

is a strictly convex entropy for (2.1).

Proof. If we write system (2.1) under the compact form:

∂_{t} (\begin{matrix}ρ \\ J \end{matrix}) + ∂_{x}F_{ε}(ρ, J) = 0 ,

we can compute the Jacobian matrix:

DF_{ε}(ρ, J) = (\begin{matrix}0 & 1 \\ \frac{1}{ε^{2}} (ψ(u) − uψ′(u)) & \frac{1}{ε}ψ′(u) \end{matrix}) .

It is now easy to check that the two eigenvalues of DF_{ε}(ρ, J) are real, distinct, and given by (2.2). The discriminant of the characteristic polynomial is positive since

ψ′(u)^{2} + 4(ψ(u) − uψ′(u)) = (ψ′(u) − 2u)^{2} + 4(ψ(u) − u^{2}) ≥ 4(ψ(u) − u^{2}) > 0 .

Thus, the system (2.1) is strictly hyperbolic. The reader will easily check that, up to a multiplicative constant, the eigenvectors of DF_{ε}(ρ, J) are given by (2.3).

That ℍ is a strictly convex entropy for (2.1) was already checked in [7], where the corresponding flux is also computed.

Observe that the characteristic speeds of (2.1) only depend on the velocity u. This property is shared by the isothermal Euler system. Another analogy with the isothermal Euler system is the following:

Lemma 2.1. Let λ_{1,2} be defined by (2.2). Then we have:

∀ u ∈ (−1, 1) , λ_{1}(u) < u < λ_{2}(u) .

The proof is an easy consequence of the definition (2.2). Thanks to Lemma 2.1, we are able to define two functions

Λ_{i}(u) := ∫_{0}^{u} \frac{dv}{λ_{i}(v) − v} , i = 1, 2 , u ∈ (−1, 1) .  (2.5)

It is clear that Λ_{1} is decreasing, while Λ_{2} is increasing. These functions are useful to construct Riemann invariants of (2.1):

Lemma 2.2. Let Λ_{1}, Λ_{2} be defined by (2.5), and let us define:

Z_{i}(ρ, J) := − ln ρ + Λ_{i}(\frac{εJ}{ρ} ) , i = 1, 2 .

Then Z_{1}, Z_{2} are Riemann invariants for (2.1):

∇Z_{1} · r_{1}^{ε} = 0 , ∇Z_{2} · r_{2}^{ε} = 0 ,

where the eigenvectors r_{1}^{ε}, r_{2}^{ε} are defined by (2.3).

The proof is again a basic application of the chain rule, and we omit it. We end this section with a result on the genuine nonlinearity of the characteristic fields:

Proposition 2.2. The function ψ is strictly convex and satisfies

lim_{u→±1} ψ(u) = 1 , lim_{u→±1} ψ′(u) = ±2 .

Furthermore λ_{1} and λ_{2} satisfy

∀ u ∈ (−1, 1) , λ′_{i}(u) > 0 , i = 1, 2 ,

and both characteristic fields of (2.1) are genuinely nonlinear.

Proof. We first show the strict convexity of ψ. Using the definition (1.7), we compute:

ψ′′(u) = 2 + \frac{𝔾′′′}{(𝔾′)^{2}} (𝔾^{−1}(u)) − \frac{(𝔾′′)^{2}}{(𝔾′)^{3}} (𝔾^{−1}(u)) .

Therefore, ψ is strictly convex if, and only if the following condition holds:

∀ β ∈ ℝ , 2 𝔾′(β)^{3} + 𝔾′(β) 𝔾′′′(β) − 𝔾′′(β)^{2} > 0 .

Using the relation 𝔾(β) = 𝔽′(β)/𝔽(β), and recalling that 𝔽(β) is positive, see (1.6), we obtain that ψ is strictly convex if, and only if, we have

∀ β ∈ ℝ , −𝔽′′(β)^{3}+𝔽(β) 𝔽′′(β) 𝔽^{(4)}(β)+2𝔽′(β) 𝔽′′(β) 𝔽′′′(β)−𝔽′(β)^{2} 𝔽^{(4)}(β)−𝔽(β) 𝔽′′′(β)^{2} > 0 .

We now use the expression (1.6) of 𝔽, and are thus reduced to proving the following inequality:

∀ β ∈ ℝ , Υ(β) := \frac{1}{β^{9}} \{(sinh β)^{3} + 2 β^{3} cosh β − 3 β^{2} sinh β − β^{4} sinh β\} > 0 .  (2.6)

A numerical evidence of this inequality is shown in Figure 1, where we have depicted the function Υ on the interval [0, 100]. For the sake of completeness, we are going to prove that Υ(β) is positive for all β, by expanding Υ in power series. (Observe that Υ can indeed be expanded in power series, and the radius of convergence is +∞). After some simplifications, we are led to

Υ(β) = ∑_{k=0}^{+∞} υ_{k} β^{2k} , υ_{k} := \frac{19683 9^{k} − 8931 − 10552 k − 4640 k^{2} − 896 k^{3} − 64 k^{4}}{4 (2k + 9)!} .

Figure 1: The function Υ on [0, 100] (with a log-log scale).

The reader will easily check that υ_{0}, and υ_{1} are positive. When k ≥ 2, one has

\frac{19683 9^{k}}{k^{4}} ≥ \frac{19683 9^{2}}{2^{4}} ≥ 64 + \frac{896}{2} + \frac{4640}{4} + \frac{10552}{8} + \frac{8931}{16}

≥ 64 + \frac{896}{k} + \frac{4640}{k^{2}} + \frac{10552}{k^{3}} + \frac{8931}{k^{4}} .

Consequently, the υ_{k}’s are nonnegative for all k ≥ 1, and Υ(β) is bounded from below by υ_{0} > 0. This proves that ψ is strictly convex.

From the definition (1.6) of the function 𝔾, we have

lim_{|β|→+∞} 𝔾′(β) = 0 , lim_{u→±1} 𝔾^{−1}(u) = ±∞ .

Since we have ψ(u) = u^{2} + 𝔾′ ∘ 𝔾^{−1}(u), we get

lim_{|u|→1} ψ(u) = 1 .

To compute the limit of ψ′(u) at ±1, we use the expression

ψ′(u) = 2u + \frac{𝔾′′}{𝔾′} [𝔾^{−1}(u)] .  (2.7)

The function 𝔾 satisfies

𝔾′(β) = \frac{1}{β^{2}} − \frac{1}{(sinh β)^{2}} , 𝔾′′(β) = \frac{−2}{β^{3}} + \frac{2 cosh β}{(sinh β)^{3}} ,

which implies

lim_{|β|→+∞} \frac{𝔾′′}{𝔾′} (β) = 0 .

This yields the desired limit for ψ′(u) at ±1.

We are now going to show that the characteristic fields are genuinely nonlinear. Using the definition (2.2), the λ_{i}’s satisfy the second order polynomial equation

∀ u ∈ (−1, 1) , λ_{i}(u)(λ_{i}(u) − ψ′(u)) − (ψ(u) − uψ′(u)) = 0 .

Differentiating this equation with respect to u, and using the relations

−[2λ_{1}(u) − ψ′(u)] = 2λ_{2}(u) − ψ′(u) = \sqrt{Δ(u)} ,

where

∀ u ∈ (−1, 1) , Δ(u) := ψ′(u)^{2} + 4(ψ(u) − uψ′(u)) > 0 ,  (2.8)

we end up with

λ′_{1}(u) \sqrt{Δ(u)} = ψ′′(u) (u − λ_{1}(u)) , λ′_{2}(u) \sqrt{Δ(u)} = ψ′′(u) (λ_{2}(u) − u) .

Using Lemma 2.1 and the strict convexity of ψ, we get λ′_{1}(u) > 0, and λ′_{2}(u) > 0. Using (2.2) and (2.3), we compute

∇λ_{i}^{ε} · r_{i}^{ε} = \frac{1}{ερ} λ′_{i}(u) (λ_{i}(u) − u) ≠ 0 ,

so both characteristic fields are genuinely nonlinear.

Because the characteristic fields are genuinely nonlinear, we can expect to solve the Riemann problem for (2.1) with a juxtaposition of shock waves and rarefaction waves (see e.g. [3, 12, 18] for the general theory). In the next section, we are going to provide some details on the Rankine-Hugoniot relations, and exhibit the 1-Lax shocks and the 2-Lax shocks (meaning the discontinuities satisfying Lax’ shock inequalities). The rarefaction curves will be computed by using the Riemann invariants of Lemma 2.2. Then we shall show that the Riemann problem for (2.1) is uniquely solvable, without any restriction on the data.

3 Global solvability of the Riemann problem

3.1 The shock curves

In this paragraph, we are going to compute the shock curves for the hyperbolic system (2.1). We fix, once and for all, a state (ρ_{−}, J_{−}) in the domain of hyperbolicity. We are going to compute the states (ρ_{+}, J_{+}) such that there exists a speed s ∈ ℝ verifying the Rankine-Hugoniot relations for the system (2.1). The Rankine-Hugoniot relations read

\begin{matrix}J_{+} − J_{−} = s (ρ_{+} − ρ_{−}) , \\ \frac{1}{ε^{2}} [ρ_{+} ψ(u_{+}) − ρ_{−} ψ(u_{−})] = s (J_{+} − J_{−}) . \end{matrix}  (3.1)

In (3.1), we have used the notations u_{±} = εJ_{±}/ρ_{±}.

If ρ_{+} = ρ_{−}, we get the obvious solution J_{+} = J_{−}. We thus assume from now on ρ_{+} ≠ ρ_{−}, and we rewrite the Rankine-Hugoniot relations (3.1) in the following way:

s = \frac{1}{ε} \frac{ρ_{+}u_{+} − ρ_{−}u_{−}}{ρ_{+} − ρ_{−}} , (ρ_{+} − ρ_{−})(ρ_{+}ψ(u_{+}) − ρ_{−}ψ(u_{−})) = (ρ_{+}u_{+} − ρ_{−}u_{−})^{2} .  (3.2)

We have the following result:

Proposition 3.1. When u_{+} ≠ u_{−}, the solutions of (3.2) consist in the two following curves, that are parametrized by u_{+}:

(ρ_{+}, J_{+}, s) = (ρ_{−} X_{1}(u_{−}, u_{+}), \frac{ρ_{−} X_{1}(u_{−}, u_{+}) u_{+}}{ε} , \frac{1}{ε} \frac{X_{1}(u_{−}, u_{+}) u_{+} − u_{−}}{X_{1}(u_{−}, u_{+}) − 1} ) ,

(ρ_{+}, J_{+}, s) = (ρ_{−} X_{2}(u_{−}, u_{+}), \frac{ρ_{−} X_{2}(u_{−}, u_{+}) u_{+}}{ε} , \frac{1}{ε} \frac{X_{2}(u_{−}, u_{+}) u_{+} − u_{−}}{X_{2}(u_{−}, u_{+}) − 1} ) ,

where u_{+} ∈ (−1, u_{−}) ∪ (u_{−}, 1), and where X_{1}(u_{−}, u_{+}), X_{2}(u_{−}, u_{+}) are the roots of the second order polynomial equation

P(X) := X^{2} (ψ(u_{+}) − u_{+}^{2}) − X (ψ(u_{+}) + ψ(u_{−}) − 2u_{−}u_{+}) + (ψ(u_{−}) − u_{−}^{2}) = 0 ,  (3.3)

and satisfy

0 < X_{2}(u_{−}, u_{+}) < 1 < X_{1}(u_{−}, u_{+}) .

Proof. Defining X := ρ_{+}/ρ_{−}, the second equation of (3.2) is found to be equivalent to (3.3). Recall that ψ(u) − u^{2} > 0, see (1.7). We thus only need to show that P has two positive real roots when u_{+} ≠ u_{−}. The discriminant δ̃(u_{−}, u_{+}) of P is given by the relation:

δ̃(u_{−}, u_{+}) = (u_{+} − u_{−})^{2} δ(u_{−}, u_{+}) ,

δ(u_{−}, u_{+}) := z^{2} − 4u_{+}z + 4ψ(u_{+}) = z^{2} − 4u_{−}z + 4ψ(u_{−}) , with z := \frac{ψ(u_{+}) − ψ(u_{−})}{u_{+} − u_{−}} . (3.4)

Since we have

δ̃(u_{−}, u_{+}) = (u_{+} − u_{−})^{2} [(z − 2u_{+})^{2} + 4(ψ(u_{+}) − u_{+}^{2})] > 0 ,

P has two distinct real roots X_{2}(u_{−}, u_{+}) < X_{1}(u_{−}, u_{+}). (The reader will understand the choice of our notations when we shall prove Lax’ shock inequalities). Note that these roots are simple when u_{+} ≠ u_{−}, therefore X_{1}, X_{2} are C^{∞} functions of (u_{−}, u_{+}) by the implicit functions Theorem. Moreover, we have

X_{1}X_{2} > 0 , (ψ(u_{+}) − u_{+}^{2})(X_{1} + X_{2}) = ψ(u_{+}) + ψ(u_{−}) − 2u_{−}u_{+} > u_{+}^{2} + u_{−}^{2} − 2u_{−}u_{+} > 0 ,

P(1) = −(u_{+} − u_{−})^{2} < 0 ,

and we thus have the inequalities

0 < X_{2} < 1 < X_{1} ,

as long as u_{+} ≠ u_{−}. The proof is complete.

We have thus constructed all the nontrivial solutions to (2.1) of the form

(ρ, J) = \begin{cases}(ρ_{−}, J_{−}) & if x < s t , \\ (ρ_{+}, J_{+}) & if x > s t . \end{cases}

We now need to characterize those discontinuous solutions that are admissible, in the sense that they satisfy Lax’ shock inequalities, see [12]. As expected from the general theory, only half of each curve defined above is admissible:

Proposition 3.2. When u_{+} < u_{−}, Lax’ shock inequalities hold, that is, we have:

\frac{1}{ε} λ_{1}(u_{+}) < \frac{1}{ε} \frac{X_{1}(u_{−}, u_{+}) u_{+} − u_{−}}{X_{1}(u_{−}, u_{+}) − 1} < \frac{1}{ε} λ_{2}(u_{+}) , \frac{1}{ε} \frac{X_{1}(u_{−}, u_{+}) u_{+} − u_{−}}{X_{1}(u_{−}, u_{+}) − 1} < \frac{1}{ε} λ_{1}(u_{−}) ,

\frac{1}{ε} λ_{1}(u_{−}) < \frac{1}{ε} \frac{X_{2}(u_{−}, u_{+}) u_{+} − u_{−}}{X_{2}(u_{−}, u_{+}) − 1} < \frac{1}{ε} λ_{2}(u_{−}) , \frac{1}{ε} λ_{2}(u_{+}) < \frac{1}{ε} \frac{X_{2}(u_{−}, u_{+}) u_{+} − u_{−}}{X_{2}(u_{−}, u_{+}) − 1} .

When u_{+} > u_{−}, Lax’ shock inequalities do not hold.

Proof. We first assume u_{+} > u_{−}. In this case, Proposition 2.2 yields λ_{i}(u_{+}) > λ_{i}(u_{−}), for i = 1, 2. Therefore, Lax’ shock inequalities cannot hold when u_{+} > u_{−}, and the discontinuity is not admissible.

From now on, we assume u_{+} ∈ (−1, u_{−}), and we want to prove Lax’ shock inequalities, that is:

λ_{1}(u_{+}) < \frac{X_{1}(u_{−}, u_{+}) u_{+} − u_{−}}{X_{1}(u_{−}, u_{+}) − 1} < λ_{2}(u_{+}) , \frac{X_{1}(u_{−}, u_{+}) u_{+} − u_{−}}{X_{1}(u_{−}, u_{+}) − 1} < λ_{1}(u_{−}) .

Recall that the velocities λ_{1}, λ_{2} are defined in Proposition 2.1, and that X_{1}(u_{−}, u_{+}) > 1. a) We first show the inequality

[X_{1}(u_{−}, u_{+}) − 1] λ_{1}(u_{+}) < X_{1}(u_{−}, u_{+})u_{+} − u_{−} .  (3.5)

From the definition of X_{1}(u_{−}, u_{+}), see (3.3) and (3.4) for the definition of δ and z, we obtain:

\frac{X_{1}(u_{−}, u_{+}) − 1}{u_{−} − u_{+}} = \frac{z − 2u_{+} + \sqrt{δ(u_{−}, u_{+})}}{2(ψ(u_{+}) − u_{+}^{2})} ,

\frac{X_{1}(u_{−}, u_{+})u_{+} − u_{−}}{u_{−} − u_{+}} = \frac{u_{+}z − 2ψ(u_{+}) + u_{+}\sqrt{δ(u_{−}, u_{+})}}{2(ψ(u_{+}) − u_{+}^{2})} .

Simplifying (3.5) with the help of these two relations, it turns out that (3.5) is equivalent to

(λ_{1}(u_{+}) − u_{+})z + 2(ψ(u_{+}) − u_{+}λ_{1}(u_{+})) < (u_{+} − λ_{1}(u_{+}))\sqrt{δ(u_{−}, u_{+})} .  (3.6)

Because u_{+} > λ_{1}(u_{+}) (see Lemma 2.1), (3.6) will be satisfied provided that we have the weaker inequality

[(λ_{1}(u_{+}) − u_{+})z + 2(ψ(u_{+}) − u_{+}λ_{1}(u_{+}))]^{2} < [u_{+} − λ_{1}(u_{+})]^{2} δ(u_{−}, u_{+}) .  (3.7)

Expanding (3.7), and using (3.4), (3.7) can be rewritten as

z (u_{+} − λ_{1}(u_{+})) − ψ(u_{+}) + λ_{1}(u_{+})^{2} > 0 .  (3.8)

Because ψ is strictly convex, and u_{+} < u_{−}, we obtain (we use (2.2) for the last equality):

z (u_{+} − λ_{1}(u_{+})) − ψ(u_{+}) + λ_{1}(u_{+})^{2} > ψ′(u_{+}) (u_{+} − λ_{1}(u_{+})) − ψ(u_{+}) + λ_{1}(u_{+})^{2} = 0 .

As a consequence, (3.8) holds, and therefore (3.5) is proved.

b) We now prove the inequality

X_{1}(u_{−}, u_{+}) u_{+} − u_{−} < [X_{1}(u_{−}, u_{+}) − 1] λ_{2}(u_{+}) .

This inequality is seen to be equivalent to

(u_{+} − λ_{2}(u_{+}))z − 2(ψ(u_{+}) − u_{+}λ_{2}(u_{+})) < (λ_{2}(u_{+}) − u_{+})\sqrt{δ(u_{−}, u_{+})} .  (3.9)

Because λ_{2}(u_{+}) > u_{+}, the right-hand side of (3.9) is positive, and (3.9) will hold provided that the left-hand side is nonpositive, that is,

(λ_{2}(u_{+}) − u_{+})z + 2(ψ(u_{+}) − u_{+}λ_{2}(u_{+})) ≥ 0 .  (3.10)

Observe that we have

(λ_{2}(u_{+}) − u_{+}) z > (λ_{2}(u_{+}) − u_{+}) ψ′(u_{+}) .

Using the definition (2.2) of λ_{2}(u_{+}), and the definition (2.8) of Δ, we can write

λ_{2}(u_{+}) = \frac{ψ′(u_{+}) + \sqrt{Δ(u_{+})}}{2} .

Then, one checks that

(λ_{2}(u_{+}) − u_{+}) ψ′(u_{+}) + 2(ψ(u_{+}) − u_{+}λ_{2}(u_{+})) = Δ(u_{+}) + (ψ′(u_{+}) − 2u_{+})\sqrt{Δ(u_{+})} ≥ 0 .

Hence (3.10) holds, and the second Lax inequality is proved.

c) Eventually, we prove the inequality

\frac{X_{1}(u_{−}, u_{+})u_{+} − u_{−}}{X_{1}(u_{−}, u_{+}) − 1} < λ_{1}(u_{−}) .  (3.11)

Differentiating with respect to u ∈ (−1, u_{−}), one checks that the function

u ⟼ \frac{X_{1}(u_{−}, u)u − u_{−}}{X_{1}(u_{−}, u) − 1}

is strictly increasing. To prove (3.11), it is therefore sufficient to show the inequality

lim_{u_{+}→u_{−}} \frac{X_{1}(u_{−}, u_{+})u_{+} − u_{−}}{X_{1}(u_{−}, u_{+}) − 1} = λ_{1}(u_{−}) .  (3.12)

From the equality

\frac{X_{1}(u_{−}, u_{+})u_{+} − u_{−}}{X_{1}(u_{−}, u_{+}) − 1} = \frac{u_{+}z − 2ψ(u_{+}) + u_{+}\sqrt{δ(u_{−}, u_{+})}}{z − 2u_{+} + \sqrt{δ(u_{−}, u_{+})}} , z = \frac{ψ(u_{+}) − ψ(u_{−})}{u_{+} − u_{−}} ,

we compute

lim_{u_{+}→u_{−}} \frac{X_{1}(u_{−}, u_{+})u_{+} − u_{−}}{X_{1}(u_{−}, u_{+}) − 1} = \frac{u_{−}ψ′(u_{−}) − 2ψ(u_{−}) + u_{−}\sqrt{Δ(u_{−})}}{ψ′(u_{−}) − 2u_{−} + \sqrt{Δ(u_{−})}}

= \frac{ψ′(u_{−}) − \sqrt{Δ(u_{−})}}{2} = λ_{1}(u_{−}) ,

and (3.12) holds. We have thus proved (3.11), and this concludes the proof of the 1-Lax shock inequalities. The 2-Lax shock inequalities are proved in a completely similar way, and we shall not reproduce the calculations.

Remark 3.1. It would be interesting to know whether the admissible shocks also satisfy the entropy inequality:

s(ℍ(ρ_{+}, J_{+}) − ℍ(ρ_{−}, J_{−})) − (q(ρ_{+}, J_{+}) − q(ρ_{−}, J_{−})) ≥ 0 ,

where ℍ is the strictly convex entropy defined by (2.4), and q is the corresponding flux. From the general theory, we know that this inequality holds for small amplitude shocks, that is, when u_{+} < u_{−}, and u_{+} close to u_{−}, see e.g. [3]. Proving the entropy inequality for shocks of arbitrarily large amplitude seems to be particularly challenging.

The following Lemma will be useful when solving the Riemann problem:

Lemma 3.1. The functions X_{1} and X_{2} satisfy

∀ u_{−} ∈ (−1, 1) , lim_{u_{+}→−1} X_{1}(u_{−}, u_{+}) = +∞ ,

∀ u_{+} ∈ (−1, 1) , lim_{u_{−}→1} X_{2}(u_{−}, u_{+}) = 0 .

Proof. When u_{+} < u_{−}, the definition of X_{1} is

X_{1}(u_{−}, u_{+}) = \frac{ψ(u_{+}) + ψ(u_{−}) − 2u_{−}u_{+} + (u_{−} − u_{+})\sqrt{δ(u_{−}, u_{+})}}{2(ψ(u_{+}) − u_{+}^{2})} ,

see Proposition 3.1, and (3.4). We have

ψ(u_{+}) + ψ(u_{−}) − 2u_{−}u_{+} + (u_{−} − u_{+})\sqrt{δ(u_{−}, u_{+})} ≥ ψ(u_{+}) + ψ(u_{−}) − 2u_{−}u_{+}

≥ ψ(u_{+}) − u_{+}^{2} + ψ(u_{−}) − u_{−}^{2} ≥ ψ(u_{−}) − u_{−}^{2} > 0 ,

and we also have

lim_{u_{+}→−1}[ψ(u_{+}) − u_{+}^{2}] = 0^{+} .

The first part of the Lemma follows.

For u_{+} < u_{−}, we also have

X_{2}(u_{−}, u_{+}) = \frac{ψ(u_{+}) + ψ(u_{−}) − 2u_{−}u_{+} − \sqrt{δ̃(u_{−}, u_{+})}}{2(ψ(u_{+}) − u_{+}^{2})} ,

with

δ̃(u_{−}, u_{+}) = [ψ(u_{+}) + ψ(u_{−}) − 2u_{−}u_{+}]^{2} − 4(ψ(u_{+}) − u_{+}^{2})(ψ(u_{−}) − u_{−}^{2}) .

We thus compute

ψ(u_{+}) + ψ(u_{−}) − 2u_{−}u_{+} − \sqrt{δ̃(u_{−}, u_{+})} → 0 , when u_{−} → 1 ,

and the proof is complete.

3.2 The rarefaction curves

In this paragraph, we construct the rarefaction curves. Recall that a i-rarefaction wave (i = 1, 2), connecting two states (ρ_{−}, J_{−}) and (ρ_{+}, J_{+}), is a solution to (2.1) of the form:

(ρ, J)(t, x) = \begin{cases}(ρ_{−}, J_{−}) & if x ≤ λ_{i}^{ε}(ρ_{−}, J_{−})t , \\ W (x/t) & if λ_{i}^{ε}(ρ_{−}, J_{−})t ≤ x ≤ λ_{i}^{ε}(ρ_{+}, J_{+})t , \\ (ρ_{+}, J_{+}) & if λ_{i}^{ε}(ρ_{+}, J_{+})t ≤ x , \end{cases}

with λ_{i}^{ε}(ρ_{−}, J_{−}) ≤ λ_{i}^{ε}(ρ_{+}, J_{+}), and where W (ξ) is a solution to the ODE:

W′(ξ) = \frac{r_{i}^{ε}(W (ξ))}{∇λ_{i}^{ε}(W (ξ)) · r_{i}^{ε}(W (ξ))} , ξ ∈ [λ_{i}^{ε}(ρ_{−}, J_{−}), λ_{i}^{ε}(ρ_{+}, J_{+})] .

Thanks to Lemma 2.2, we know that the integral curves of the vector field r_{1}^{ε} (resp. r_{2}^{ε}) are nothing but the level sets of the Riemann invariant Z_{1} (resp. Z_{2}). Moreover, the constraint λ_{i}^{ε}(ρ_{−}, J_{−}) ≤ λ_{i}^{ε}(ρ_{+}, J_{+}) is equivalent to u_{−} ≤ u_{+}, because the functions λ_{i} are increasing (see Lemma 2.2). We thus have:

Proposition 3.3. Two states (ρ_{−}, J_{−}) and (ρ_{+}, J_{+}) are connected by a 1-rarefaction wave if and only if

u_{−} ≤ u_{+} , and ρ_{+} = ρ_{−} exp ∫_{u_{−}}^{u_{+}} \frac{dw}{λ_{1}(w) − w} .

Two states (ρ_{−}, J_{−}) and (ρ_{+}, J_{+}) are connected by a 2-rarefaction wave if and only if

u_{−} ≤ u_{+} , and ρ_{+} = ρ_{−} exp ∫_{u_{−}}^{u_{+}} \frac{dw}{λ_{2}(w) − w} .

Moreover, we have the following asymptotic behavior for the rarefaction curves:

Lemma 3.2. We have

∫_{−1}^{0} \frac{dw}{w − λ_{1}(w)} = ∫_{0}^{1} \frac{dw}{w − λ_{1}(w)} = ∫_{−1}^{0} \frac{dw}{λ_{2}(w) − w} = ∫_{0}^{1} \frac{dw}{λ_{2}(w) − w} = +∞ .

Proof. From (2.2) and (1.7), we have

2(u − λ_{1}(u)) = 2u − ψ′(u) + \sqrt{Δ(u)} = \sqrt{Δ(u)} − \frac{𝔾′′}{𝔾′} (𝔾^{−1}(u)) ,

and Δ(u) = 4 𝔾′(𝔾^{−1}(u)) + (\frac{𝔾′′}{𝔾′} (𝔾^{−1}(u)))^{2} .

We thus have

2(u − λ_{1}(u)) 𝔾′(𝔾^{−1}(u)) = [\sqrt{4(𝔾′)^{3} + (𝔾′′)^{2}} − 𝔾′′] (𝔾^{−1}(u)) .

Recall the expression

𝔾(β) = \frac{1}{tanh β} − \frac{1}{β} ,

so the function 𝔾 satisfies

lim_{β→+∞} β^{2} 𝔾′(β) = 1 ,

lim_{β→+∞} β^{3} [\sqrt{4𝔾′(β)^{3} + 𝔾′′(β)^{2}} − 𝔾′′(β)] = 2\sqrt{2} + 2 .

We also compute

lim_{β→+∞} β (1 − 𝔾(β)) = 1 , that is, lim_{u→1} 𝔾^{−1}(u) (1 − u) = 1 .

Eventually, we obtain the following asymptotic behavior:

\frac{1}{u − λ_{1}(u)} ∼ \frac{1}{(1 + \sqrt{2})(1 − u)} , when u → 1^{−} ,

and the result for the first integral follows. One shows that the three other integrals diverge in a similar way.

3.3 Global solvability of the Riemann problem

In this section, we show how to solve the Riemann problem for (2.1), that is, we solve (2.1) with an initial datum of the form:

(ρ_{0}, J_{0})(x) = \begin{cases}(ρ_{l}, J_{l}) & if x < 0, \\ (ρ_{r}, J_{r}) & if x > 0, \end{cases}

with ρ_{l}, ρ_{r} > 0, ε|J_{l}| < ρ_{l}, and ε|J_{r}| < ρ_{r}. From the general theory, see e.g. [3, 18], we know that the solution is a 1-wave (either a shock or a rarefaction) connecting (ρ_{l}, J_{l}) to an intermediate state (ρ_{m}, J_{m}), followed by a 2-wave connecting (ρ_{m}, J_{m}) to (ρ_{r}, J_{r}). In order to determine the intermediate state (ρ_{m}, J_{m}), we define the following functions for i = 1, 2:

∀ (u_{−}, u_{+}) ∈ (−1, 1) × (−1, 1) , L_{i}(u_{−}, u_{+}) := \begin{cases}X_{i}(u_{−}, u_{+}) & if u_{+} < u_{−}, \\ exp ∫_{u_{−}}^{u_{+}} \frac{dw}{λ_{i}(w) − w} & if u_{+} ≥ u_{−}. \end{cases} (3.13)

To prove some monotonicity properties for L_{1} and L_{2}, we shall need the following result:

Lemma 3.3. • One has 𝔾′′(β) > 0 for β < 0, and 𝔾′′(β) < 0 for β > 0. Moreover, one has ψ′(u) > 2u for u < 0, and ψ′(u) < 2u for u > 0.

• For all β > 0, one has

2 𝔾′(β)^{2} + 𝔾′′(β) (1 − 𝔾(β)) > 0 .

Moreover, when 0 < u_{+} < u_{−}, one has

4(ψ(u_{+}) − u_{+}^{2}) + (2u_{+} − ψ′(u_{+}))(2u_{+} − z) > 0 .

Proof. We begin with the first point of Lemma 3.3. With the definition (1.6) of 𝔾, we compute

𝔾′′(β) = \frac{1}{2[β sinh(β)]^{3}} [4 β^{3} cosh(β) + 3 sinh(β) − sinh(3β)] .

Of course, 𝔾′′ is odd. Expanding in power series (with infinite radius of convergence), we get

4 β^{3} cosh(β) + 3 sinh(β) − sinh(3β) = ∑_{p=3}^{+∞} a_{p} β^{2p+1} , a_{p} := \frac{3 − 8 p + 32 p^{3} − 3 9^{p}}{(2p + 1)!} < 0 .

This yields the first part of the Lemma. Using (1.7), we compute

ψ′(u) − 2u = \frac{𝔾′′}{𝔾′} [𝔾^{−1}(u)] .

and 𝔾^{−1}(u) has the same sign as u. We have thus proved the first point of Lemma 3.3. We now turn to the second point of Lemma 3.3. Starting from the relation

2 𝔾′(β)^{2} + 𝔾′′(β) (1 − 𝔾(β))

= \frac{1 − tanh(β)}{β^{3} sinh(β)^{2} tanh(β)} [2β^{3} + β^{2} (exp(2β) + 1) − 2β (exp(2β) − 1) + cosh(2β) − 1] ,

and expanding in power series, we obtain

2β^{3} + β^{2} (exp(2β) + 1) − 2β (exp(2β) − 1) + cosh(2β) − 1 = ∑_{p=6}^{+∞} μ_{p} β^{p} ,

with

μ_{p} = \begin{cases}\frac{2^{p−2} (p − 1)(p − 4)}{p!} & if p is even, \\ \frac{2^{p−1} (p − 5)}{(p − 1)!} & if p is odd. \end{cases}

We can deduce

∀ β > 0 , 2 𝔾′(β)^{2} + 𝔾′′(β) (1 − 𝔾(β)) > 0 .

Eventually, for 0 < u_{+} < u_{−}, we define

A := 4(ψ(u_{+}) − u_{+}^{2}) + (2u_{+} − ψ′(u_{+}))(2u_{+} − z) .

Then setting β_{+} := 𝔾^{−1}(u_{+}) > 0, and using (1.7), we obtain

A = 4 𝔾′(β_{+}) − 2 𝔾(β_{+}) \frac{𝔾′′(β_{+})}{𝔾′(β_{+})} + z \frac{𝔾′′(β_{+})}{𝔾′(β_{+})} .

Because ψ is strictly convex, we have z < ψ′(u_{−}) < 2, see Proposition 2.2, and we also have 𝔾′′(β_{+}) < 0 thanks to the first point of Lemma 3.3. Therefore, we have

A > 4 𝔾′(β_{+}) + 2 (1 − 𝔾(β_{+})) \frac{𝔾′′(β_{+})}{𝔾′(β_{+})} > 0 ,

and the proof is complete.

With the help of Lemma 3.3, we are going to show the following regularity and monotonicity properties for the functions L_{1}, L_{2}:

Proposition 3.4. The functions L_{1} and L_{2} are positive, and within C^{2}((−1, 1) × (−1, 1)). Moreover, one has

\frac{∂L_{1}}{∂u_{+}} (u_{−}, u_{+}) < 0 , \frac{∂L_{2}}{∂u_{−}} (u_{−}, u_{+}) < 0 ,

and

lim_{u_{+}→−1} L_{1}(u_{−}, u_{+}) = +∞ , lim_{u_{+}→1} L_{1}(u_{−}, u_{+}) = 0 ,

lim_{u_{−}→−1} L_{2}(u_{−}, u_{+}) = +∞ , lim_{u_{−}→1} L_{2}(u_{−}, u_{+}) = 0 .

Proof. The L_{i}’s are clearly positive, and they are also C^{∞} on either side of the diagonal {u_{+} = u_{−}}. To prove the C^{2} smoothness, one only needs to show that the partial derivatives, up to order 2, can be extended by continuity on either side of the diagonal, and do not experience a jump. This can be checked directly from the definitions (3.13). First, it is clear that the L_{i}’s are continuous. Now we compute the partial derivative of L_{1} with respect to u_{+}:

\frac{∂L_{1}}{∂u_{+}} = \frac{1}{λ_{1}(u_{+}) − u_{+}} L_{1} , if u_{+} ≥ u_{−} ,

while

2(ψ(u_{+}) − u_{+}^{2}) \frac{∂L_{1}}{∂u_{+}} + 2(ψ′(u_{+}) − 2u_{+}) L_{1} = ψ′(u_{+}) − 2u_{−} − \sqrt{δ} + (u_{−} − u_{+})\frac{∂\sqrt{δ}}{∂u_{+}} , if u_{+} < u_{−} .

The partial derivatives thus have finite limits on either side of the diagonal, and we compute

\frac{∂L_{1}}{∂u_{+}} (u_{−}, u_{−} + 0) = \frac{2}{ψ′(u_{−}) − 2u_{−} − \sqrt{Δ(u_{−})}} ,

\frac{∂L_{1}}{∂u_{+}} (u_{−}, u_{−} − 0) = \frac{ψ′(u_{−}) − 2u_{−} + \sqrt{Δ(u_{−})}}{−2(ψ(u_{−}) − u_{−}^{2})} .

These two quantities are equal thanks to (2.8). Therefore L_{1} has a continuous partial derivative with respect to u_{+}. For the partial derivative with respect to u_{−}, the calculations are similar, as for the second order partial derivatives. We shall leave the end of the calculations to the interested reader.

The regularity property of L_{1}, L_{2} is no suprise, since it is predicted by the general theory, see [3, 12, 18]. Moreover, the limits stated in the Proposition are easily deduced from Lemmas 3.1 and 3.2. The only thing left to prove is the monotonicity property of L_{1}, L_{2}. We only detail the proof for L_{1}. For u_{+} ≥ u_{−}, one has

\frac{∂L_{1}}{∂u_{+}} (u_{−}, u_{+}) = \frac{1}{λ_{1}(u_{+}) − u_{+}} L_{1}(u_{−}, u_{+}) < 0 ,

see Lemma 2.1. From now on, we thus restrict to the case u_{+} < u_{−}, and we wish to prove

∀ u_{+} < u_{−} , \frac{∂X_{1}}{∂u_{+}} (u_{−}, u_{+}) < 0 .

Using (3.3), and differentiating with respect to u_{+}, we obtain

[2(ψ(u_{+}) − u_{+}^{2})X_{1} − (ψ(u_{+}) + ψ(u_{−}) − 2u_{−}u_{+})]\frac{∂X_{1}}{∂u_{+}} = (ψ′(u_{+}) − 2u_{−})X_{1} − (ψ′(u_{+}) − 2u_{+})X_{1}^{2} .

(3.14)

Using the discriminant (3.4) of the polynomial P, the equality (3.14) also reads

\frac{(u_{−} − u_{+})\sqrt{δ}}{X_{1}} \frac{∂X_{1}}{∂u_{+}} = (1 − X_{1})(ψ′(u_{+}) − 2u_{+}) + 2(u_{+} − u_{−}) .  (3.15)

We thus want to prove that the right-hand side of (3.15) is negative when u_{+} < u_{−}. Using the first result of Lemma 3.3, we already find that the right-hand side of (3.15) is negative when u_{+} is nonpositive. We can thus assume 0 < u_{+} < u_{−}, and (ψ′(u_{+}) − 2u_{+}) < 0. In this case, the right-hand side of (3.15) is negative if and only if:

\frac{X_{1} − 1}{u_{−} − u_{+}} < \frac{2}{2u_{+} − ψ′(u_{+})} .

Using the expression of X_{1}, this is found to be equivalent to

\sqrt{δ} < \frac{4(ψ(u_{+}) − u_{+}^{2})}{2u_{+} − ψ′(u_{+})} − (z − 2u_{+}) .  (3.16)

By the second point of Lemma 3.3, the right-hand side of (3.16) is positive, and (3.16) is therefore equivalent to:

(2u_{+} − z)^{2} + 4(ψ(u_{+}) − u_{+}^{2}) < [\frac{4(ψ(u_{+}) − u_{+}^{2})}{2u_{+} − ψ′(u_{+})} − (z − 2u_{+})]^{2} .  (3.17)

Developping and making a few simplifications, we find that (3.17) is equivalent to:

4ψ(u_{+}) − ψ′(u_{+})^{2} + 2z(ψ′(u_{+}) − 2u_{+}) > 0 .

We now use the expression (3.4) of z, and we thus want to prove that

Θ(u_{−}) := (4ψ(u_{+}) − ψ′(u_{+})^{2})(u_{−} − u_{+}) + 2(ψ′(u_{+}) − 2u_{+})(ψ(u_{−}) − ψ(u_{+})) > 0 , (3.18)

as long as 0 < u_{+} < u_{−} < 1. We consider that u_{+} is fixed, and that u_{−} is a variable in the open interval (u_{+}, 1). We compute:

Θ′′(u_{−}) = 2(ψ′(u_{+}) − 2u_{+})ψ′′(u_{−}) < 0 .

Hence Θ is a strictly concave function, and we get:

∀ u_{−} ∈ (u_{+}, 1) , Θ(u_{−}) > min(Θ(u_{+}), Θ(1)) = min(0, Θ(1)) .

Eventually, using the notation β_{+} := 𝔾^{−1}(u_{+}), we compute

Θ(1) = (2 − ψ′(u_{+}))[2ψ(u_{+}) + ψ′(u_{+}) − 2u_{+} − u_{+}ψ′(u_{+})]

= \frac{2 − ψ′(u_{+})}{𝔾′(β_{+})} [2𝔾′(β_{+})^{2} + 𝔾′′(β_{+})(1 − 𝔾(β_{+}))] .

Thanks to Lemma 3.3, this last term is positive, and we can conclude that (3.18) holds. Going back to the beginning, this means that the right-hand side of (3.15) is always negative. This completes the proof.

Using the monotonicity properties of the functions L_{1}, L_{2}, we can solve the Riemann problem without any restriction on the data. Indeed, let ρ_{l}, ρ_{r} > 0 and J_{l}, J_{r} be such that ε|J_{l}| < ρ_{l} and ε|J_{r}| < ρ_{r}. To solve the Riemann problem, we need to find an intermediate state (ρ_{m}, J_{m}) such that (ρ_{l}, J_{l}) is connected to (ρ_{m}, J_{m}) by a 1-wave, and (ρ_{m}, J_{m}) is connected to (ρ_{r}, J_{r}) by a 2-wave. This means that ρ_{m}, and u_{m} = εJ_{m}/ρ_{m} must solve:

\frac{ρ_{m}}{ρ_{l}} = L_{1}(u_{l}, u_{m}) , \frac{ρ_{r}}{ρ_{m}} = L_{2}(u_{m}, u_{r}) .

Thanks to Proposition 3.4, we know that there exists a unique solution u_{m} ∈ (−1, 1) to the equation:

\frac{ρ_{r}}{ρ_{l}} = L_{1}(u_{l}, u_{m}) L_{2}(u_{m}, u_{r}) .

Consequently, the Riemann problem admits a unique solution.

When the Riemann data contain vacuum, the solution to the Riemann problem is a rarefaction wave, as follows directly from Proposition 3.3 and Lemma 3.2. For instance, if the Riemann data take the form:

(ρ, J) = \begin{cases}(0, 0) & if x < 0, \\ (ρ_{r}, J_{r}) & if x > 0, \end{cases}

with ρ_{r} > 0, ε |J_{r}| < ρ_{r}, then the solution to the Riemann problem is given by:

(ρ, J) = \begin{cases}(0, 0) & if x ≤ λ_{2}(−1)t/ε, \\ W (x/t) & if x ∈ [λ_{2}(−1)t/ε, λ_{2}(u_{r})t/ε], \\ (ρ_{r}, J_{r}) & if x ≥ λ_{2}(u_{r})t/ε. \end{cases}

Observe that the definition (2.2) yields λ_{2}(−1) = −1. When the right state for the Riemann problem is the vacuum, the corresponding solution is a 1-rarefaction wave.

We now prove the invariance of some regions for the Riemann problem. Let ρ_{0} > 0, and J_{0} ∈ ℝ be such that ε|J_{0}| < ρ_{0}. Define the following set:

I_{ρ_{0},J_{0}} := {(ρ, J) ∈ (0, ρ_{0}] × ℝ/ε|J| < ρ , Z_{i}(ρ, J) ≥ Z_{i}(ρ_{0}, J_{0}) i = 1, 2} .  (3.19)

We want to show that I_{ρ_{0},J_{0}} is an invariant region for the Riemann problem. (Recall that Z_{1}, Z_{2} are defined in Lemma 2.2). Thanks to the characterization of [10, Corollary 3.7], it is enough to prove that the set I_{ρ_{0},J_{0}} is convex (in the (ρ, J) plane). Using (2.5), and Lemma 3.2, we first get that Λ_{1} is a decreasing diffeomorphism from (−1, 1) onto ℝ, while Λ_{2} is an increasing diffeomorphism from (−1, 1) onto ℝ. Moreover, we can rewrite I_{ρ_{0},J_{0}} as

I_{ρ_{0},J_{0}} = {(ρ, J) ∈ (0, ρ_{0}) × ℝ/

Λ_{2}^{−1}(Λ_{2}(u_{0}) + ln ρ − ln ρ_{0}) ≤ \frac{εJ}{ρ} ≤ Λ_{1}^{−1}(Λ_{1}(u_{0}) + ln ρ − ln ρ_{0})\} .

Differentiating with the chain rule, one checks that the derivative of the function

J_{2} : ρ ∈ (0, ρ_{0}) ⟼ ρ Λ_{2}^{−1}(Λ_{2}(u_{0}) + ln ρ − ln ρ_{0}) ,

is given by

J′_{2}(ρ) = λ_{2} ∘ Λ_{2}^{−1}(Λ_{2}(u_{0}) + ln ρ − ln ρ_{0}) .

Hence J′_{2} is increasing, because λ_{2} and Λ_{2}^{−1} are increasing, and J_{2} is convex. In the same way, one shows that the function

J_{1} : ρ ∈ (0, ρ_{0}) ⟼ ρ Λ_{1}^{−1}(Λ_{1}(u_{0}) + ln ρ − ln ρ_{0})

is concave. Therefore, I_{ρ_{0},J_{0}} is convex, and invariant for the Riemann problem.

4 Positively invariant regions for the viscous approximation

In this section, we examine the existence of positively invariant regions for (1.8). Let ρ_{0} > 0, and consider the region I_{ρ_{0},0} defined in (3.19) (that is, for J_{0} = 0). If there was no relaxation term −J in the right-hand side of (1.8), the invariance of I_{ρ_{0},0} would simply follow from the result of [6]. For the system (1.8), we follow [1, page 66], and we thus only need to check that at each point of the boundary ∂I_{ρ_{0},0}, the vector (0, −J) points into the interior of I_{ρ_{0},0}. The same type of analysis was used in [17] for a closely related problem.

With the exception of the points (ρ_{0}, 0) and (0, 0), the boundary of I_{ρ_{0},0} is the union of two curves:

∂I_{ρ_{0},0} = {(ρ, J) ∈ (0, ρ_{0}) × ℝ/εJ = ρ Λ_{2}^{−1}( ln ρ − ln ρ_{0})}

⋃ {(ρ, J) ∈ (0, ρ_{0}) × ℝ/εJ = ρ Λ_{1}^{−1}( ln ρ − ln ρ_{0})} .

On the first of these curves, an outgoing normal vector is given by:

n := \frac{N}{|N|} , N := (\frac{1}{ε} λ_{2} ∘ Λ_{2}^{−1}(ln ρ − ln ρ_{0}), −1) ,

so that

n · (0, −J) = \frac{J}{|N|} = \frac{ρ Λ_{2}^{−1}(ln ρ − ln ρ_{0})}{|N|} < 0 ,

the last inequality being valid because Λ_{2}^{−1} is increasing and vanishes at 0.

On the second curve, one can also compute

n := \frac{N}{|N|} , N := (−\frac{1}{ε} λ_{1} ∘ Λ_{1}^{−1}(ln ρ − ln ρ_{0}), 1) ,

so that

n · (0, −J) = \frac{−J}{|N|} = \frac{−ρ Λ_{1}^{−1}(ln ρ − ln ρ_{0})}{|N|} < 0 ,

the last inequality being valid because Λ_{1}^{−1} is decreasing and vanishes at 0. Consequently, the argument of [1] shows that I_{ρ_{0},0} is positively invariant for (1.8). This completes the proof of Theorem 1.1.

The existence of an invariant region for (1.8) is a crucial step towards proving the convergence of the sequence of solutions to (1.8), because it yields L^{∞} bounds that are uniform with respect to ν. (Observe that these bounds are also uniform with respect to ε for the density ρ and for the velocity u). However, the region I_{ρ_{0},0} reaches the boundary of the domain of hyperbolicity, which yields many technical difficulties in showing the convergence of the vanishing viscosity method. Let us conclude with the following remark on the entropies for (1.1):

Remark 4.1. It is rather easy to compute the entropies for the system (1.1). Indeed, choosing the variables (ρ, u) instead of the conservative variables (ρ, J), the equation for the entropies is found to be:

∂_{ρρ}η = \frac{ψ(u) − u^{2}}{ρ^{2}} ∂_{uu}η + \frac{2u − ψ′(u)}{ρ} ∂_{ρu}η + \frac{ψ′(u) − 2u}{ρ^{2}} ∂_{u}η .  (4.1)

Note that when ψ(u) = 1 + u^{2}, one recovers the classical equation of the entropies for the isothermal Euler equations, see [13]. Note also that the coefficients of (4.1) involve the velocity u, which is a consequence of the non-Galilean invariance of (1.1).

The equation (4.1) can be simplified a little if one uses the variable β = 𝔾^{−1}(u) ∈ ℝ instead of u ∈ (−1, 1):

∂_{ρρ}η = \frac{1}{ρ^{2} 𝔾′(β)} ∂_{ββ}η − \frac{𝔾′′(β)}{ρ 𝔾′(β)} ∂_{ρβ}η .

Acknowledgments Research of the authors was supported by the EU financed network HYKE, HPRN-CT-2002-00282.

References

[1] C. Bardos. Introduction aux problèmes hyperboliques non linéaires. volume 1047 of Lecture Notes in Math., pages 1–74. Springer, 1984.

[2] C. Bardos, R. Santos, and R. Sentis. Diffusion approximation and computation of the critical size. Trans. Amer. Math. Soc., 284(2):617–649, 1984.

[3] A. Bressan. Hyperbolic systems of conservation laws. Oxford University Press, 2000.

[4] G.-Q. Chen, P. G. LeFloch. Compressible Euler equations with general pressure law. Arch. Ration. Mech. Anal., 153(3):221–259, 2000.

[5] G.-Q. Chen, P. G. LeFloch. Existence theory for the isentropic Euler equations. Arch. Ration. Mech. Anal., 166(1):81–98, 2003.

[6] K. N. Chueh, C. C. Conley, and J. A. Smoller. Positively invariant regions for systems of nonlinear diffusion equations. Indiana Univ. Math. J., 26(2):373–392, 1977.

[7] J.-F. Coulombel, F. Golse, and T. Goudon. Diffusion approximation and entropy-based moment closure for kinetic equations. Asymptotic Analysis, 45(1,2):1–39, 2005.

[8] X. Ding, G.-Q. Chen, and P. Luo. Convergence of the fractional step Lax-Friedrichs scheme and Godunov scheme for the isentropic system of gas dynamics. Comm. Math. Phys., 121(1):63–84, 1989.

[9] R. J. DiPerna. Convergence of the viscosity method for isentropic gas dynamics. Comm. Math. Phys., 91(1):1–30, 1983.

[10] D. Hoff. Invariant regions for systems of conservation laws. Trans. Amer. Math. Soc., 289(2):591–610, 1985.

[11] H. Huang, Z. Wang. Convergence of viscosity solutions for isothermal gas dynamics. SIAM J. Math. Anal., 34(3):595–610, 2002.

[12] P. D. Lax. Hyperbolic systems of conservation laws. II. Comm. Pure Appl. Math., 10:537– 566, 1957.

[13] P. G. LeFloch, V. Shelukhin. Symmetries and global solvability of the isothermal gas dynamics equations. Arch. Ration. Mech. Anal., 175(3):389–430, 2005.

[14] C. D. Levermore. Moment closure hierarchies for kinetic theories. J. Statist. Phys., 83(5-6):1021–1065, 1996.

[15] P.-L. Lions, B. Perthame, and P. E. Souganidis. Existence and stability of entropy solutions for the hyperbolic systems of isentropic gas dynamics in Eulerian and Lagrangian coordinates. Comm. Pure Appl. Math., 49(6):599–638, 1996.

[16] P.-L. Lions, B. Perthame, and E. Tadmor. Kinetic formulation of the isentropic gas dynamics and p-systems. Comm. Math. Phys., 163(2):415–431, 1994.

[17] P. Marcati, A. Milani. The one-dimensional Darcy’s law as the limit of a compressible Euler flow. J. Differential Equations, 84(1):129–147, 1990.

[18] D. Serre. Systems of conservation laws. 1. Cambridge University Press, 1999.