# Analysis of large amplitude shock profiles for non-equilibrium radiative hydrodynamics: formation of Zeldovich spikes

Jean-François Coulombel, Thierry Goudon, Pauline Lafitte, Chunjin Lin

► To cite this version:

Jean-François Coulombel, Thierry Goudon, Pauline Lafitte, Chunjin Lin. Analysis of large amplitude shock profiles for non-equilibrium radiative hydrodynamics: formation of Zeldovich spikes. Shock Waves, 2012, 22 (3), pp.181-197. 10.1007/s00193-012-0368-9 . hal-00593491

HAL Id: hal-00593491

https://hal.science/hal-00593491v1

Submitted on 16 May 2011

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
ANALYSIS OF LARGE AMPLITUDE SHOCK PROFILES FOR
NON-EQUILIBRIUM RADIATIVE HYDRODYNAMICS: FORMATION OF
ZELDOVICH SPIKES

JEAN-FRANÇOIS COULOMBEL¹,², THIERRY GOUDON², PAULINE LAFITTE², CHUNJIN LIN³

¹ CNRS

² Project Team SIMPAF, INRIA Lille - Nord Europe Research Centre
& Laboratoire Paul Painlevé (UMR CNRS 8524) CNRS-Université Lille 1
59655 Villeneuve d'Ascq Cedex, France

³ College of science, Hohai University, Nanjing, 210098, P.R. China

ABSTRACT. We consider a model for the interaction of a gas with photons. In the article [10], smooth traveling wave solutions called shock profiles have been constructed under a suitable smallness assumption between the asymptotic states. In this work, we construct piecewise smooth traveling wave solutions that connect two asymptotic states with a large jump. In particular, we give a rigorous mathematical justification to the formation of the so-called Zeldovich spike.

# 1. INTRODUCTION

This paper is concerned with the analysis of the effects of radiant heat exchanges on the shock structure of gas dynamics equations. Indeed, it is well-known that the Euler equations of gas dynamics admit discontinuous solutions, admissible discontinuities obeying the Rankine-Hugoniot conditions and the entropy inequality. These singular solutions are smoothed out by considering viscosity and/or thermal diffusion [5], producing shock profiles that join smoothly the distinct asymptotic states. Here, we wish to investigate the influence of another physical effect, dealing with models where hydrodynamics is coupled to a radiation field through energy exchanges.

The effect of the coupling with radiation on the shock structure has been explained on physical grounds in the seminal work of Zeldovich and Raizer [18, Chapter VII, §4]. We can summarize the physical picture as follows, referring to Figure 1 below for a typical density and temperature profile. At the discontinuity, radiation heats the colder fluid, creating the so-called precursor: a zone in front of the shock where the temperature increases smoothly. Owing to hydrodynamics, the upstream heating can overshoot the temperature behind the shock. The phenomena is also subject to relaxation because the radiant flux cools down the fluid in the post-shock region, equilibrating with the far end conditions. Of course, the details of these intricate mechanisms highly depend on the models chosen to describe the energy (and possibly momentum) exchanges (and it seems not reasonable to infer from the analysis of simplified models what could be the behavior of systems containing a more complex physics). The question has many important consequences, and has motivated the design of approximate solutions, see [15, 18]. These derivations often use the assumption that the upstream state is the vacuum. More elaborate solutions are designed in [13, 12], which can be used to evaluate the performances of numerical codes for radiative hydrodynamics. In this work we wish to make the following assertion, issued from [15, p. 563], mathematically rigorous "We can consider radiation exchange as a mechanism for producing partly or completely dispersed shocks. If the shock is not too strong [this mechanism gives] a completely dispersed continuous solution. With increasing shock strength, a regime is reached in which a continuous solution is not possible; instead there is a temperature discontinuity at the front, followed by a significant temperature overshoot." Actually, the first part of this claim, which is concerned with shocks of small amplitude has been investigated in [10]. By using reasoning introduced in [7] (see also [17]), it is proved that the Euler

<footnote n="1991">1991 Mathematics Subject Classification. Primary: 35L67; Secondary: 34C37.</footnote>

<footnote>Key words and phrases. Radiating gases, traveling waves, Zeldovich spike.</footnote>
system coupled to an elliptic equation for the radiation admits continuous shock profiles, when the strength of the shock is below a certain threshold, that depends on the equation of state. For this case, the temperature profile is monotone and does not exhibit any overshoot behind the shock. Moreover, the smaller the shock, the smoother the solution. The stability of these profiles is investigated in [11]. Further analysis of the problem can be found in [8, 16]. Here we address the question of the behavior of the profile when the amplitude of the shock becomes large: we shall show the occurrence of discontinuous and non monotone temperature profiles. Owing to the design of specific shock-capturing methods, numerical evidence of both smooth and non-monotone profiles in radiative hydrodynamics is discussed in [3].

To be more specific, we consider the following model:

$$ egin{cases} 
\partial_t \rho + \partial_x (\rho u) = 0, \\
\partial_t (\rho u) + \partial_x (\rho u^2 + p) = 0, \\
\partial_t (\rho E) + \partial_x (\rho E u + p u) = \sigma_a (n - \sigma \theta^4), \\
-\frac{1}{\sigma_s} \partial_{xx} n = -\sigma_a (n - \sigma \theta^4). 
\end{cases} 
\quad (1) $$

The material is described as a single-species fluid: the evolution of the density $\rho$, the velocity $u$ and the total energy $E$ is governed by the non-relativistic Euler system. The coupling with the radiation field, $n$ being the specific intensity of radiation, is treated with the non equilibrium diffusion approximation. Here and below $\sigma$ stands for the Stefan-Boltzmann constant, and $\sigma_a, \sigma_s$ are constants that measure the influence of absorption and scattering phenomena. With $\theta$ the temperature of the fluid, the coupling term $\sigma_a (n - \sigma \theta^4)$ describes the transfer of energy between the radiation and the material, due to emission and absorption processes. This system, which already appears in Zeldovich's discussion [18, sp. Eq. (7.40)-(7.43)], is certainly the simplest non equilibrium radiative hydrodynamics model. The scaling issues and the derivation of the diffusion approximation, starting from a kinetic equation for the photons distribution function, are presented e. g. in [2, 14, 6, 10], where more complete hierarchies of equations can be found. Compared to more detailed models, we observe in particular that the radiative pressure is neglected in (1), together with any relativistic corrections. It corresponds to the so-called "low energy density approximation" (because it corresponds to a regime where the typical radiative energy $\sigma\theta_0^4$ is small compared to the typical material $\rho_0 v_0^2$, with $v_0$ the sound speed). We point out that taking into account these additional terms is certainly harmful on the structure of the shock profiles (see Remark 1 below). The pressure $p$ of the fluid is assumed to satisfy the polytropic perfect gas law

$$ p = (\gamma - 1) \rho e = R \rho \theta, \quad (2) $$

where $e$ is the specific internal energy, $\gamma > 1$ is the adiabatic constant and $R > 0$ is the perfect gas constant. Eventually, the specific total energy $E$ in (1) is the sum of the specific internal energy and of the kinetic energy, that is

$$ E = e + \frac{u^2}{2}. $$

System (1) is supplemented with the entropy inequality

$$ \partial_t (\rho s) + \partial_x (\rho s u) \geq \frac{\sigma_a (n - \sigma \theta^4)}{\theta}, \quad (3) $$

where $s$ is the specific entropy, which is defined by the law of thermodynamics

$$ \theta \, ds = \, de + p \, d\frac{1}{\rho}. $$

Let us observe that (3) holds with an equality sign in regions where the solution is smooth. However we shall consider here some solutions to (1) that are only piecewise smooth and for which the entropy balance law (3) holds as a strict inequality. We recall that the inequality in (3) is determined by the vanishing viscosity principle as in classical gas dynamics. We refer for instance to [4] for some discussions on this topic, see [2, 14] in the context of radiating gases.
The paper is organized as follows. First, we need to introduce in Section 2 a few definitions, which allow to state the main result of this work. It justifies, for strong enough shocks, the existence of discontinuous and non monotone radiative profiles. We find it convenient to rewrite System (1) in dimensionless form. This is performed in Section 3 so that we can identify a single parameter, depending on the physical constants of the problem, which governs the threshold. We also show that the analysis of traveling wave solutions for (1) leads to consider a certain non linear system of ODEs. Therefore, the construction of shock profiles can be reduced to the construction of heteroclinic solutions of the differential system. This is studied in Section 4.

## 2. STATEMENT OF THE MAIN RESULT

We are interested here in traveling wave solutions to (1), (3) that only depend on the real variable $\xi = x - ct$, where c is a given velocity, and that connect some asymptotic states at $\pm\infty$, where material is at radiative equilibrium. We first define two classes of piecewise smooth functions with discontinuity of the first kind.

**Definition 1 (Piecewise smooth functions).** We let $E_1$ denote the set of functions $f: \mathbb{R} \to \mathbb{R}$ such that $f \in C^1(\mathbb{R} \setminus \{0\})$, and $f$ as well as its first derivative $f'$ have finite limits at $0^+$ and $0^-$. We let $E_2$ denote the set of functions $f \in E_1$ such that $f'$, which is well-defined on $\mathbb{R} \setminus \{0\}$, also belongs to $E_1$.

Admissible traveling waves for (1) are defined as follows.

**Definition 2 (Admissible traveling waves).** Let $(\rho_\ell, u_\ell, p_\ell)$, $(\rho_r, u_r, p_r) \in \mathbb{R}^3$ with $\rho_\ell, \rho_r, p_\ell, p_r > 0$ be two given states and let $c \in \mathbb{R}$ be a velocity. We shall say that $\xi \mapsto (\rho, u, p)(\xi) \in (E_2)^3$ defines an admissible traveling wave solution to (1) connecting $(\rho_\ell, u_\ell, p_\ell)$ and $(\rho_r, u_r, p_r)$ with velocity $c$ if the following conditions hold:

• $\lim_{\xi \to -\infty}(\rho, u, p)(\xi) = (\rho_\ell, u_\ell, p_\ell)$, $\lim_{\xi \to +\infty}(\rho, u, p)(\xi) = (\rho_r, u_r, p_r)$,

• $\inf_{\mathbb{R}^\pm} \rho > 0$, $\inf_{\mathbb{R}^\pm} p > 0$,

• for all $\varphi \in C^1(\mathbb{R})$ with compact support, there holds:

$$ \int_{\mathbb{R}} (u(\xi) - c) \rho(\xi) \varphi'(\xi) \, d\xi = 0, $$

$$ \int_{\mathbb{R}} ((u(\xi) - c) \rho(\xi) u(\xi) + p(\xi)) \varphi'(\xi) \, d\xi = 0, $$

$$ \int_{\mathbb{R}} ((u(\xi) - c) \rho(\xi) E(\xi) + p(\xi) u(\xi)) \varphi'(\xi) \, d\xi + \int_{\mathbb{R}} \sigma_a (n(\xi) - \sigma \theta(\xi)^4) \varphi(\xi) \, d\xi = 0, $$

where the function $n$ is given by

$$ \forall \xi \in \mathbb{R}, \quad n(\xi) := \frac{\sqrt{\sigma_a \sigma_s}}{2} \int_{\mathbb{R}} \exp(-\sqrt{\sigma_a \sigma_s} |\xi - \eta|) \sigma \theta(\eta)^4 \, d\eta, \quad (4) $$

• for all nonnegative $\varphi \in C^1(\mathbb{R})$ with compact support, there holds:

$$ \int_{\mathbb{R}} (u(\xi) - c) \rho(\xi) s(\xi) \varphi'(\xi) \, d\xi + \int_{\mathbb{R}} \frac{\sigma_a (n(\xi) - \sigma \theta(\xi)^4)}{\theta(\xi)} \varphi(\xi) \, d\xi \le 0. $$

An admissible traveling wave solution is called a shock profile if furthermore the function

$$ \begin{cases} (\rho_\ell, u_\ell, p_\ell), & \text{if } x < ct, \\ (\rho_r, u_r, p_r), & \text{if } x > ct, \end{cases} $$

is an entropic shock wave for the Euler equations.

Let us just make a few comments before stating our main result. First of all, functions in $E_2$ that have finite limits at $\pm\infty$ are necessarily bounded. The restrictions $\inf_{\mathbb{R}^\pm} \rho > 0$ and $\inf_{\mathbb{R}^\pm} p > 0$ in Definition 2 ensure that the temperature $\theta = p/(R \rho)$ is well-defined, belongs to $E_2$, and is bounded from below by some positive constant. (This property excludes the close to vacuum solutions studied
in [18, eq. (7.44) & (7.45)].) In particular, we can then define the energy density $n$ of photons as the unique solution in $S'(\mathbb{R})$ to the elliptic equation

$$n - \frac{1}{\sigma_a \sigma_s} \partial_{\xi\xi} n = \sigma \theta^4.$$

The solution $n$ to the latter equation is given by the convolution formula (4) in Definition 2. The definition of $n$ and the bound from below for $\theta$ make the entropy inequality in Definition 2 meaningful in the sense that all integrals involve only integrable functions. Eventually, if the traveling wave solution $(\rho, u, p)$ is smooth, that is when $\rho, u, p$ are continuous at 0, the entropy inequality in Definition 2 becomes an equality that follows from algebraic manipulations. This is the reason why the entropy inequality never appeared in the previous work [10] since the results in [10] only dealt with smooth traveling waves.

Let us now state our main result.

**Theorem 1.** Let $\gamma \in ]1, 5/3[$, and let $(\rho_\ell, u_\ell, p_\ell) \in \mathbb{R}^3$ with $\rho_\ell, p_\ell > 0$. Assume moreover that the following condition is satisfied:

$$\frac{\sigma_a \sigma^2}{\sigma_s R} \frac{\theta_\ell^7}{p_\ell^2} > \frac{\gamma + 1}{24 (\gamma - 1)^2}. \qquad (5)$$

Then there exists a constant $M > 0$ such that for all states $(\rho_r, u_r, p_r, c) \in \mathbb{R}^4$ satisfying the following conditions:

* the function

$$\begin{cases} (\rho_\ell, u_\ell, p_\ell), & \text{if } x < ct, \\ (\rho_r, u_r, p_r), & \text{if } x > ct, \end{cases}$$

is an entropic 1-shock wave for the Euler equations,

* $\|(\rho_r, u_r, p_r) - (\rho_\ell, u_\ell, p_\ell)\| \ge M,$

there exists a shock profile solution to (1) that connects $(\rho_\ell, u_\ell, p_\ell)$ to $(\rho_r, u_r, p_r)$ with velocity $c$, and $\rho, u, p, \theta$ are all discontinuous at 0. Moreover, the density profile $\rho$, the velocity profile $u$ and the pressure profile $p$ are monotone, while the temperature profile $\theta = \frac{p}{R\rho}$ is increasing on $\mathbb{R}^-$ and decreasing on $\mathbb{R}^+$. In particular, the maximum of $\theta$ on $\mathbb{R}$ is attained at $0^+$, and there holds $\theta(0^-) < \theta_r$.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Left: typical density profile. Right: typical temperature profile.</figcaption></figure>

The definition of 1-shock wave is recalled in Section 3.2 below. We refer to Figure 1 for a schematic representation of the density and temperature profiles. The peak of $\theta$ at $0^+$ is the so-called Zeldovich spike [18]. Note also that the maximal temperature of the precursor $\theta(0^-)$ always remains smaller than $\theta_r$. According to the terminology introduced in [18], the profiles we are obtaining are thus always subcritical. We believe that the restrictions on $\gamma$ and (5) in Theorem 1 are probably not optimal, but they give a good control of the behavior of a solution to a system of ordinary differential equations. Meanwhile, it is not clear that the result of Theorem 1 should hold for all value of $\gamma$.
**Remark 1.** For extremely strong shocks the model (1) becomes questionable because the underlying assumptions of the low energy density approximation do not apply to such a situation. It is likely that taking into account the radiative pressure in the momentum equation, its influence in the energy equation, as well as relativistic corrections, can restore the continuity of the profile. The analysis of the system which incorporates those terms is beyond the scope of the present work; instead we refer to the arguments developed in [15, p. 563 & pp. 579 ff.] and [18, Chapter VII, §18], and more recently to the discussion in [12].

### 3. PRELIMINARY REDUCTIONS

**3.1. The equations in dimensionless form.** In order to simplify the calculations below, we first reduce Equations (1) to a dimensionless form. To do so, we consider a fixed density $ρ_0 > 0$ and a fixed pressure $p_0 > 0$. We also introduce the constants

$$u_0 := \sqrt{\frac{p_0}{\rho_0}}, \quad \theta_0 := \frac{p_0}{R \rho_0}, \quad n_0 := \sigma \theta_0^4.$$

The new unknowns are

$$\tilde{\rho} := \frac{\rho}{\rho_0}, \quad \tilde{u} := \frac{u}{u_0}, \quad \tilde{p} := \frac{p}{p_0}, \quad \tilde{\theta} := \frac{\theta}{\theta_0}, \quad \tilde{n} := \frac{n}{n_0}.$$

Observing that the third equation in (1) reads

$$\partial_t \left( \frac{p}{\gamma - 1} + \frac{1}{2} \rho u^2 \right) + \partial_x \left( \frac{\gamma p u}{\gamma - 1} + \frac{1}{2} \rho u^3 \right) = \sigma_a (n - \sigma \theta^4),$$

we find that (1) equivalently reads

$$\left\{
\begin{aligned}
& \partial_t \tilde{\rho} + u_0 \partial_x (\tilde{\rho} \tilde{u}) = 0, \\
& \partial_t (\tilde{\rho} \tilde{u}) + u_0 \partial_x (\tilde{\rho} \tilde{u}^2 + \tilde{p}) = 0, \\
& \partial_t \left( \frac{\tilde{p}}{\gamma - 1} + \frac{1}{2} \tilde{\rho} \tilde{u}^2 \right) + u_0 \partial_x \left( \frac{\gamma \tilde{p} \tilde{u}}{\gamma - 1} + \frac{1}{2} \tilde{\rho} \tilde{u}^3 \right) = \frac{\sigma_a \sigma \theta_0^4}{p_0} (\tilde{n} - \tilde{\theta}^4), \\
& - \frac{1}{\sigma_a \sigma_s} \partial_{xx} \tilde{n} = -(\tilde{n} - \tilde{\theta}^4).
\end{aligned}
\right.
\quad (6)$$

We now introduce the characteristic length and time

$$x_0 := \frac{1}{\sqrt{\sigma_a \sigma_s}}, \quad t_0 := \frac{x_0}{u_0}.$$

In the rescaled space and time variables

$$y := \frac{x}{x_0}, \quad \tau := \frac{t}{t_0},$$

the system (6) reads

$$\left\{
\begin{aligned}
& \partial_{\tau} \tilde{\rho} + \partial_{y}(\tilde{\rho} \tilde{u}) = 0, \\
& \partial_{\tau}(\tilde{\rho} \tilde{u}) + \partial_{y}(\tilde{\rho} \tilde{u}^2 + \tilde{p}) = 0, \\
& \partial_{\tau} \left( \frac{\tilde{p}}{\gamma - 1} + \frac{1}{2} \tilde{\rho} \tilde{u}^2 \right) + \partial_{y} \left( \frac{\gamma \tilde{p} \tilde{u}}{\gamma - 1} + \frac{1}{2} \tilde{\rho} \tilde{u}^3 \right) = \frac{\sigma_a \sigma \theta_0^4 t_0}{p_0} (\tilde{n} - \tilde{\theta}^4), \\
& - \partial_{yy} \tilde{n} = -(\tilde{n} - \tilde{\theta}^4).
\end{aligned}
\right.
\quad (7)$$

The new constitutive relations for the fluid are

$$\tilde{p} = (\gamma - 1) \tilde{\rho} \tilde{e} = \tilde{\rho} \tilde{\theta}, \quad \tilde{E} = \tilde{e} + \frac{\tilde{u}^2}{2}. \quad (8)$$

In other words, we have reduced to the case $R = \sigma = \sigma_a \sigma_s = 1$ in (2). There is only one remaining dimensionless parameter in (7), which we define as

$$\alpha := \left( \frac{\sigma_a \sigma \theta_0^4 t_0}{p_0} \right)^2. \quad (9)$$
The entropy inequality (3) is reduced to

$$ \partial_{\tau}(\tilde{\rho}\tilde{s}) + \partial_{y}(\tilde{\rho}\tilde{s}\tilde{u}) \geq \frac{\sigma_a \sigma \theta_0^4 t_0}{p_0} \frac{\tilde{n} - \tilde{\theta}^4}{\tilde{\theta}}, \quad (10) $$

by introducing the rescaled entropy $\tilde{s} := s/R$ that satisfies

$$ \tilde{\theta} d\tilde{s} = d\tilde{e} + \tilde{p} d\frac{1}{\tilde{\rho}}. $$

Eventually, we observe that the choice of the constants $\rho_0$ and $p_0$ is completely arbitrary. For instance, in view of Theorem 1, we can choose $\rho_0 := \rho_\ell$ and $p_0 := p_\ell$. Let us observe that the corresponding temperature $\theta_0$ coincides with $\theta_\ell$ but $u_0$ does not necessarily coincide with $u_\ell$ (in particular because $u_0$ is positive while $u_\ell$ has no reason to be positive). We adopt this convention from now on for $\rho_0$ and $p_0$.

The notion of admissible traveling wave solution for (7), (10) is entirely analogous to Definition 2. Up to rescaling the new time and space variables $\tau, y$ by the dimensionless factor $\sqrt{\alpha}$, we have the following Definition.

**Definition 3** (Admissible traveling waves for the dimensionless system). Let $(\tilde{\rho}_\ell, \tilde{u}_\ell, \tilde{p}_\ell)$, $(\tilde{\rho}_r, \tilde{u}_r, \tilde{p}_r) \in \mathbb{R}^3$ with $\tilde{\rho}_\ell, \tilde{p}_\ell, \tilde{p}_\ell, \tilde{p}_r > 0$ be two given states and let $\tilde{c} \in \mathbb{R}$. We shall say that $\xi \mapsto (\tilde{\rho}, \tilde{u}, \tilde{p})(\xi) \in (\mathbb{E}_2)^3$ defines an admissible traveling wave solution to (7) connecting $(\tilde{\rho}_\ell, \tilde{u}_\ell, \tilde{p}_\ell)$ and $(\tilde{\rho}_r, \tilde{u}_r, \tilde{p}_r)$ with velocity $\tilde{c}$ if the following conditions hold:

* $\lim_{\xi \to -\infty} (\tilde{\rho}, \tilde{u}, \tilde{p})(\xi) = (\tilde{\rho}_\ell, \tilde{u}_\ell, \tilde{p}_\ell)$, $\lim_{\xi \to +\infty} (\tilde{\rho}, \tilde{u}, \tilde{p})(\xi) = (\tilde{\rho}_r, \tilde{u}_r, \tilde{p}_r)$,

* $\inf_{\mathbb{R}^\pm} \tilde{\rho} > 0$, $\inf_{\mathbb{R}^\pm} \tilde{p} > 0$,

* for all $\varphi \in C^1(\mathbb{R})$ with compact support, there holds:

$$ \int_{\mathbb{R}} (\tilde{u}(\xi) - \tilde{c}) \tilde{\rho}(\xi) \varphi'(\xi) d\xi = 0, $$

$$ \int_{\mathbb{R}} ((\tilde{u}(\xi) - \tilde{c}) \tilde{\rho}(\xi) \tilde{u}(\xi) + \tilde{p}(\xi)) \varphi'(\xi) d\xi = 0, $$

$$ \int_{\mathbb{R}} ((\tilde{u}(\xi) - \tilde{c}) \tilde{\rho}(\xi) \tilde{E}(\xi) + \tilde{p}(\xi) \tilde{u}(\xi)) \varphi'(\xi) d\xi + \int_{\mathbb{R}} (\tilde{n}(\xi) - \tilde{\theta}(\xi)^4) \varphi(\xi) d\xi = 0, $$

where the function $\tilde{n}$ is given by

$$ \tilde{n}(\xi) := \frac{1}{2\sqrt{\alpha}} \int_{\mathbb{R}} \exp\left(-\frac{|\xi - \eta|}{\sqrt{\alpha}}\right) \tilde{\theta}(\eta)^4 d\eta, $$

* for all nonnegative $\varphi \in C^1(\mathbb{R})$ with compact support, there holds:

$$ \int_{\mathbb{R}} (\tilde{u}(\xi) - \tilde{c}) \tilde{\rho}(\xi) \tilde{s}(\xi) \varphi'(\xi) d\xi + \int_{\mathbb{R}} \frac{\tilde{n}(\xi) - \tilde{\theta}(\xi)^4}{\tilde{\theta}(\xi)} \varphi(\xi) d\xi \leq 0. $$

An admissible traveling wave solution is called a shock profile if furthermore the function

$$ \begin{cases} (\tilde{\rho}_\ell, \tilde{u}_\ell, \tilde{p}_\ell), & \text{if } y < \tilde{c}\tau, \\ (\tilde{\rho}_r, \tilde{u}_r, \tilde{p}_r), & \text{if } y > \tilde{c}\tau, \end{cases} $$

is an entropic shock wave for the Euler equations in dimensionless form.

Given a shock profile $(\tilde{\rho}, \tilde{u}, \tilde{p}) \in \mathbb{E}_2^3$ in the sense of Definition 3, the corresponding weak solution to (7) is $(\tilde{\rho}, \tilde{u}, \tilde{p})(\sqrt{\alpha}(y - \tilde{c}\tau))$.

Of course, constructing a shock profile for the dimensionless equations (7) is entirely equivalent to constructing a shock profile for the original equations (1) since the transformation from (1) to (7) can be performed in either way. With our convention for $\rho_0$ and $p_0$ above, we are therefore reduced to constructing a shock profile for the dimensionless system (7), (8) with $\tilde{\rho}_\ell = \tilde{p}_\ell = 1$. In this case, we observe that the constant $\alpha$ in (9) coincides with the parameter on the left-hand side of condition (5) in Theorem 1. We can thus equivalently rephrase Theorem 1 as follows.
**Theorem 2.** Let $\gamma \in ]1, 5/3[$, $\alpha > 0$, and let $(\tilde{\rho}_\ell, \tilde{u}_\ell, \tilde{p}_\ell) \in \mathbb{R}^3$ with $\tilde{\rho}_\ell = \tilde{p}_\ell = 1$. Assume moreover that the following condition is satisfied:

$$ \alpha > \frac{\gamma + 1}{24(\gamma - 1)^2}. \qquad (11) $$

Then there exists a constant $\tilde{M} > 0$ such that for all states $(\tilde{\rho}_r, \tilde{u}_r, \tilde{p}_r, \tilde{c}) \in \mathbb{R}^4$ satisfying

• the function

$$ \begin{cases} (\tilde{\rho}_\ell, \tilde{u}_\ell, \tilde{p}_\ell), & \text{if } y < \tilde{c}\tau, \\ (\tilde{p}_r, \tilde{u}_r, \tilde{p}_r), & \text{if } y > \tilde{c}\tau, \end{cases} $$

is an entropic 1-shock wave for the Euler equations in dimensionless form,

$$ \bullet \quad \|(\tilde{\rho}_r, \tilde{u}_r, \tilde{p}_r) - (\tilde{\rho}_\ell, \tilde{u}_\ell, \tilde{p}_\ell)\| \ge \tilde{M}, $$

there exists a shock profile solution to (7) that connects $(\tilde{\rho}_\ell, \tilde{u}_\ell, \tilde{p}_\ell)$ to $(\tilde{\rho}_r, \tilde{u}_r, \tilde{p}_r)$ with velocity $\tilde{c}$, and $\tilde{\rho}, \tilde{u}, \tilde{p}$ are all discontinuous at 0. Moreover, the density profile $\tilde{p}$, the velocity profile $\tilde{u}$ and the pressure profile $\tilde{p}$ are monotone, while the temperature profile $\tilde{\theta}$ is increasing on $\mathbb{R}^-$ and decreasing on $\mathbb{R}^+$. In particular, the maximum of $\tilde{\theta}$ on $\mathbb{R}$ is attained at $0^+$, and moreover there holds $\tilde{\theta}(0^-) < \tilde{\theta}_r$.

**3.2. Parameterization of the shock curves.** We consider the Euler equations in dimensionless form that we have just derived in Paragraph 3.1. Dropping the tildas for convenience, we thus consider the system of conservation laws

$$ \left\{ \begin{aligned} \partial_{\tau} \rho + \partial_{y}(\rho u) &= 0, \\ \partial_{\tau}(\rho u) + \partial_{y}(\rho u^{2} + p) &= 0, \\ \partial_{\tau}(\rho E) + \partial_{y}(\rho E u + p u) &= 0, \end{aligned} \right. \qquad (12) $$

with the constitutive relations

$$ p = (\gamma - 1) \rho e = \rho \theta, \quad E = e + \frac{u^2}{2}. $$

The system (12) is supplemented with the entropy inequality

$$ \partial_{\tau}(\rho s) + \partial_{y}(\rho s u) \geq 0, \qquad (13) $$

where $s$ denotes the entropy that satisfies

$$ \theta \, ds = \, dc + p \, d\frac{1}{\rho}. $$

An entropic shock wave is a piecewise smooth solution

$$ (\rho, u, p) = \begin{cases} (\rho_{\ell}, u_{\ell}, p_{\ell}), & \text{if } y < c\tau, \\ (\rho_{r}, u_{r}, p_{r}), & \text{if } y > c\tau, \end{cases} \qquad (14) $$

that satisfies (12) in the sense of distributions as well as the entropy inequality (13). The velocity $c$ of the shock should also be distinct from the characteristic velocities of the Euler equations on either side of the shock (noncharacteristic discontinuity). Following Lax [9], we shall say that (14) is a 1-shock if the inequalities

$$ u_r - \sqrt{\frac{\gamma p_r}{\rho_r}} < c < u_r, \quad c < u_\ell - \sqrt{\frac{\gamma p_\ell}{\rho_\ell}}, $$

hold, and we shall say that (14) is a 3-shock if the inequalities

$$ u_r + \sqrt{\frac{\gamma p_r}{\rho_r}} < c, \quad u_\ell < c < u_\ell + \sqrt{\frac{\gamma p_\ell}{\rho_\ell}}, $$

hold.

We recall that for polytropic perfect gases, a piecewise smooth noncharacteristic solution to (12) of the form (14) satisfies the entropy inequality (13) if and only if it satisfies either the 1-shock inequalities or the 3-shock inequalities. In other words, the entropy admissibility criterion is equivalent to Lax shock inequalities. We refer to [1, chapter 13] for a general overview of admissibility
criteria for shock waves and for further references. The following Lemma follows from the analysis of [1, chapter 13], see also [4] for similar discussions.

**Lemma 1.** Let $(\rho_{\ell}, u_{\ell}, p_{\ell}) \in \mathbb{R}^3$ with $\rho_{\ell} = p_{\ell} = 1$. Then the set of vectors $(\rho_r, u_r, p_r, c) \in \mathbb{R}^4$ such that (14) defines a 1-shock is a smooth curve that can be parameterized by $p_r \in ]1, +\infty[$. For each of these 1-shocks, we can define the quantities

$$ \mathbf{j} := u_{\ell} - c = \rho_{r} (u_{r} - c), \quad \mathbf{C} := u_{\ell} - c + \frac{1}{\mathbf{j}} = u_{r} - c + \frac{p_{r}}{\mathbf{j}}, \quad \mathbf{a} := \frac{u_{\ell} - u_{r}}{2}. $$

When the parameter $p_r$ tends to $+\infty$, the following asymptotic expansions hold:

$$ \rho_r = \frac{1}{\mu} - \frac{1 - \mu^2}{\mu^2 p_r} + O(p_r^{-2}), \quad \mu := \frac{\gamma - 1}{\gamma + 1} \in ]0, 1[, $$

$$ \mathbf{j} = \sqrt{\frac{p_r}{1 - \mu}} \left( 1 + \frac{\mu}{2p_r} + O(p_r^{-2}) \right), $$

$$ \mathbf{C} = \sqrt{\frac{p_r}{1 - \mu}} \left( 1 + \frac{2 - \mu}{2p_r} + O(p_r^{-2}) \right), $$

$$ 2\mathbf{a} = \sqrt{(1 - \mu) p_r} \left( 1 - \frac{2 + \mu}{2p_r} + O(p_r^{-2}) \right). $$

*Proof of Lemma 1.* We briefly sketch the argument (see [1, chapter 13.4] for more details). The function (14) is a weak solution to the Euler equations (12) if and only if the Rankine-Hugoniot jump conditions hold. For noncharacteristic discontinuities, the Rankine-Hugoniot jump conditions yield the well-known relation

$$ e_r - e_\ell + \frac{p_r + p_\ell}{2} \left( \frac{1}{\rho_r} - \frac{1}{\rho_\ell} \right) = 0. $$

Using the constitutive relation $p = (\gamma - 1) \rho e$, we end up with

$$ \rho_r = \rho_\ell \frac{p_r + \mu p_\ell}{\mu p_r + p_\ell}, \quad \mu := \frac{\gamma - 1}{\gamma + 1}. $$

The asymptotic expansion of $\rho_r$ follows immediately. The expansion of $\mathbf{j}$ follows from the relation<fnref n="1" />

$$ \mathbf{j} = \sqrt{\frac{p_r - p_\ell}{\rho_r - \rho_\ell}} \rho_r \rho_\ell. $$

The end of the proof of Lemma 1 follows from some elementary computations. □

The (positive) quantity $p_r - 1$ gives a measure of the shock amplitude. When $p_r$ is close to 1, $(\rho_r, u_r, p_r)$ is close to $(\rho_\ell, u_\ell, p_\ell)$ and the velocity $c$ is close to the characteristic velocity $u_\ell - \sqrt{\gamma p_\ell / \rho_\ell}$. For large shocks, that is when $p_r$ tends to $+\infty$, all the quantities $\mathbf{j}, \mathbf{C}, \mathbf{a}$ behave like $\sqrt{p_r}$, and the velocity $c$ tends to $-\infty$. We also observe that the velocity $u_r$ varies monotonically along the shock curve, and therefore the quantity $\mathbf{a}$ determines the shock in a unique way. The asymptotic expansions of Lemma 1 will be used in Section 4 to solve the traveling wave equation that we shall derive in the following paragraph.

From now on, we consider a fixed left state $(\rho_\ell, u_\ell, p_\ell) \in \mathbb{R}^3$ with $\rho_\ell = p_\ell = 1$, and a 1-shock wave of the form (14). The analysis for 3-shocks is entirely equivalent. The expression “for large enough amplitude” means “when $p_r$ is sufficiently large” or equivalently “when $\mathbf{a}$ is sufficiently large”. Since we restrict to 1-shocks, all constants $\mathbf{j}, \mathbf{C}, \mathbf{a}$ are positive.

3.3. **The traveling wave equation and the jump conditions.** To study the shock profiles, we use the constants defined in Lemma 1 above and introduce the following system of differential equations

$$ \begin{cases} V' = W, \\ V W' = -W^2 - f(V) W + \frac{1}{2\alpha} (V^2 - \mathbf{a}^2), \end{cases} \qquad (15) $$

<footnote n="1">¹We recall that the mass flux $\mathbf{j}$ across the discontinuity is positive for a 1-shock and negative for a 3-shock.</footnote>
where the function $f$ is defined as

$$f(V) := \frac{4\mu}{\mathbf{j}} \left( \frac{\mathbf{C}}{\gamma + 1} - V \right)^3 \left( V + \frac{\gamma \mathbf{C}}{\gamma + 1} \right)^3 (2V + \mu \mathbf{C}), \quad (16)$$

and $\mu$ is defined in Lemma 1. (It will be convenient from time to time to use the parameter $\mu$ instead of $\gamma$ in order to simplify some expressions.) The function $f$ depends on the shock wave (14) through the constants $\mathbf{j}$ and $\mathbf{C}$. This is in sharp contrast with the situation in [17, 7] where the system to be solved corresponds to (15) with $f \equiv 1$ and $\alpha = 1$. We shall see below that the behavior of the function $f$ induces some technical difficulties for the construction of shock profiles, and that it completely destroys the symmetry arguments used in [7]. In the meantime, it is precisely this lack of symmetry that makes the occurrence of Zeldovich spikes possible as we shall also see below. The goal of this paragraph is to show the following result.

**Proposition 1.** Let (14) define a 1-shock wave for the Euler equations in dimensionless form. Then there exists a shock profile solution to (7) connecting $(\rho_{\ell}, u_{\ell}, p_{\ell})$ to $(\rho_r, u_r, p_r)$ with velocity $c$ and with a monotone velocity profile $u$, if and only if there exists $(V, W) \in E_2 \times E_1$ that satisfy

(i) $(V, W)$ satisfies the differential equation (15) on $\mathbb{R} \setminus \{0\}$,

(ii) $W$ takes negative values,

(iii) $\lim_{-\infty} (V, W) = (\mathbf{a}, 0)$, $\lim_{+\infty} (V, W) = (-\mathbf{a}, 0)$,

(iv) $V(0^{-}) \ge 0$, $V(0^{+}) = -V(0^{-})$, and

$$\frac{\mathbf{j}}{\mu} (V(0^{+}) W(0^{+}) - V(0^{-}) W(0^{-})) = \left(\frac{\mathbf{C}}{\gamma+1} - V(0^{+})\right)^4 \left(V(0^{+}) + \frac{\gamma \mathbf{C}}{\gamma+1}\right)^4 - \left(\frac{\mathbf{C}}{\gamma+1} - V(0^{-})\right)^4 \left(V(0^{-}) + \frac{\gamma \mathbf{C}}{\gamma+1}\right)^4.$$

*Proof of Proposition 1.* Let us first assume that there exists a shock profile solution to (7) that connects $(\rho_{\ell}, u_{\ell}, p_{\ell})$ to $(\rho_r, u_r, p_r)$ with velocity $c$. Then if we define the shifted velocity $U := u - c \in E_2$, we find that the derivative (in the sense of distributions) of both $\rho U$ and $\rho U^2 + p$ vanishes. This means that these piecewise smooth functions are constant, and examining the limit at $\pm\infty$, we obtain the relations

$$\forall \xi \in \mathbb{R} \setminus \{0\}, \quad \rho(\xi) = \frac{\mathbf{j}}{U(\xi)}, \quad \theta(\xi) = (\mathbf{C} - U(\xi)) U(\xi).$$

Now we observe that $\theta \in E_2$ is bounded, so the photons energy density $n$ given in Definition 3 by the convolution formula is a $C^1$ function and satisfies

$$\forall \xi \in \mathbb{R}, \quad n'(\xi) = -\frac{1}{2\alpha} \int_{\mathbb{R}} \operatorname{sgn}(\xi - \eta) \exp\left(-\frac{|\xi - \eta|}{\sqrt{\alpha}}\right) \theta(\eta)^4 d\eta,$$

where $\operatorname{sgn}$ denotes the sign function. Furthermore, $n$ is twice continuously differentiable on $\mathbb{R} \setminus \{0\}$ and $n''$ has finite limits at $0^{\pm}$. Let us now define the radiative heat flux as

$$q(\xi) := -\alpha n'(\xi) = \frac{1}{2} \int_{\mathbb{R}} \operatorname{sgn}(\xi - \eta) \exp\left(-\frac{|\xi - \eta|}{\sqrt{\alpha}}\right) \theta(\eta)^4 d\eta. \quad (17)$$

Then $q$ is continuous over $\mathbb{R}$ (as the convolution of a function in $L^1$ by a function in $L^{\infty}$), and $q$ tends to 0 at $\pm\infty$ (change variables $\eta \to \xi - \eta$ and apply the dominated convergence Theorem). Furthermore, $q$ satisfies the following relations where, from now on, $\partial_{\xi}$ always denotes the derivative in the sense of distributions:

$$n - \theta^4 = -\partial_{\xi} q, \quad q - \alpha \partial_{\xi} \xi q = -\alpha \partial_{\xi} \theta^4.$$

We can thus rewrite the third relation in Definition 3 as

$$\int_{\mathbb{R}} (\mathbf{j} E(\xi) + p(\xi) u(\xi) + q(\xi)) \varphi'(\xi) d\xi = 0,$$

meaning that the function $\mathbf{j}E + p u + q$ is constant over $\mathbb{R}$. After a few simplifications, see [10] for similar computations, we obtain that the radiative heat flux $q$ satisfies

$$\forall \xi \in \mathbb{R} \setminus \{0\}, \quad \frac{2\mu}{\mathbf{j}} q(\xi) = (U(\xi) - U_{\ell}) (U(\xi) - U_r).$$
Let us now introduce the shifted velocity

$$V := U - \frac{U_{\ell} + U_{r}}{2}.$$

Then we derive the relations

$$\forall \xi \in \mathbb{R} \setminus \{0\}, \quad \frac{2\mu}{\mathbf{j}} q(\xi) = V(\xi)^2 - \mathbf{a}^2, \quad \theta(\xi) = \left( \frac{\mathbf{C}}{\gamma+1} - V(\xi) \right) \left( V(\xi) + \frac{\gamma \mathbf{C}}{\gamma+1} \right).$$

Since we know that $q$ is continuous at $0$, we derive the relation $V(0^+) = \pm V(0^-)$. Furthermore, if we differentiate the latter expression of $q$ twice in the sense of distributions, we obtain

$$\frac{2\mu}{\mathbf{j}} (q - \alpha \partial_{\xi\xi} q) = V^2 - \mathbf{a}^2 - 2\alpha (V'^2 + V V'') - 2\alpha (V(0^+) V'(0^+) - V(0^-) V'(0^-)) \delta,$$

where $'$ stands for the classical derivative on $\mathbb{R} \setminus \{0\}$ and $\delta$ denotes the Dirac mass at $0$. It remains to equate the latter expression with

$$\frac{2\mu}{\mathbf{j}} (q - \alpha \partial_{\xi\xi} q) = -\frac{2\mu\alpha}{\mathbf{j}} \partial_{\xi}\theta^4 = -\frac{8\mu\alpha}{\mathbf{j}} \theta^3 \theta' - \frac{2\mu\alpha}{\mathbf{j}} (\theta(0^+)^4 - \theta(0^-)^4) \delta,$$

and the result of Proposition 1 almost follows. The only thing left to prove is that $V(0^-)$ should be nonnegative while $V(0^+)$ should be nonpositive. This last fact is obtained by considering the entropy inequality in Definition 3. Integrating by parts on $\mathbb{R}^-$ and $\mathbb{R}^+$ (recall that we consider profiles with only one possible discontinuity of the first kind at $0$), we obtain the inequality $\mathbf{j}(s(0^+) - s(0^-)) \ge 0$. There are two possible cases. Either $V$ is continuous at $0$ and in this case, we can always perform a shift and set $V(0) = 0$. Otherwise, $V$ is discontinuous at $0$ and we have $s(0^+) \ge s(0^-)$. This implies that the states at $0^-$ and $0^+$ satisfy the Rankine-Hugoniot relations with the speed $c$ as well as the entropy inequality. In particular, see [1, 4], $u(0^-)$ should be larger than $u(0^+)$ and this implies $V(0^-) > 0 > V(0^+) = -V(0^-)$.

To complete the proof of Proposition 1, we now assume that $(V, W) \in \mathbb{E}_2 \times \mathbb{E}_1$ is a solution to (15) on $\mathbb{R} \setminus \{0\}$ that satisfies conditions (ii), (iii), (iv). We define the functions

$$U := V + \frac{U_r + U_{\ell}}{2}, \quad u := U + c, \quad \rho := \frac{\mathbf{j}}{U}, \quad p := \mathbf{j}(\mathbf{C} - U), \qquad (18)$$

which all belong to $\mathbb{E}_2$. The function $\rho$ is well-defined and is bounded from below by a positive constant. Indeed, $V$ is bounded from below by $-\mathbf{a}$ so $U$ is bounded from below by $u_r - c > 0$. In the same way, $p$ is bounded from below by $p_\ell > 0$ because $p$ is monotone. Let us also define the radiative heat flux $q$ by the formula

$$\forall \xi \in \mathbb{R} \setminus \{0\}, \quad \frac{2\mu}{\mathbf{j}} q(\xi) = V(\xi)^2 - \mathbf{a}^2.$$

Since $V(0^+)^2 = V(0^-)^2$, we obtain that $q \in \mathbb{E}_2$ is continuous on $\mathbb{R}$ and tends to $0$ at $\pm\infty$. Moreover, the definition of $q$ and $\theta$ as well as the jump conditions at $0$ yield the relation

$$q - \alpha \partial_{\xi\xi} q = -\alpha \partial_{\xi} \theta^4.$$

If we define the photons energy density $n$ by the convolution formula of Definition 3, we obtain that $q$ and $-\alpha \partial_{\xi} n$ are two temperate distributions that satisfy the same elliptic equation with constant coefficients. This gives $q = -\alpha \partial_{\xi} n$, and therefore $-\partial_{\xi} q = n - \theta^4$. The relations in Definition 3 that should be satisfied by an admissible traveling wave follow from the definition (18) of the functions $\rho, u, p$. (We feel free to skip the details.) $\square$

In the following Section, we construct a piecewise smooth shock profile by analyzing the system of differential equations (15). Our goal is to construct a “heteroclinic orbit” of (15) that connects $(\mathbf{a}, 0)$ to $(-\mathbf{a}, 0)$ with an admissible jump at $0$. The jump conditions are those given in Proposition 1. The analysis relies on three steps. We first construct the heteroclinic orbit on $\mathbb{R}^-$, then on $\mathbb{R}^+$, and in the end we try to match them at $0$. In a final paragraph, we justify the monotonicity properties of the velocity, density, pressure and temperature profiles.
## 4. CONSTRUCTION OF THE SHOCK PROFILE

4.1. **Construction of the heteroclinic orbit on $\mathbb{R}^-$**. The first main task is to construct a solution to (15) on $\mathbb{R}^-$ that tends to $(\mathbf{a}, 0)$ at $-\infty$. We follow the method developed in [7], which was already adapted to the system (15) in [10] in the case of weak shocks. More precisely, let us introduce the auxiliary system of ordinary differential equations:

$$ \begin{cases} V' = V W, \\ W' = -W^2 - f(V) W + \frac{1}{2\alpha} (V^2 - \mathbf{a}^2), \end{cases} \qquad (19) $$

where the function $f$ is still defined by (16). The transformation from (15) to (19) amounts to multiplying the vector field in the right-hand side of (15) by $V$. Our first main result for the modified system (19) is the following.

**Proposition 2.** Let $\gamma \in ]1, 5/3[$, and let $\alpha > 0$ satisfy (11). Then for a shock wave of large enough amplitude, the following results hold:

(i) $f(0)^2 > 2\mathbf{a}^2/\alpha$ and we thus define the quantity

$$ W_0 := \frac{-f(0) + \sqrt{f(0)^2 - 2\mathbf{a}^2/\alpha}}{2}. $$

(ii) There exists a unique $-\text{up to a shift- solution } (V_1, W_1)$ to (19) that is defined on $\mathbb{R}$ and satisfies

$$ \lim_{\xi \to -\infty} (V_1, W_1)(\xi) = (\mathbf{a}, 0), \quad \lim_{\xi \to +\infty} (V_1, W_1)(\xi) = (0, W_0), $$

where the convergence is exponential. Moreover $V_1$ is a decreasing diffeomorphism, and $-f(\mathbf{a})/2 < W_1 < 0$.

The proof of Proposition 2 is based on a precise description of the phase portrait of (19) when the shock amplitude is large. The crucial point is to obtain a factorization of the second order polynomial in $W$ that appears in (19). Let us observe right now that Point (i) in Proposition 2 is clear since, by using the asymptotic expansions of Lemma 1, we have

$$ f(0) \sim \kappa p_r^3, \quad \mathbf{a}^2 \sim \kappa p_r, $$

for large shocks. Here and from now on, $\kappa$ denotes a generic positive constant. In particular, $W_0 \sim -\kappa p_r^{-2}$. Before proving point (ii) of Proposition 2, we establish a preliminary result.

**Lemma 2.** Let $\gamma \in ]1, 5/3[$, and let $\alpha > 0$ satisfy (11). If the shock amplitude is sufficiently large, then the polynomial

$$ W^2 + f(V) W + \frac{1}{2\alpha} (\mathbf{a}^2 - V^2), $$

evaluated at $W = -f(\mathbf{a})/2$ takes negative values for all $V \in [0, \mathbf{a}]$.

**Proof of Lemma 2.** Using the asymptotic expansions of Lemma 1, we first compute the equivalent $f(\mathbf{a}) \sim 4\mu > 0$ for large shocks. In particular, $f(\mathbf{a})$ is positive for large shocks. Therefore, we wish to show the inequality:

$$ \forall V \in [0, \mathbf{a}], \quad D(V) := f(V) - \frac{f(\mathbf{a})}{2} + \frac{1}{\alpha f(\mathbf{a})} (V^2 - \mathbf{a}^2) > 0. $$

Using Lemma 1, we find $D(0) \sim \kappa p_r^3 \gg 1$ and $D(\mathbf{a}) = f(\mathbf{a})/2 \sim 2\mu$ for large shocks. We are going to show that, for sufficiently large shocks, $D$ has the following behavior: $D$ is increasing on an interval $[0, V]$ and is decreasing on the interval $[V, \mathbf{a}]$, where $V \in ]0, \mathbf{a}]$ depends on $\gamma, \alpha$ and on the shock. This fact will imply that the minimum of $D$ on the interval $[0, \mathbf{a}]$ is attained for $V = \mathbf{a}$ and this minimum is a positive value (close to $2\mu$). The behavior of $D$ follows from a lengthy, though elementary, analysis.

We introduce the rescaled variable

$$ X := \frac{\gamma + 1}{C} V \in [0, X_{\text{max}}], \quad X_{\text{max}} := \frac{\gamma + 1}{C} \mathbf{a}, $$
and the shifted variable

$$Y := X + \frac{\gamma - 1}{2} \in \left[ \frac{\gamma - 1}{2}, Y_{\max} \right], \quad Y_{\max} := X_{\max} + \frac{\gamma - 1}{2}.$$

We compute

$$D(V) = \frac{4 \mu \mathbf{C}^7}{(\gamma + 1)^7 \mathbf{j}} F(X) - \frac{f(\mathbf{a})}{2} + \frac{\mathbf{C}^2}{(\gamma + 1)^2 \alpha f(\mathbf{a})} (X^2 - X_{\max}^2) =: D_1(X), \quad (20)$$

where we use the notation

$$F(X) := (1 - X)^3 (X + \gamma)^3 (2X + \gamma - 1). \quad (21)$$

The expression of $F$ simplifies with the introduction of the shifted variable $Y$, and we obtain

$$F(X) = 2Y(\Gamma^2 - Y^2)^3 =: G(Y), \quad \Gamma := \frac{\gamma + 1}{2}.$$

Furthermore, the asymptotic expansions of Lemma 1 give

$$Y_{\max} = \Gamma - \frac{2}{p_r} + O(p_r^{-2}),$$

for large shocks. In particular, $Y_{\max} < \Gamma$ for large shocks. We are now going to prove that for sufficiently large shocks, and under the assumptions of Lemma 2, the polynomial function $D_1$ defined in (20) is increasing on an interval $[0, \underline{X}]$ and is decreasing on the interval $[\underline{X}, X_{\max}]$.

We compute

$$\frac{(\gamma + 1)^2}{2 \mathbf{C}^2} D'_1(X) = \frac{2 \mu \mathbf{C}^5}{(\gamma + 1)^5 \mathbf{j}} G'(Y) + \frac{1}{\alpha f(\mathbf{a})} X,$$

$$\frac{(\gamma + 1)^2}{2 \mathbf{C}^2} D''_1(X) = \frac{2 \mu \mathbf{C}^5}{(\gamma + 1)^5 \mathbf{j}} G''(Y) + \frac{1}{\alpha f(\mathbf{a})},$$

$$\frac{(\gamma + 1)^2}{2 \mathbf{C}^2} D'''_1(X) = \frac{2 \mu \mathbf{C}^5}{(\gamma + 1)^5 \mathbf{j}} G'''(Y),$$

and

$$G'(Y) = -2 (Y^2 - \Gamma^2)^2 (7Y^2 - \Gamma^2),$$

$$G''(Y) = 12Y(\Gamma^2 - Y^2)(7Y^2 - 3\Gamma^2),$$

$$G'''(Y) = -12(35Y^4 - 30\Gamma^2Y^2 + 3\Gamma^4).$$

Now we examine the variations of the functions $D_1'', D_1', D_1$.

• Let us first look at the sign of $D_1'''(X)$, or equivalently at the sign of $G'''$. We have

$$G'''(Y) = -420 (Y - Y_1 \Gamma) (Y - Y_2 \Gamma) (Y + Y_1 \Gamma) (Y + Y_2 \Gamma),$$

with

$$Y_1 := \left( \frac{15 - 2\sqrt{30}}{35} \right)^{1/2}, \quad Y_2 := \left( \frac{15 + 2\sqrt{30}}{35} \right)^{1/2}.$$

The factors $Y+Y_{1,2}\Gamma$ take positive values when $Y$ belongs to the interval $[(\gamma-1)/2, Y_{\max}]$. Moreover, the factors $Y - Y_{1,2}\Gamma$ are positive for $Y = Y_{\max}$ when the shock is sufficiently large (use $Y_{\max} \sim \Gamma$ and $Y_1 < Y_2 < 1$), and they are negative for $Y = (\gamma - 1)/2$ (use $\gamma \in ]1, 5/3[$ and $1/3 < Y_1 < Y_2$). It follows that for sufficiently large shocks, $D_1'''(0)$ is negative, $D_1'''(X_{\max})$ is also negative, and $D_1'''$ has exactly two simple roots $X_1 < X_2$ in the open interval $]0, X_{\max}]$. These roots are given by

$$X_1 := Y_1 \Gamma - \frac{\gamma - 1}{2}, \quad X_2 := Y_2 \Gamma - \frac{\gamma - 1}{2}.$$

• For large shocks, the sign of $D_1''(0)$ is the same as the sign of $G''((\gamma - 1)/2)$ because the factor $\mathbf{C}^5/\mathbf{j}$ tends to $+\infty$ while $f(\mathbf{a})$ tends to $4\mu$. For $\gamma \in ]1, 5/3[$, we obtain that $D_1''(0)$ is negative. Moreover, $D_1''$ is decreasing on $[0, X_1]$, is increasing on $[X_1, X_2]$ and is decreasing on $[X_2, X_{\max}]$. Using Lemma 1, we find that $D_1''(X_{\max})$ is positive for large shocks. We deduce that $D_1''$ has a unique root $X_3 \in ]X_1, X_2[$ in the interval $[0, X_{\max}]$, that $D_1''$ takes negative values on $[0, X_3[$ and that $D_1''$ takes positive values on $]X_3, X_{\max}]$.
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Phase portrait of (19) in the region $V > 0, W < 0$. The dashed curve represents the nullcline of the second coordinate of the vector field in (19).</figcaption></figure>

• For large shocks, $D'_1(0)$ is positive (use $\gamma \in ]1, 5/3[$), while we compute

$$ \frac{(\gamma + 1)^2 f(\mathbf{a})}{2 \mathbf{C}^2} D'_1(X_{\max}) = \frac{1}{\alpha} - \frac{24 (\gamma - 1)^2}{\gamma + 1} + O(p_r^{-1}). $$

If (11) is satisfied, we obtain $D'_1(X_{\max}) < 0$ for large shocks. Moreover, $D'_1$ is decreasing on $[0, X_3]$ and is increasing on $[X_3, X_{\max}]$ so $D'_1$ has a unique root, which we call $\underline{X}$, in the interval $[0, X_{\max}]$, $D'_1$ takes positive values on $[0, \underline{X}]$ and $D'_1$ takes negative values on $[\underline{X}, X_{\max}]$. Consequently, $D_1$ is first increasing on the interval $[0, \underline{X}]$, then it is decreasing on $[\underline{X}, X_{\max}]$. As explained above, this fact is sufficient to prove that the minimum of $D$ on the interval $[0, \mathbf{a}]$ is achieved for $V = \mathbf{a}$ when the amplitude of the shock is sufficiently large. The proof of Lemma 2 is thus complete. □

Lemma 2 immediately implies the following result.

**Lemma 3.** Let $\gamma \in ]1, 5/3[$, and let $\alpha > 0$ satisfy (11). If the shock amplitude is sufficiently large, then there exist two $C^{\infty}$ functions $W_1, W_2$ on $[0, \mathbf{a}]$ that satisfy

$$ \forall V \in [0, \mathbf{a}], \quad W^2 + f(V)W + \frac{1}{2\alpha}(\mathbf{a}^2 - V^2) = (W - W_1(V))(W - W_2(V)), \\ W_2(V) < -\frac{f(\mathbf{a})}{2} < W_1(V) \le 0. $$

In particular, $W_1$ and $W_2$ satisfy

$$ W_1(0) = W_0 \sim -\kappa p_r^{-2}, \quad W_2(0) \sim -\kappa p_r^3, \quad W_1(\mathbf{a}) = 0, \quad W_2(\mathbf{a}) = -f(\mathbf{a}), $$

where $W_0$ is defined in Proposition 2.

We refer to Figure 2 for a representation of the graphs of $W_1$ and $W_2$.

**Proof of Lemma 3.** Lemma 2 directly implies that for all $V \in [0, \mathbf{a}]$, the discriminant of the polynomial (in $W$)

$$ W^2 + f(V)W + \frac{1}{2\alpha}(\mathbf{a}^2 - V^2), $$

is positive. The existence of the roots $W_{1,2}(V)$ as well as the smooth dependance on $V$ follows. The other properties are immediate (use Lemma 1 for the asymptotic expansions). □

Let us now turn to the proof of Proposition 2.

**Proof of Proposition 2.** Let us consider the differential equation (19). The point $(\mathbf{a}, 0)$ is a saddle point because the Jacobian matrix of the vector field at this point equals

$$ \begin{pmatrix} 0 & \mathbf{a} \\ \mathbf{a}/\alpha & -f(\mathbf{a}) \end{pmatrix}. $$
The latter matrix has a unique positive eigenvalue $\lambda_+$ (the other eigenvalue is negative), that is given by

$$ \lambda_+ := \frac{-f(\mathbf{a}) + \sqrt{f(\mathbf{a})^2 + 4\mathbf{a}^2/\alpha}}{2}. $$

In particular, $\lambda_+ \sim \mathbf{a}/\sqrt{\alpha} \gg 1$ for large shocks. An eigenvector associated with $\lambda_+$ is $e_+ := (1, \lambda_+/\mathbf{a})$, so we find $e_+ \sim (1, 1/\sqrt{\alpha})$ for large shocks. Applying standard results of ODE theory, we know that there exists a unique -up to a shift- maximal solution $(V_1, W_1)$ to (19) that is defined on an interval of the form $]-\infty, \xi_1[$, and such that when $\xi$ tends to $-\infty$, the orbit $(V_1, W_1)$ is tangent to the half-line $(\mathbf{a}, 0) + \mathbb{R}^- e_+$. In particular, for all $\xi$ in a neighborhood of $-\infty$, we have

$$ V_1(\xi) < \mathbf{a}, \quad W_1(\xi) < 0. $$

Let us show that the solution $(V_1, W_1)$ satisfies the properties of Proposition 2.

• We first prove that the orbit $(V_1, W_1)$ remains within the open rectangle $[0, \mathbf{a}] \times ]-f(\mathbf{a})/2, 0[$. This rectangle is represented in Figure 2. Let $K$ denote the closed rectangle $[0, \mathbf{a}] \times [-f(\mathbf{a})/2, 0]$. We decompose the boundary of $K$ into three parts

$$ \partial K = \{(\mathbf{a}, 0)\} \cup \{0\} \times [-f(\mathbf{a})/2, 0] \cup \mathcal{K}. $$

Of course, the orbit $(V_1, W_1)$ cannot pass through the stationary point $(\mathbf{a}, 0)$. Moreover, the vertical axis $\{0\} \times \mathbb{R}$ is invariant by the flow of (19). Indeed, the restriction of the flow of (19) to the vertical axis coincides with the flow of the scalar ODE

$$ W' = -W^2 - f(0)W - \frac{\mathbf{a}^2}{2\alpha}. \qquad (22) $$

Since $(V_1, W_1)$ tends to $(\mathbf{a}, 0)$ at $-\infty$, the orbit $(V_1, W_1)$ cannot pass through any point of the segment $\{0\} \times [-f(\mathbf{a})/2, 0]$. If $(V_1, W_1)$ passes through a point of $\partial K$, then this point belongs to the region $\mathcal{K}$ in which the vector field of (19) points toward the interior of $K$. We can conclude that $(V_1, W_1)$ does not pass through any point of $\partial K$, and therefore that it belongs to the interior of $K$.

• Since the orbit $(V_1, W_1)$ remains within the compact subset $K$ of $\mathbb{R}^2$, the solution is globally defined. Moreover, we have $V'_1 < 0$ so $V_1$ has a finite limit $\bar{V} < \mathbf{a}$ as $\xi$ tends to $+\infty$. The function $W'_1$ may change sign, depending on the position of the orbit with respect to the graph of $W_1$, so we can not directly conclude that $W_1$ has a limit at $+\infty$. A detailed analysis is necessary to show that $(V_1, W_1)$ converges to $(0, W_0)$ at $+\infty$.

The $\omega$-limit set of $(V_1, W_1)$ is a nonempty compact connected subset of the rectangle $K$. Since $V_1$ has a finite limit $\bar{V}$ at $+\infty$, the $\omega$-limit set of $(V_1, W_1)$ is of the form $\{\bar{V}\} \times I$ where $I$ is a closed subinterval of $[-f(\mathbf{a})/2, 0]$. Moreover, we know that the $\omega$-limit set $\{\bar{V}\} \times I$ is globally invariant by the flow of (19), and $\bar{V} \in [0, \mathbf{a}]$.

Taking a closer look at the vector field of (19), we can conclude first that $\bar{V}$ equals 0. Indeed, if $\bar{V}$ is positive, then $I$ cannot be restricted to one single point because, in such a case, the $\omega$-limit set of $(V_1, W_1)$ would be a stationary point, but (19) has no stationary point with $\bar{V} > 0$ as a first coordinate. Since we know that $I$ contains at least two points, there exists $\bar{W} \in I$ with $\bar{W} < 0$. However, the solution to (19) passing through $(\bar{V}, \bar{W})$ does not satisfy $V \equiv \bar{V}$ so $\{\bar{V}\} \times I$ is not invariant by the flow. We are led to a contradiction, hence $\bar{V}$ equals 0.

It remains to show that $I$ is the singleton $\{W_0\}$. On the vertical axis, the flow of (19) coincides with the flow of the ODE (22). Hence, $I \subset [-f(\mathbf{a})/2, 0]$ is globally invariant by the flow of (22). We use Lemma 3, and observe that the two stationary points of (22) are $W_0 \sim -\kappa p_r^{-2} \in [-f(\mathbf{a})/2, 0]$, and $W_2(0) \sim -\kappa p_r^3 \in [-f(\mathbf{a})/2, 0]$. (We recall that $f(\mathbf{a})$ tends to $4\mu$ for large shocks.) Consequently, the only compact subinterval of $[-f(\mathbf{a})/2, 0]$ that is globally invariant by the flow of (22) is $\{W_0\}$. We have thus proved that the $\omega$-limit set of $(V_1, W_1)$ is $\{(0, W_0)\}$. In other words, $(V_1, W_1)(\xi)$ tends to $(0, W_0)$ as $\xi$ tends to $+\infty$. A direct calculation shows that $(0, W_0)$ is an attractive point (the eigenvalues of the Jacobian matrix are $W_0$ and $-2W_0 - f(0)$) so the convergence is exponential. $\square$

Following the method of [7], Proposition 2 gives the main result of this paragraph.
**Proposition 3.** Let $\gamma \in ]1, 5/3[, \text{ and let } \alpha > 0$ satisfy (11). Then for a shock wave of large enough amplitude, there exists a unique -up to a shift- solution $(V_b, W_b)$ to (15) on $]-\infty, 0[$ that satisfies

$$ \lim_{\xi \to -\infty} (V_b, W_b)(\xi) = (\mathbf{a}, 0), \quad \lim_{\xi \to 0^{-}} (V_b, W_b)(\xi) = (0, W_0). $$

Moreover $V_b$ is a decreasing diffeomorphism, and $-f(\mathbf{a})/2 < W_b < 0$.

**Proof of Proposition 3.** Using the exponential convergence of $(V_1, W_1)$ towards $(0, W_0)$ at $+\infty$, we can introduce the function

$$ \forall \xi \in \mathbb{R}, \quad \Xi_1(\xi) := - \int_{\xi}^{+\infty} V_1(\eta) \, d\eta < 0. $$

Then $\Xi_1$ is an increasing $C^\infty$ diffeomorphism from $\mathbb{R}$ to $]-\infty, 0[$, and the solution $(V_b, W_b)$ to (15) is defined as $(V_1, W_1) \circ \Xi_1^{-1}$. $\square$

**4.2. Construction of the orbit on $\mathbb{R}^+$.** In this paragraph we are going to construct a solution to (15) on $\mathbb{R}^+$ that tends to $(-\mathbf{a}, 0)$ at $+\infty$. Similarly to what we have done in the preceding paragraph, we are first going to construct a solution to the auxiliary system (19) with suitable properties. Then a change of variables will yield the solution to (15). Our second main result for the auxiliary system (19) reads as follows.

**Proposition 4.** Let $\gamma \in ]1, 5/3[, \text{ and let } \alpha > 0$ satisfy (11). Then for a shock wave of large enough amplitude, there exists a unique -up to a shift- solution $(V_2, W_2)$ to (19) that is defined on an interval of the form $]-\infty, \xi_2[$ ($\xi_2$ possibly equals $+\infty$) and that satisfies

$$ \lim_{\xi \to -\infty} (V_2, W_2)(\xi) = (-\mathbf{a}, 0), \quad \lim_{\xi \to \xi_2} (V_2, W_2)(\xi) = (V_{\max}, -\infty), $$

for some $V_{\max} \in ]-\mathbf{a}, 0]$. Moreover $V_2$ is an increasing diffeomorphism, $W_2$ is decreasing, and if we consider $W_2$ as a function of $V \in ]-\mathbf{a}, V_{\max}$[, the inequality (28) below holds.

The proof of Proposition 4 is based on an accurate description of the vector field of (19) in the quarter plane $\{V < 0, W < 0\}$. We begin with two preliminary results.

**Lemma 4.** Let $\gamma \in ]1, 5/3[$ and let $\alpha > 0$. Then if the shock amplitude is sufficiently large, and if $W_0$ is defined as in Proposition 2, there holds

$$ \forall (V, W) \in ]-\mu \mathbf{C}/2, 0[ \times ]-\infty, W_1(0)[, \quad -f'(V)W + \frac{V}{\alpha} > 0. \quad (23) $$

**Proof of Lemma 4.** The proof relies on simple though lengthy computations. Recalling the definition (21) of the function $F$, we compute

$$ \begin{aligned} f'(V) &= \frac{4 \mu \mathbf{C}^6}{\mathbf{j}(\gamma + 1)^6} F'(X), \\ F'(X) &= -14 (X - 1)^2 (X + \gamma)^2 \left( X^2 + (\gamma - 1)X + \frac{3\gamma^2 - 8\gamma + 3}{14} \right), \\ F''(X) &= 42 (1 - X)(X + \gamma)(2X + \gamma - 1) \left\{ \left( X + \frac{\gamma - 1}{2} \right)^2 - \frac{3(\gamma + 1)^2}{28} \right\}. \end{aligned} $$

Using $\gamma \in ]1, 5/3[$, we obtain $F'' \le 0$ on the interval $[-(\gamma - 1)/2, 0]$ so $F' \ge F'(0)$ on this interval. When $V$ belongs to $]-\mu \mathbf{C}/2, 0[$, there holds $X \in ]-(\gamma - 1)/2, 0[$, so we already have the lower bound

$$ f'(V) \ge \frac{4 \mu \mathbf{C}^6 \gamma^2}{\mathbf{j}(\gamma + 1)^6} (8\gamma - 3 - 3\gamma^2) > 0. $$

Using the latter inequality, we can conclude that in the region $V \in ]-\mu \mathbf{C}/2, 0[$, $W < W_1(0)$, there holds

$$ -f'(V)W + \frac{V}{\alpha} > -\frac{4\mu \mathbf{C}^6 W_1(0) \gamma^2}{\mathbf{j}(\gamma+1)^6} (8\gamma - 3 - 3\gamma^2) - \frac{\mu \mathbf{C}}{2\alpha}, \quad (24) $$
and it only remains to prove that the right hand side of (24) is positive for sufficiently large shocks. Equivalently, it is sufficient to prove

$$ \frac{8 \gamma^2 \mathbf{C}^5 |W_0| \alpha}{\mathbf{j}(\gamma+1)^6} (8\gamma - 3 - 3\gamma^2) > 1, $$

for large shocks. We use the definition of $W_0$ and Lemma 1 to obtain

$$ \frac{8 \gamma^2}{(\gamma + 1)^6} |W_0| \alpha \sim \frac{\mathbf{j} \mathbf{a}^2}{\gamma \mu^2 \mathbf{C}^7}, $$

and consequently

$$ \frac{8 \gamma^2 \mathbf{C}^5 |W_0| \alpha}{\mathbf{j}(\gamma+1)^6} (8\gamma - 3 - 3\gamma^2) \sim \frac{(1-\mu)^2}{4\gamma \mu^2} (8\gamma - 3 - 3\gamma^2), $$

for large shocks. It is now a simple exercise to check that this quantity is larger than 1 for $\gamma \in ]1, 5/3[$ (and $\mu \in ]0, 1/4[$). $\square$

**Lemma 5.** Let $\gamma \in ]1, 5/3[$ and let $\alpha > 0$. If the shock amplitude is sufficiently large, we have the following results:

(i) There exists a unique pair of real numbers $(\underline{V}, \overline{V})$ satisfying $-\mathbf{a} < \underline{V} < -\mu \mathbf{C}/2 < \overline{V} < 0$ and

$$ \Delta(V) := f(V)^2 + \frac{2}{\alpha}(V^2 - \mathbf{a}^2) \begin{cases} > 0, & \text{if } V \in [-\mathbf{a}, \underline{V}[ \cup [\overline{V}, 0], \\ < 0, & \text{if } V \in [\underline{V}, \overline{V}]. \end{cases} $$

(ii) For any $(V, W) \in [-\mathbf{a}, -\mu \mathbf{C}/2] \times \mathbb{R}^-$, there holds

$$ W^2 + f(V)W + \frac{\mathbf{a}^2 - V^2}{2\alpha} > 0. \qquad (25) $$

(iii) There exist two continuous functions $W_1, W_2$ on $[\overline{V}, 0]$, that are of class $C^\infty$ on $[\overline{V}, 0]$ and that satisfy

$$ \forall V \in [\overline{V}, 0], \quad W^2 + f(V)W + \frac{1}{2\alpha}(\mathbf{a}^2 - V^2) = (W - W_1(V))(W - W_2(V)), \\ W_2(V) < W_1(V) < 0. $$

(iv) For all $(V, W) \in [-\mu \mathbf{C}/2, 0[ \times ]-\infty, W_2(0)[$, the inequality (25) holds.

*Proof of Lemma 5.* • Let us begin with the proof of point (i). As in the proof of Lemma 2, we introduce the rescaled variable $X := (\gamma + 1) V/\mathbf{C}$, and get

$$ \Delta(V) = \frac{16 \mu^2 \mathbf{C}^{14}}{(\gamma + 1)^{14} \mathbf{j}^2} \left[ F(X)^2 + \frac{(\gamma + 1)^{12} \mathbf{j}^2}{8 \mu^2 \mathbf{C}^{12} \alpha} (X^2 - X_{\max}^2) \right] =: \frac{16 \mu^2 \mathbf{C}^{14}}{(\gamma + 1)^{14} \mathbf{j}^2} H(X), \quad (26) $$

where the polynomial function $F$, and the parameter $X_{\max}$ were already defined in the proof of Lemma 2:

$$ F(X) = (1 - X)^3 (X + \gamma)^3 (2X + \gamma - 1), \quad X_{\max} = 1 - \frac{2}{p_r} + O(p_r^{-2}). $$

Observe that for sufficiently large shocks, there holds

$$ \frac{(\gamma + 1)^{12} \mathbf{j}^2}{8 \mu^2 \mathbf{C}^{12} \alpha} \sim \kappa p_r^{-5} \ll 1. $$

Furthermore we know that the holomorphic function $F^2$ has exactly one double root in the closed disk

$$ \{X \in \mathbb{C} / (\operatorname{Im} X)^2 + \operatorname{Re} X (\operatorname{Re} X + 1) \le 0\}. $$

Applying Rouché's Theorem, we find that the function $H$, defined in (26), has two (complex) roots inside the same closed disk provided that the shock amplitude is sufficiently large. Since $H$ is real valued when $X$ is real and using

$$ H(-X_{\max}) > 0, \quad H\left(-\frac{\gamma - 1}{2}\right) < 0, \quad H(0) > 0, $$
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Phase portrait of (19) in the domain $V < 0$, $W < 0$. The dashed curves represent the nullclines of the second coordinate of the vector field.</figcaption></figure>

we can conclude that $H$ has exactly two real distinct (and simple) roots in the interval $[-X_{\max}, 0]$, and these roots are located in the intervals $]-X_{\max}, -(\gamma-1)/2[$ and $]-(\gamma-1)/2, 0[$. Finally, going back to the variable $V$, we get property (i).

• We now prove the property (ii). According to (i), we already know that the second order polynomial equation in $W$

$$W^2 + f(V)W + \frac{1}{2\alpha}(\mathbf{a}^2 - V^2) = 0, \qquad (27)$$

has no real root for $V \in ]\bar{V}, -\mu C/2]$. In particular, property (ii) holds for $V$ in the interval $]\bar{V}, -\mu C/2]$, so we now assume $V \in ] -\mathbf{a}, \bar{V}]$. In that case, (i) shows that Equation (27) has two real roots, and it is not difficult to see that these roots are positive. Consequently, for any $(V, W) \in ] -\mathbf{a}, \bar{V}] \times \mathbb{R}^-$, the inequality (25) holds.

• The property (iii) is also a direct consequence of (i). For $V \in [\bar{V}, 0]$, Equation (27) has two real roots. These roots are simple for $V \neq \bar{V}$ and they are necessarily negative. We refer to Figure 3 for a schematic picture of the graphs of $W_1$ and $W_2$. The functions $W_1, W_2$ are nothing but the extension to $[\bar{V}, 0]$ of the same functions that were previously defined on $[0, \mathbf{a}]$, see Proposition 2.

• For $V \in [-\mu C/2, \bar{V}]$, property (i) shows that Equation (27) has no real root. In particular, (25) holds for all $W \in \mathbb{R}$ and not only for $W < W_2(0)$. Let us now consider the case $V \in [\bar{V}, 0]$. Our goal is to show that the roots $W_{1,2}(V)$ are both larger than $W_2(0)$, see Figure 3. The proof relies on Lemma 4 and splits in two steps.

**Step 1:** the function $W_1$ is increasing on $[\bar{V}, 0]$. Indeed, for $V > \bar{V}$, we have $\sqrt{\Delta(V)} W_1'(V) = -f'(\bar{V}) W_1(V) + V/\alpha$. In particular, there holds $W_1'(0) > 0$ for sufficiently large shocks (use $f'(0) > 0$ and $W_0 < 0$). Let us assume that there exists $V^* \in ]\bar{V}, 0[$ such that $W_1'(V^*) \le 0$. Then there is no loss of generality in assuming that $V^*$ is the largest zero of $W_1'$ in the interval $]\bar{V}, 0[$, and $W_1'$ is therefore positive on the interval $]V^*, 0]$. Since $W_1(V^*) < W_1(0)$, we can apply Lemma 4. Inequality (23) yields $\sqrt{\Delta(V^*)} W_1'(V^*) > 0$, hence a contradiction. Consequently we get $W_1'(V) > 0$ for all $V \in ]\bar{V}, 0[$ and $W_1$ is increasing.

**Step 2:** $W_2$ is decreasing on $[\bar{V}, 0]$. Indeed, we have $W_2(V) < W_1(V) < W_1(0)$ for all $V \in ]\bar{V}, 0[$ (use step 1). We then use the relation $\sqrt{\Delta(V)} W_2'(V) = f'(V) W_2(V) - V/\alpha < 0$ where the last inequality follows from Lemma 4 again. In particular, for all $V \in [\bar{V}, 0[$, the roots $W_1(V), W_2(V)$ to Equation (27) are larger than $W_2(0)$, so (25) holds for all $W < W_2(0)$. The proof of (iv) is now complete. □

Let us now prove Proposition 4.
*Proof of Proposition 4.* We follow more or less the proof of Proposition 2. The stationary point $(-\mathbf{a}, 0)$ is a saddle point for the system (19). The Jacobian of the vector field at $(-\mathbf{a}, 0)$ equals the matrix

$$ \begin{pmatrix} 0 & -\mathbf{a} \\ -\mathbf{a}/\alpha & -f(-\mathbf{a}) \end{pmatrix}, $$

whose only positive eigenvalue (the other eigenvalue is negative) is

$$ \lambda_+ := \frac{-f(-\mathbf{a}) + \sqrt{f(-\mathbf{a})^2 + 4\mathbf{a}^2/\alpha}}{2}. $$

The correspondent eigenvector is $e_+ := (1, -\lambda_+/\mathbf{a})$. Since $f(-\mathbf{a}) \sim -\kappa p_r^3$ for large shocks, there holds $\lambda_+ \sim \kappa p_r^3$, $e_+ \sim (1, -\kappa p_r^{5/2})$. By standard ODE theory, there exists a unique -up to a shift-solution $(V_2, W_2)$ to (19) that is defined on an interval of the form $]-\infty, \xi_2[$ and such that for $\xi \to -\infty$, the orbit tends to $(-\mathbf{a}, 0)$ with the tangent line $(-\mathbf{a}, 0) + \mathbb{R}^+ e_+$. In particular, for $\xi$ in a neighborhood of $-\infty$, there holds

$$ V_2(\xi) > -\mathbf{a}, \quad W_2(\xi) < 0. $$

Using the same arguments as in the proof of Proposition 2, we can already conclude that the orbit $(V_2, W_2)$ remains within the region $]-\mathbf{a}, 0[ \times ]-\infty, 0[$. (This property is easily proved by simply looking at the behavior of the vector field on the boundary of this region.) In particular, $V_2$ is increasing.

Let us define $V_{\max} := \sup_{\xi < \xi_2} V_2(\xi)$. (This quantity is independent of the possible shift in the definition of $(V_2, W_2)$.) Since $V_2$ is a diffeomorphism from $]-\infty, \xi_2[$ to $]-\mathbf{a}, V_{\max}[$, we feel free to consider $W_2$ as a function of $V \in ]-\mathbf{a}, V_{\max}[$, the orbit $(V_2, W_2)$ being the graph of $W_2$. We obviously have $W_2(-\mathbf{a}) = 0$. The crucial point for proving Proposition 4 is to show that $W_2$ must diverge as $V$ tends to $V_{\max}$. This is proved as follows. From the system (19), we get

$$ \frac{dW_2}{dV} = - \frac{W_2^2 + f(V)W_2 + (\mathbf{a}^2 - V^2)/(2\alpha)}{V W_2} \leq - \frac{W_2 + f(V)}{V}, $$

which gives

$$ \forall V \in ]-\mathbf{a}, V_{\max}[; \quad -\frac{d(V W_2)}{dV} \leq f(V). $$

We integrate this inequality from $-\mathbf{a} + \varepsilon$ to $V \in ]-\mathbf{a}, V_{\max}[$ and then pass to the limit $\varepsilon \to 0$. The function $f$ is integrated by using the relation

$$ f(V) = -\frac{\mu}{\mathbf{j}} \frac{d}{dV} \left[ \left( \frac{\mathbf{C}}{\gamma+1} - V \right)^4 \left( \frac{\gamma \mathbf{C}}{\gamma+1} + V \right)^4 \right]. $$

We thus obtain

$$ \forall V \in ]-\mathbf{a}, V_{\max}[; $$

$$ W_2(V) \leq \frac{\mu}{\mathbf{j}V} \left[ \left( \frac{\mathbf{C}}{\gamma+1} - V \right)^4 \left( \frac{\gamma \mathbf{C}}{\gamma+1} + V \right)^4 - \left( \frac{\mathbf{C}}{\gamma+1} + \mathbf{a} \right)^4 \left( \frac{\gamma \mathbf{C}}{\gamma+1} - \mathbf{a} \right)^4 \right] =: G(V). \quad (28) $$

To complete the proof of Proposition 4, the only remaining thing to prove is $\lim_{V \to V_{\max}} W_2(V) = -\infty$. There are two possible cases according to the position of $V_{\max}$ with respect to $-\mu \mathbf{C}/2$.

**Case 1:** $V_{\max} \leq -\mu \mathbf{C}/2$. In this case the orbit $(V_2, W_2)$ remains within the region $]-\mathbf{a}, -\mu \mathbf{C}/2[ \times \mathbb{R}^-$, and (ii) shows that $W_2$ is decreasing. If $W_2$ is bounded from below, then the $\omega$-limit set of $(V_2, W_2)$ is a stationary point. However there is no stationary point of (19) with $V = V_{\max}$ and $W < 0$, so $W_2$ is not bounded from below. We thus conclude $\lim_{\xi \to \xi_2} W_2(\xi) = -\infty$.

**Case 2:** $V_{\max} > -\mu \mathbf{C}/2$. Using the inequality (28), we have $W_2(-\mu \mathbf{C}/2) \leq G(-\mu \mathbf{C}/2)$, and Lemma 1 gives

$$ G\left(-\frac{\mu \mathbf{C}}{2}\right) \sim 2(1-\mu)\left(\mu^4 - \frac{1}{256(1-\mu)^4}\right) p_r^3, $$
for large shocks. On the other hand, Lemma 1 also gives

$$W_2(0) \sim -\frac{\mu^2 (1+\mu)^3}{16} p_r^3,$$

for large shocks, so we have $W_2(-\mu C/2) \le G(-\mu C/2) < W_2(0)$. This fact relies on the inequality

$$\frac{\mu^2 (1+\mu)^3}{16} < 2(1-\mu) \left( \frac{1}{256(1-\mu)^4} - \mu^4 \right),$$

which holds for $0 < \mu < 1/4$. (This can be proved by studying the variations of these quantities with respect to $\mu$.) The proof of the latter inequality is one of the crucial points where we use the assumption $1 < \gamma < 5/3$. In particular, once the orbit $(V_2, W_2)$ passes through the point $(-\mu C/2, W_2(-\mu C/2))$, it enters the region $]-\mu C/2, 0[ \times ]-\infty, W_2(0]$ and the property (iii) shows that the orbit cannot escape from this region. In particular, $W_2$ is a decreasing function and the same argument as above shows that $W_2$ tends to $-\infty$ as $\xi$ tends to $\xi_2$. $\square$

**Proposition 5.** There exists a unique -up to a shift- solution $(V_\#(\xi), W_\#(\xi))$ to (15) that is defined on an interval of the form $]-\infty, +\infty[$ and that satisfies

$$\lim_{\xi \to +\infty} (V_\#(\xi), W_\#(\xi)) = (-\mathbf{a}, 0), \quad \lim_{\xi \to \xi_\#} (V_\#(\xi), W_\#(\xi)) = (V_{\max}, -\infty),$$

for some $V_{\max} \in ]-\mathbf{a}, 0]$. Furthermore, $V_\#(\xi)$ is a decreasing diffeomorphism and $W_\#(\xi)$ is increasing.

**Proof of Proposition 5.** The proof is entirely analogous to the proof of Proposition 3. There is no loss of generality in assuming $\xi_2 > 0$. We introduce the function

$$\forall \xi \in ]-\infty, \xi_2[$, $\Xi_2(\xi) := -\int_{\xi}^{0} V_2(\eta) \, d\eta,$$

which is a decreasing $C^\infty$ diffeomorphism from $]-\infty, \xi_2[$ to an open interval of the form $]-\infty, +\infty[$. The desired solution $(V_\#(\xi), W_\#(\xi))$ to (15) is defined as $(V_\#(\xi), W_\#(\xi)) := (V_2, W_2) \circ \Xi_2^{-1}$. $\square$

**4.3. Matching the orbits at 0.** In this paragraph, up a shift in $\xi$, we need to determine the points $(V(0^-), W(0^-))$ and $(V(0^+), W(0^+))$ of the piecewise solution to (15). According to Propositions 3 and 5, both functions $V_\#(\xi)$ and $W_\#(\xi)$ are decreasing diffeomorphisms so we feel free to consider $W_\#(\xi)$ and $W_\#(\xi)$ as functions of $V$. More precisely, $W_\#(\xi)$ can be considered as a function of $V \in [0, \mathbf{a}]$ while $W_\#(\xi)$ can be considered as a function of $V \in [-\mathbf{a}, V_{\max}]$.

Let us introduce the following two continuous functions:

$$K: V \in [-\mathbf{a}, V_{\max}] \mapsto W_\#(V) + W_\#(-V),$$

and

$$L: V \in [-\mathbf{a}, V_{\max}] \mapsto \frac{\mu}{jV} \left[ \left( \frac{\mathbf{C}}{\gamma+1} - V \right)^4 \left( \frac{\gamma \mathbf{C}}{\gamma+1} + V \right)^4 - \left( \frac{\mathbf{C}}{\gamma+1} + V \right)^4 \left( \frac{\gamma \mathbf{C}}{\gamma+1} - V \right)^4 \right].$$

Lemma 1 shows that $L(-\mathbf{a})$ is a large negative number for large shocks while $K$ vanishes at $-\mathbf{a}$, hence $(K-L)(-\mathbf{a}) > 0$. Let us first consider the case $V_{\max} \le -\mathbf{a} + \mu \mathbf{C}$ (with the assumption $1 < \gamma < 5/3$, Lemma 1 shows that $-\mathbf{a} + \mu \mathbf{C}$ is negative for large shocks). Then $L$ remains bounded for $V$ tending to $V_{\max}$ while $K$ tends to $-\infty$. Consequently, the function $K-L$ vanishes at some point, denoted $V^+$, in the interval $]-\mathbf{a}, V_{\max}[ \subset ]-\mathbf{a}, -\mathbf{a} + \mu \mathbf{C}]$.

Let us now consider the case $V_{\max} > -\mathbf{a} + \mu \mathbf{C}$. We observe that with the above definition for $L$ and with the function $G$ defined in (28), there holds $L(-\mathbf{a} + \mu \mathbf{C}) = G(-\mathbf{a} + \mu \mathbf{C})$, so we get

$$K(-\mathbf{a} + \mu \mathbf{C}) - L(-\mathbf{a} + \mu \mathbf{C}) = \underbrace{W_\#(\mathbf{a} - \mu \mathbf{C})}_{<0} + \underbrace{W_\#(-\mathbf{a} + \mu \mathbf{C}) - G(-\mathbf{a} + \mu \mathbf{C})}_{\le 0} < 0,$$

where we have used (28). We can again conclude that the function $K-L$ vanishes at some point, denoted $V^+$, in the interval $]-\mathbf{a}, -\mathbf{a} + \mu \mathbf{C}]$.
Up to using a shift in $\xi$, there is no loss of generality in assuming $V_b(0) = -V^+$, $V_\#(0) = V^+$. Then the piecewise solution to (15) is defined by

$$ \forall \xi \neq 0, \quad (V, W)(\xi) := \begin{cases} (V_b, W_b)(\xi), & \text{if } \xi < 0, \\ (V_\#, W_\#)(\xi), & \text{if } \xi > 0. \end{cases} $$

It is straightforward to check that $(V, W)$ satisfies properties (i)-(iv) of Proposition 1. To complete the proof of Theorem 2, it remains to check the monotonicity properties of the density, temperature profiles etc. This will be done in the next paragraph. We emphasize that with our construction, the limit $V(0^+)$ belongs to the interval $]-a, -a + \mu C[$.

### 4.4. Justification of the Zeldovich spike.
The goal of this last paragraph is to show the monotonicity properties for the velocity, density, pressure and temperature profiles. We recall that under the assumption $1 < \gamma < 5/3$ made in Theorem 2, there holds $0 < \mu < 1/4$ and $-a + \mu C < -\mu C/2$ for large enough shocks (use Lemma 1).

In the preceding paragraph, we have shown how to construct a solution $(V, W)$ to (15) satisfying properties (i), (ii), (iii), (iv) of Proposition 1. We have already seen in the proof of Proposition 1 that the corresponding velocity, density and pressure profiles are given by (18). In particular, $u$ is decreasing (because $V$ is decreasing) while $\rho$ and $p$ are increasing. The temperature profile is given by

$$ \theta(\xi) = \left( \frac{\mathbf{C}}{\gamma + 1} - V(\xi) \right) \left( \frac{\gamma \mathbf{C}}{\gamma + 1} + V(\xi) \right), \quad (29) $$

so we have

$$ \forall \xi \neq 0, \quad \theta'(\xi) = -(2V(\xi) + \mu \mathbf{C}) V'(\xi). $$

In particular, $\theta'(\xi)$ has the same sign as $2V(\xi) + \mu \mathbf{C}$. If $\xi < 0$, $V(\xi)$ is positive since it belongs to $]V(0^-), \mathbf{a}]$, so $\theta'(\xi)$ is positive and $\theta$ is increasing on $\mathbb{R}^-$. If $\xi > 0$, we have shown in the preceding paragraph that $V(\xi)$ belongs to the interval $]-\mathbf{a}, -\mathbf{a} + \mu \mathbf{C}]$, so $2V(\xi) + \mu \mathbf{C}$ is negative and $\theta$ is decreasing on $\mathbb{R}^+$.

To complete the proof of Theorem 2, we need to show $\theta(0^-) < \theta_r$. Using the expression (29), we have

$$ \theta(0^-) - \theta_r = (\mathbf{a} + V(0^-))(\mathbf{a} - \mu \mathbf{C} - V(0^-)), $$

and since $V(0^-) = -V(0^+) > \mathbf{a} - \mu \mathbf{C}$, we get the desired inequality. The proof of Theorem 2 is thus complete.

Similar manipulations can be used to estimate the temperature overshoot. Coming back to (29), we obtain

$$ \theta(0^+) - \theta_r = (\mathbf{a} + V(0^+))(\mathbf{a} - \mu \mathbf{C} - V(0^+)) < \mu \mathbf{C}(2\mathbf{a} - \mu \mathbf{C}), $$

since $-\mathbf{a} < V(0^+) < \mu \mathbf{C} - \mathbf{a}$. Then, bearing in mind $\frac{2}{\gamma+1} = 1 - \mu$ and $\frac{2\gamma}{\gamma+1} = 1 + \mu$, we use the asymptotic behavior detailed in Lemma 1. It yields

$$ \theta_r = \left( \frac{\mathbf{C}}{\gamma + 1} + \mathbf{a} \right) \left( \frac{\gamma \mathbf{C}}{\gamma + 1} - \mathbf{a} \right) \sim \mu p_r $$

and

$$ \mu \mathbf{C}(2\mathbf{a} - \mu \mathbf{C}) \sim \mu p_r \left(1 - \frac{\mu}{1 - \mu}\right) = \mu p_r \frac{3 - \gamma}{2}. $$

Therefore, for large shocks, $\frac{5-\gamma}{2}\theta_r$ can be considered as a bound from above of the maximal temperature $\theta(0^+)$.

## REFERENCES

[1] S. Benzoni-Gavage, D. Serre. Multidimensional hyperbolic partial differential equations. Oxford Mathematical Monographs. Oxford University Press, 2007.

[2] C. Buet, B. Després. Asymptotic analysis of fluid models for the coupling of radiation and hydrodynamics. J. Quant. Spectrosc. Radiat. Transfer, 85:385-418, 2004.

[3] J.-F. Coulombel and P. Lafitte. Computation of shock profiles in radiative hydrodynamics. CiCP, 6:1118-1136, 2009.

[4] R. Courant, K. O. Friedrichs. Supersonic flow and shock waves. Springer-Verlag, 1976. Applied Mathematical Sciences, Vol. 21.
[5] D. Gilbarg. The existence and limit behavior of the one-dimensional shock layer. Amer. J. Math., 73:256-274, 1951.

[6] P. Godillon-Lafitte and T. Goudon. A coupled model for radiative transfer: Doppler effects, equilibrium and non-equilibrium diffusion asymptotics. SIAM Multiscale Modeling and Simulation, 4(4):1245-1279, 2005.

[7] S. Kawashima, S. Nishibata. Shock waves for a model system of the radiating gas. SIAM J. Math. Anal., 30(1):95-117, 1999.

[8] C. Lattanzio, C. Mascia, and D. Serre. Shock waves for radiative hyperbolic-elliptic systems. Indiana Univ. Math. J., 56(5):2601-2640, 2007.

[9] P. D. Lax. Hyperbolic systems of conservation laws and the mathematical theory of shock waves. Society for Industrial and Applied Mathematics, 1973. Conference Board of the Mathematical Sciences Regional Conference Series in Applied Mathematics, No. 11.

[10] C. Lin, J.-F. Coulombel, and T. Goudon. Shock profiles for non-equilibrium radiating gases. Phys. D, 218(1):83-94, 2006.

[11] C. Lin, J.-F. Coulombel, and T. Goudon. Asymptotic stability of shock profiles in radiative hydrodynamics. CRAS, 345(11):625-628, 2007.

[12] R. Lowrie and J. D. Edwards. Radiative shock solutions with grey nonequilibrium diffusion. Shock Waves, 18(2):129-143, 2008.

[13] R. Lowrie and R. M. Rauenzahn. Radiative shock solutions in the equilibrium diffusion limit. Shock Waves, 16(6):445-453, 2007.

[14] R. B. Lowrie, J. E. Morel, and J. A. Hittinger. The coupling of radiation and hydrodynamics. The Astrophysical J., 521:432-450, 1999.

[15] D. Mihalas, B. W. Mihalas. Foundations of radiation hydrodynamics. Oxford University Press, 1984.

[16] T. Nguyen, R. G. Plaza, and K. Zumbrun. Stability of radiative shock profiles for hyperbolic-elliptic coupled systems. Phys. D, 239(8):428-453, 2010.

[17] S. Schochet, E. Tadmor. The regularized Chapman-Enskog expansion for scalar conservation laws. Arch. Rational Mech. Anal., 119(2):95-107, 1992.

[18] Ya. B. Zel'dovich, Yu. P. Raizer. Physics of shock waves and high-temperature hydrodynamic phenomena. Dover Publications, 2002.

E-mail address: jean-francois.coulombel@math.univ-lille1.fr

E-mail address: thierry.goudon@inria.fr

E-mail address: pauline.lafitte@math.univ-lille1.fr

E-mail address: chunjin.lin@gmail.com