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

We consider a model for the interaction of a gas with photons. In the article (Lin et al. Phys D 218:83–94, 2006), 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.


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 [19,Chapter VII,Section 4]. We can summarize the physical picture as follows, referring to Fig. 1 below for a typical density and temperature profile. At the discontinuity, radiation heats the colder fluid, creating the so-called precursor: a zone ahead of the shock where the temperature increases smoothly. Owing to hydrodynamics, the upstream heating can overshoot the temperature behind the shock. The phenomenon 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 does not seem 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 e.g. [16,19]. These derivations often use the assumption that the upstream state is the vacuum. More elaborate solutions are designed in [12,13], 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 [16, 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 already been investigated in [10]. By using reasoning introduced in [7], see also [18], it was proved in [10] that the Euler 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 asymptotic stability of these profiles with respect to small perturbations is investigated in [11]. Further analysis of the problem in an abstract framework can be found in [8,17]. 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 shockcapturing 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: (1) The material is described as a single-species fluid: the evolution of the density ρ, the velocity u and the total energy E is governed by the non-relativistic Euler system. Coupling with the radiation field, n being the specific intensity of radiation, is treated with the non-equilibrium diffusion approximation. Here and below σ stands for the Stefan-Boltzmann constant, and σ a , σ s are constants that measure the influence of absorption and scattering phenomena. With θ the temperature of the fluid, the coupling term σ a (n − σ θ 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 [19, 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,6,10,14], 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 σ θ 4 0 is small compared to the typical kinetic energy ρ 0 v 2 0 , 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 where e is the specific internal energy, γ > 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 System (1) is supplemented with the entropy inequality where s is the specific entropy, which is defined by the law of thermodynamics 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 Sect. 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 Sect. 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 Ordinary Differential Equations. Therefore, the construction of shock profiles can be reduced to the construction of heteroclinic solutions of the differential system. This is studied in Sect. 4.

Statement of the main result
We are interested here in traveling wave solutions to (1), (3) that only depend on the real variable ξ = x − c t, where c is a given velocity, and that connect some asymptotic states at ±∞, 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 : R → R such that f ∈ C 1 (R\{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 ∈ E 1 such that f , which is well defined on R\{0}, also belongs to E 1 .
Admissible traveling waves for (1) are defined as follows.
Definition 2 (Admissible traveling waves) Let (ρ , u , p ), (ρ r , u r , p r ) ∈ R 3 with ρ , ρ r , p , p r > 0 be two given states and let c ∈ R be a velocity. We shall say that ξ → (ρ, u, p)(ξ ) ∈ (E 2 ) 3 defines an admissible traveling wave solution to (1) connecting (ρ , u , p ) and (ρ r , u r , p r ) with velocity c if the following conditions hold: where the function n is given by • for all nonnegative ϕ ∈ C 1 (R) with compact support, there holds: An admissible traveling wave solution is called a shock profile if furthermore the function 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 ±∞ are necessarily bounded. The restrictions inf R ± ρ > 0 and inf R ± p > 0 in Definition 2 ensure that the temperature θ = p/(R ρ) is well defined, belongs to E 2 , and is bounded from below by some positive constant. 1 In particular, we can then define the energy density n of photons as the unique solution in S (R) to the elliptic equation 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 θ make the entropy inequality in Definition 2 meaningful in the sense that all integrals involve only integrable functions. Eventually, if the traveling wave solution (ρ, u, p) is smooth, that is when ρ, 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.
Assume moreover that the following condition is satisfied: Then there exists a constant M > 0 such that for all states (ρ r , u r , p r , c) ∈ R 4 satisfying the following conditions: • the function is an entropic 1-shock wave for the Euler equations, there exists a shock profile solution to (1) that connects (ρ , u , p ) to (ρ r , u r , p r ) with velocity c, and ρ, u, p, θ are all discontinuous at 0. Moreover, the density profile ρ, the velocity profile u and the pressure profile p are monotone, while the temperature profile θ = p Rρ is increasing on R − and decreasing on R + . In particular, the maximum of θ on R is attained at 0 + , and there holds θ(0 − ) < θ r .
The definition of 1−shock wave is recalled in Sect. 3.2 below. We refer to Fig. 1  the density and temperature profiles. The peak of θ at 0 + is the so-called Zeldovich spike [19]. Note also that the maximal temperature of the precursor θ(0 − ) always remains smaller than the far field temperature θ r . According to the terminology introduced in [19], the profiles we are obtaining are thus always subcritical. We believe that the restrictions on γ 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 ODEs. Meanwhile, it is not clear that the result of Theorem 1 should hold for all values of γ .

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, could restore the continuity of the profile. The analysis of the system which incorporates those terms is beyond the scope of the present work and would involve deeper algebraic manipulations; instead we refer to the arguments developed in [16, p. 563 & pp. 579 ff.] and [19,Chapter VII,Section 18], and more recently to the discussion in [12].

The equations in dimensionless form
In order to simplify the calculations below, we first reduce (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 The new unknowns are Observing that the third equation in (1) reads we find that (1) equivalently reads We now introduce the characteristic length and time In the rescaled space and time variables The new constitutive relations for the fluid are In other words, we have reduced to the case R = σ = σ a σ s = 1 in (2). There is only one remaining dimensionless parameter in (7), which we define as The entropy inequality (3) is reduced to by introducing the rescaled entropy s := s/R that satisfies Eventually, we observe that the choice of the constants ρ 0 and p 0 is completely arbitrary. For instance, in view of Theorem 1, we can choose ρ 0 := ρ and p 0 := p . Let us observe that the corresponding temperature θ 0 coincides with θ but u 0 does not necessarily coincide with u (in particular because u 0 is positive while u has no reason to be positive). We adopt this convention from now on for ρ 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 τ, y by the dimensionless factor √ α, we have the following Definition.
Definition 3 (Admissible traveling waves for the dimensionless system) Let ( ρ , u , p ), ( ρ r , u r , p r ) ∈ R 3 with ρ , ρ r , p , p r > 0 be two given states and let c ∈ R. We shall say that ξ → ( ρ, u, p)(ξ ) ∈ (E 2 ) 3 defines an admissible traveling wave solution to (7) connecting ( ρ , u , p ) and ( ρ r , u r , p r ) with velocity c if the following conditions hold: • for all ϕ ∈ C 1 (R) with compact support, there holds: where the function n is given by • for all nonnegative ϕ ∈ C 1 (R) with compact support, there holds: An admissible traveling wave solution is called a shock profile if furthermore the function is an entropic shock wave for the Euler equations in dimensionless form.
Given a shock profile ( ρ, u, p ) ∈ E 3 2 in the sense of Definition 3, the corresponding weak solution to (7) is 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 ρ 0 and p 0 above, we are therefore reduced to constructing a shock profile for the dimensionless system (7), (8) with ρ = p = 1. In this case, we observe that the constant α 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 γ ∈ ]1, 5/3[, α > 0, and let ( ρ , u , p ) ∈ R 3 with ρ = p = 1. Assume moreover that the following condition is satisfied: Then there exists a constant M > 0 such that for all states is an entropic 1-shock wave for the Euler equations in dimensionless form, there exists a shock profile solution to (7) that connects ( ρ , u , p ) to ( ρ r , u r , p r ) with velocity c, and ρ, u, p are all discontinuous at 0. Moreover, the density profile ρ, the velocity profile u and the pressure profile p are monotone, while the temperature profile θ is increasing on R − and decreasing on R + . In particular, the maximum of θ on R is attained at 0 + , and moreover there holds θ(0 − ) < θ r .
In the scaling, we suppose that the material temperature scales like the radiation temperature. Then, in the rescaled system the dimensionless coefficient √ α measures the strength of the coupling between the material and radiation: it penalizes the coupling term n − θ 4 in the energy balance. Hence (11) assumes that the coupling is strong enough.

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 with the constitutive relations The system (12) is supplemented with the entropy inequality where s denotes the entropy that satisfies An entropic shock wave is a piecewise smooth solution 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 hold, and we shall say that (14) is a 3-shock if the inequalities 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.
Then the set of vectors (ρ r , u r , p r , c) ∈ R 4 such that (14) defines a 1-shock is a smooth curve that can be parameterized by p r ∈ ]1, +∞[. For each of these 1-shocks, we can define the quantities When the parameter p r tends to +∞, the following asymptotic expansions hold: Proof (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 Using the constitutive relation p = (γ − 1) ρ e, we end up with The asymptotic expansion of ρ r follows immediately. The expansion of j follows from the relation 2 The end of the proof of Lemma 1 follows from some elementary computations. 2 We recall that the mass flux j across the discontinuity is positive for a 1-shock and negative for a 3-shock.
The (positive) quantity p r − 1 gives a measure of the shock amplitude. When p r is close to 1, (ρ r , u r , p r ) is close to (ρ , u , p ) and the velocity c is close to the characteristic velocity u − √ γ p /ρ . For large shocks, that is when p r tends to +∞, all the quantities j, C, a behave like √ p r , and the velocity c tends to −∞. We also observe that the velocity u r varies monotonically along the shock curve, and therefore the quantity a determines the shock in a unique way. The asymptotic expansions of Lemma 1 will be used in Sect. 4 to solve the traveling wave equation that we shall derive in the following paragraph. From now on, we consider a fixed left state (ρ , u , p ) ∈ R 3 with ρ = p = 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 a is sufficiently large". Since we restrict to 1-shocks, all constants j, C, a are positive.

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 where the function f is defined as and μ is defined in Lemma 1. (It will be convenient from time to time to use the parameter μ instead of γ in order to simplify some expressions.) The function f depends on the shock wave (14) through the constants j and C. This is in sharp contrast with the situation in [7,18] where the system to be solved corresponds to (15) with f ≡ 1 and α = 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.
Proof (Proof of Proposition 1) Let us first assume that there exists a shock profile solution to (7) that connects (ρ , u , p ) to (ρ r , u r , p r ) with velocity c. Then if we define the shifted velocity U := u − c ∈ E 2 , we find that the derivative (in the sense of distributions) of both ρ U and ρ U 2 + p vanishes. This means that these piecewise smooth functions are constant, and examining the limit at ±∞, we obtain the relations Now we observe that θ ∈ 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 where sgn denotes the sign function. Furthermore, n is twice continuously differentiable on R \{0} and n has finite limits at 0 ± . Let us now define the radiative heat flux as Then q is continuous over R (as the convolution of a function in L 1 by a function in L ∞ ), and q tends to 0 at ±∞ (change variables η → ξ − η and apply the dominated convergence Theorem). Furthermore, q satisfies the following relations where, from now on, ∂ ξ always denotes the derivative in the sense of distributions: We can thus rewrite the third relation in Definition 3 as meaning that the function j E + p u + q is constant over R.
After a few simplifications, see [10] for similar computations, we obtain that the radiative heat flux q satisfies Let us now introduce the shifted velocity Then we derive the relations Since we know that q is continuous at 0, we derive the relation V (0 + ) = ±V (0 − ). Furthermore, if we differentiate the latter expression of q twice in the sense of distributions, we where stands for the classical derivative on R\{0} and δ denotes the Dirac mass at 0. It remains to equate the latter expression with 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 R − and R + (recall that we consider profiles with only one possible discontinuity of the first kind at 0), we obtain the inequality j (s(0 + ) − s(0 − )) ≥ 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 + ) ≥ 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 To complete the proof of Proposition 1, we now assume that (V, W ) ∈ E 2 × E 1 is a solution to (15) on R\{0} that satisfies conditions (ii), (iii), (iv). We define the functions which all belong to E 2 . The function ρ is well-defined and is bounded from below by a positive constant. Indeed, V is bounded from below by −a so U is bounded from below by u r − c > 0. In the same way, p is bounded from below by p > 0 because p is monotone. Let us also define the radiative heat flux q by the formula 2 , we obtain that q ∈ E 2 is continuous on R and tends to 0 at ±∞. Moreover, the definition of q and θ as well as the jump conditions at 0 yield the relation If we define the photons energy density n by the convolution formula of Definition 3, we obtain that q and −α ∂ ξ n are two temperate distributions that satisfy the same elliptic equation with constant coefficients. This gives q = −α ∂ ξ n, and therefore −∂ ξ q = n − θ 4 . The relations in Definition 3 that should be satisfied by an admissible traveling wave follow from the definition (19) of the functions ρ, u, p. (We feel free to skip the details.) 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 (a, 0) to (−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 R − , then on 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.

Construction of the heteroclinic orbit on R −
The first main task is to construct a solution to (15) on R − that tends to (a, 0) at −∞. 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: where the function f is still defined by (16). The transformation from (15) to (20) amounts to multiplying the vector field in the right-hand side of (15) by V . Our first main result for the modified system (20) is the following.
Then for a shock wave of large enough amplitude, the following results hold: (i) f (0) 2 > 2 a 2 /α and we thus define the quantity (ii) There exists a unique -up to a shift-solution (V 1 , W 1 ) to (20) that is defined on R and satisfies where the convergence is exponential. Moreover V 1 is a decreasing diffeomorphism, and − f (a)/2 < W 1 < 0.
The proof of Proposition 2 is based on a precise description of the phase portrait of (20) when the shock amplitude is large. The crucial point is to obtain a factorization of the second order polynomial in W that appears in (20). Let us observe right now that Point (i) in Proposition 2 is clear since, using the asymptotic expansions of Lemma 1, we have for large shocks. Here and from now on, κ denotes a generic positive constant. In particular, W 0 ∼ −κ p −2 r . Before proving point (ii) of Proposition 2, we establish a preliminary result. Lemma 2 Let γ ∈ ]1, 5/3[, and let α > 0 satisfy (11). If the shock amplitude is sufficiently large, then the polynomial Proof (Proof of Lemma 2) Using the asymptotic expansions of Lemma 1, we first compute the equivalent f (a) ∼ 4 μ > 0 for large shocks. In particular, f (a) is positive for large shocks. Therefore, we wish to show the inequality: Using Lemma 1, we find D(0) ∼ κ p 3 r 1 and D(a) = f (a)/2 ∼ 2 μ 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 , a], where V ∈ ]0, a[ depends on γ, α and on the shock. This fact will imply that the minimum of D on the interval [0, a] is attained for V = a and this minimum is a positive value (close to 2 μ). The behavior of D follows from a lengthy, though elementary, analysis.
We introduce the rescaled variable and the shifted variable We compute where we use the notation The expression of F simplifies with the introduction of the shifted variable Y , and we obtain Furthermore, the asymptotic expansions of Lemma 1 give for large shocks. In particular, Y max < Γ 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 (21) is increasing on an interval [0, X ] and is decreasing on the interval [X , X max ]. We compute 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 when the shock is sufficiently large (use Y max ∼ Γ and Y 1 < Y 2 < 1), and they are negative for Y = (γ − 1)/2 (use γ ∈ ]1, 5/3[ and 1/3 < Y 1 < Y 2 ). It follows that for sufficiently large shocks, 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 • For large shocks, the sign of D 1 (0) is the same as the sign of G ((γ − 1)/2) because the factor C 5 /j tends to +∞ while f (a) tends to 4 μ.
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 while we compute If (11) In particular, W 1 and W 2 satisfy where W 0 is defined in Proposition 2.
We refer to Fig. 2 for a representation of the graphs of W 1 and W 2 .
Proof (Proof of Lemma 3) Lemma 2 directly implies that for all V ∈ [0, a], the discriminant of the polynomial (in W ) The latter matrix has a unique positive eigenvalue λ + (the other eigenvalue is negative), that is given by In particular, λ + ∼ a/ √ α 1 for large shocks. An eigenvector associated with λ + is e + := (1, λ + /a), so we find e + ∼ (1, 1/ √ α) 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 (20) that is defined on an interval of the form ] − ∞, ξ 1 [, and such that when ξ tends to −∞, the orbit (V 1 , W 1 ) is tangent to the half-line (a, 0) + R − e + . In particular, for all ξ in a neighborhood of −∞, we have Let us show that the solution (V 1 , W 1 ) satisfies the properties of Proposition 2.
• We first prove that the orbit (V Of course, the orbit (V 1 , W 1 ) cannot pass through the stationary point (a, 0). Moreover, the vertical axis {0}×R is invariant by the flow of (20). Indeed, the restriction of the flow of (20) to the vertical axis coincides with the flow of the scalar ODE Since passes through a point of ∂ K , then this point belongs to the region K in which the vector field of (20) points toward the interior of K . We can conclude that (V 1 , W 1 ) does not pass through any point of ∂ 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 R 2 , the solution is globally defined. Moreover, we have V 1 < 0 so V 1 has a finite limit V < a as ξ tends to +∞. 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 +∞.
A detailed analysis is necessary to show that (V 1 , W 1 ) converges to (0, W 0 ) at +∞.
The ω-limit set of (V 1 , W 1 ) is a nonempty compact connected subset of the rectangle K . Since V 1 has a finite limit V at +∞, the ω-limit set of (V 1 , W 1 ) is of the form {V } × I where I is a closed subinterval of [− f (a)/2, 0]. Moreover, we know that the ω-limit set {V } × I is globally invariant by the flow of (20), and V ∈ [0, a[. Taking a closer look at the vector field of (20), we can conlude first that V equals 0. Indeed, if V is positive, then I cannot be restricted to one single point because, in such a case, the ω-limit set of (V 1 , W 1 ) would be a stationary point, but (20) has no stationary point with V > 0 as a first coordinate. Since we know that I contains at least two points, there exists W ∈ I with W < 0. However, the solution to (20) passing through (V , W ) does not satisfy V ≡ V so {V } × I is not invariant by the flow. We are led to a contradiction, hence V equals 0.
It remains to show that I is the singleton {W 0 }. On the vertical axis, the flow of (20) coincides with the flow of the ODE (23). Hence, I ⊂ [− f (a)/2, 0] is globally invariant by the flow of (23). We use Lemma 3, and observe that the two stationary points of (23) are (We recall that f (a) tends to 4 μ for large shocks.) Consequently, the only compact subinterval of [− f (a)/2, 0] that is globally invariant by the flow of (23) is {W 0 }. We have thus proved that the ωlimit set of (V 1 , W 1 ) is {(0, W 0 )}. In other words, (V 1 , W 1 )(ξ ) tends to (0, W 0 ) as ξ tends to +∞. A direct calculation shows that (0, W 0 ) is an attractive point (the eigenvalues of the Jacobian matrix are W 0 and −2 W 0 − f (0)) so the convergence is exponential.

Then for a shock wave of large enough amplitude, there exists a unique-up to a shift-solution
Moreover V is a decreasing diffeomorphism, and − f (a)/2 < W < 0.
Proof (Proof of Proposition 3) Using the exponential convergence of (V 1 , W 1 ) towards (0, W 0 ) at +∞, we can introduce the function Then Ξ 1 is an increasing C ∞ diffeomorphism from R to ] − ∞, 0[, and the solution (V , W ) to (15) is defined as

Construction of the orbit on R +
In this paragraph, we are going to construct a solution to (15) on R + that tends to (−a, 0) at +∞. Similarly to what we have done in the preceeding paragraph, we are first going to construct a solution to the auxiliary system (20) with suitable properties. Then a change of variables will yield the solution to (15). Our second main result for the auxiliary system (20) reads as follows.
Then for a shock wave of large enough amplitude, there exists a unique -up to a shift-solution (V 2 , W 2 ) to (20) that is defined on an interval of the form ] − ∞, ξ 2 [ (ξ 2 possibly equals +∞) and that satisfies for some V max ∈ ] − a, 0]. Moreover V 2 is an increasing diffeomorphism, W 2 is decreasing, and if we consider W 2 as a function of V ∈ ]−a, V max [, the inequality (29) below holds.
The proof of Proposition 4 is based on an accurate description of the vector field of (20) in the quarter plane {V < 0 , W < 0}. We begin with two preliminary results.
Proof (Proof of Lemma 4) The proof relies on simple though lengthy computations. Recalling the definition (22) of the function F, we compute Using γ ∈ ]1, 5/3[, we obtain F ≤ 0 on the interval [−(γ − 1)/2, 0] so F ≥ F (0) on this interval. When V belongs to ]−μ C/2, 0[, there holds X ∈ ]−(γ −1)/2, 0[, so we already have the lower bound Using the latter inequality, we can conclude that in the and it only remains to prove that the right hand side of (25) is positive for sufficiently large shocks. Equivalently, it is sufficient to prove for large shocks. We use the definition of W 0 and Lemma 1 to obtain and consequently for large shocks. It is now a simple exercise to check that this quantity is larger than 1 for γ ∈ ]1, 5/3[ (and μ ∈ ]0, 1/4[).
Lemma 5 Let γ ∈ ]1, 5/3[ and let α > 0. If the shock amplitude is sufficiently large, we have the following results: (i) There exists a unique pair of real numbers (V , V ) satisfying −a < V < −μ C/2 < V < 0 and Proof (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 := (γ + 1) V /C, and get where the polynomial function F, and the parameter X max were already defined in the proof of Lemma 2: Observe that for sufficiently large shocks, there holds Furthermore we know that the holomorphic function F 2 has exactly one double root in the closed disk Applying Rouché's Theorem, we find that the function H , defined in (27), 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 has no real root for V ∈ ]V , −μ C/2]. In particular, property (ii) holds for V in the interval ]V , −μ C/2], so we now assume V ∈ ] − a, V ]. In that case, (i) shows that (28) has two real roots, and it is not difficult to see that these roots are positive. Consequently, for any (V, W ) ∈ ] − a, V ] × R − , the inequality (26) holds. • The property (iii) is also a direct consequence of (i). For V ∈ [V , 0], Equation (28) has two real roots. These roots are simple for V = V and they are necessarily negative. We refer to Fig. 3 Fig. 3. The proof relies on Lemma 4 and splits in two steps.
Let us now prove Proposition 4.

Proof (Proof of Proposition 4)
We follow more or less the proof of Proposition 2. The stationary point (−a, 0) is a saddle point for the system (20). The Jacobian of the vector field at (−a, 0) equals the matrix , whose only positive eigenvalue (the other eigenvalue is negative) is The correspondent eigenvector is e + := (1, −λ + /a). Since f (−a) ∼ −κ p 3 r for large shocks, there holds λ + ∼ κ p 3 r , e + ∼ (1, −κ p 5/2 r ). By standard ODE theory, there exists a unique -up to a shift-solution (V 2 , W 2 ) to (20) that is defined on an interval of the form ] − ∞, ξ 2 [ and such that for ξ → −∞, the orbit tends to (−a, 0) with the tangent line (−a, 0)+ R + e + . In particular, for ξ in a neighborhood of −∞, there holds 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 ]−a, 0[ × ]−∞, 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 ξ<ξ 2 V 2 (ξ ). (This quantity is independent of the possible shift in the definition of (V 2 , W 2 ).) Since V 2 is a diffeomorphism from ] − ∞, ξ 2 [ to ] − a, V max [, we feel free to consider W 2 as a function of V ∈ ] − a, V max [, the orbit (V 2 , W 2 ) being the graph of W 2 . We obviously have W 2 (−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 (20), we get which gives We integrate this inequality from −a + ε to V ∈ ] − a, V max [ and then pass to the limit ε → 0. The function f is integrated by using the relation We thus obtain and, for further discussion, we denote the right hand side G(V ).
To complete the proof of Proposition 4, the only remaining thing to prove is lim V →V max W 2 (V ) = −∞. There are two possible cases according to the position of V max with respect to −μ C/2. so we feel free to consider W and W as functions of V . More precisely, W can be considered as a function of V ∈ [0, a] while W can be considered as a function of V ∈ [−a, V max [. Let us introduce the following two continuous functions: Lemma 1 shows that L(−a) is a large negative number for large shocks while K vanishes at −a, hence (K − L)(−a) > 0. Let us first consider the case V max ≤ −a + μ C (with the assumption 1 < γ < 5/3, Lemma 1 shows that −a + μ C is negative for large shocks). Then L remains bounded for V tending to V max while K tends to −∞. Consequently, the function K − L vanishes at some point, denoted Let us now consider the case V max > −a + μ C. We observe that with the above definition for L and with the function G defined in (29), there holds L(−a + μ C) = G(−a + μ C), so we deduce that is negative, where we have used (29). We can again conclude that the function K − L vanishes at some point, denoted V + , in the interval ] − a, −a + μ C[.
Up to using a shift in ξ , there is no loss of generality in assuming V (0) = −V + , V (0) = V + . Then the piecewise solution to (15) is defined by 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 + μ C[.

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 < γ < 5/3 made in Theorem 2, there holds 0 < μ < 1/4 and −a + μ C < −μ C/2 for large enough shocks (use Lemma 1).
In the preceeding 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 (19). In particular, u is decreasing (because V is decreasing) while ρ and p are increasing. The temperature profile is given by so we have In particular, θ (ξ ) has the same sign as 2 V (ξ ) + μ C. If ξ < 0, V (ξ ) is positive since it belongs to ]V (0 − ), a[, so θ (ξ ) is positive and θ is increasing on R − . If ξ > 0, we have shown in the preceeding paragraph that V (ξ ) belongs to the interval ] − a, −a + μ C[, so 2 V (ξ ) + μ C is negative and θ is decreasing on R + . To complete the proof of Theorem 2, we need to show θ(0 − ) < θ r . Using the expression (30), we have θ(0 − ) − θ r = (a + V (0 − )) (a − μ C − V (0 − )), and since V (0 − ) = −V (0 + ) > a − μ C, we get the desired inequality. The proof of Theorem 2 is thus complete.

Concluding remarks
In this paper, we prove the existence of discontinuous shock profiles for a PDEs system of radiative hydrodynamics. We deal with the simplest non equilibrium model corresponding to the low energy density approximation where radiations are described by a mere diffusion equation and the coupling reduces to energy exchanges due to emission and absorption; thus it is proportional to the difference of the emitted material energy and the radiative energy σ θ 4 − n. The gas is assumed non-relativistic and follows a polytropic law with adiabatic constant 1 < γ < 5/3; the asymptotic states are supposed at radiative equilibrium. It is remarkable that for the profiles we obtain, the temperature is non monotone and exceeds the far end temperature, a structure referred to as the Zeldovich spike. The occurrence of such structure assumes that the coupling between the material and the radiation is strong enough and the strength of the shock, namely the difference between the far end states, should be large enough.
Our analysis leaves open a couple of important questions: • It restricts to a quite simple model, which neglects several physical effects. Any extrapolation to a different situation should be considered with extreme caution. The role of additional coupling terms is far from clear; it is likely that modifying the model can restore the continuity of the profile for large shocks (as claimed in [19] under some physical assumptions). The arguments are highly dependent on the structure of the PDEs system and changing the system would certainly need a specific study. • The proof is partly based on asymptotic arguments and it does not provide explicit estimates on the strength of the shock that produces the temperature overshoot. It also heavily relies on the explicit (and simple) expression of the pressure law. By contrast to the situation investigated here, it is known that close enough far end states can be connected by smooth profiles [10]. Again the estimate of the shock strength is non explicit. Finding effective thresholds that discriminate between the smooth and singular behaviors can be an issue. • With the simple model investigated here, and using the scaling assumptions in Sect. 3.1, everything is encoded by the strength of the shock and the coupling parameter α. However, a description of intermediate behaviorslow coupling, moderate shock-is still lacking. As far as we have tried, this seems a difficult algebraic problem. However, since the completion of this work, a significant breaktrough in this direction has been achieved by Mascia [15].