Vincent Giovangigli

Zai-Bao Yang

Wen-An Yong

Relaxation Limit and Initial-Layers for a Class of Hyperbolic-Parabolic Systems

We consider a class of hyperbolic-parabolic systems with small diffusion terms and stiff sources. Existence of solutions to the Cauchy problem with ill prepared initial data is established by using composite expansions including initial-layer correctors and a convergence-stability lemma. New multitime expansions are introduced and lead to second-order error estimates between the composite expansions and the solution. Reduced equilibrium systems of second-order accuracy are also investigated as well as initial-layers of Chapman-Enskog expansions.

Introduction

Relaxation is an ubiquitous phenomenon of natural sciences typically modeled by systems of partial differential equations with stiff damping sources. This is a strong motivation for investigating the mathematical properties of these systems as well as the zero relaxation limit. Such systems of balance laws have notably been investigated mathematically in a hyperbolic framework as well as in a hyperbolic-parabolic framework [START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF][START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Natalini | Convergence to equilibrium for the relaxation approximation of conservation laws[END_REF][START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System with Stiff Source Terms[END_REF][START_REF] Yong | Existence of relaxation shock profiles for hyperbolic conservation laws[END_REF][START_REF] Dafermos | Hyperbolic Conservation Laws in Continuum Physics[END_REF][START_REF] Yong | Basic aspects of hyperbolic relaxation system, in advances in the theory of shock waves[END_REF][START_REF] Lattanzio | Hyperbolic-parabolic singular limits for first-order nonlinear systems[END_REF][START_REF] Kawashima | Dissipative structure and entropy for hyperbolic systems of conservation laws[END_REF][START_REF] Donatelli | Convergence and singular limits for multi-d semilinear hyperbolic systems to parabolic systems[END_REF][START_REF] Tzavaras | Relative entropy in hyperbolic relaxation[END_REF]. We study in this work existence of solutions and asymptotic expansions for hyperbolic-parabolic systems with small diffusion terms and stiff sources with ill prepared initial data.

We consider a nonlinear second-order system of partial differential equations in the form (1.1)

∂ t u + j∈D ∂ j F j (u) -ε i,j∈D ∂ i B ij (u)∂ j u = Q(u) ε ,
where u denotes the state vector of basic physical variables, t the time, ∂ i the derivative operator in the ith spatial direction, D = {1, . . . , d} the indexing set of spatial directions, d the spatial dimension, ε a small positive relaxation parameter, F j = F j (u), j ∈ D, the convective fluxes, B ij = B ij (u), i, j ∈ D, the diffusion matrices, and Q = Q(u) the source term. It is assumed that the convective fluxes, the source term, and the diffusion matrices are smooth functions of u ∈ O u where O u is an open convex set of R n . The parameter ε is typically the ratio of a characteristic relaxation time and a characteristic convective-diffusive time. The source term Q is assumed to be in quasilinear form and there exists a fixed linear space E of R n of dimension n e termed the slow manifold or the equilibrium manifold such that Q(u) ∈ E ⊥ . Denoting by Π e the matrix formed by column vectors of a basis of E, the projection u = Π t e u is then the natural conservative slow variable. Such systems of partial differential equations (1.1) notably arise from various applications like nonequilibrium twotemperature gases with fast relaxation of translational and internal temperatures [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF][START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation: Local convergence results, Nonlinear Analysis[END_REF], radiation hydrodynamics [START_REF] Yong | An interesting class of partial differential equations[END_REF], or reactive mixtures with fast chemistry [START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF][START_REF] Giovangigli | Solutions for Models of Chemically Reacting Compressible Mixtures[END_REF]. Existence results with uniform a priori estimates have been obtained for such system with stiff sources for well prepared initial data [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF][START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation: Local convergence results, Nonlinear Analysis[END_REF]. Denoting by u ε the solution to the Cauchy problem for the system (1.1) and u ε its initial value, well prepared initial data are such that Q( u ε ) = 0. The method of proof has been based on direct a priori estimates that couple standard hyperbolic-parabolic estimates with new contributions arising from stiff sources [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation: Local convergence results, Nonlinear Analysis[END_REF][START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF].

Our aim in this paper is to extend such previous existence results for hyperbolic-parabolic systems to the case of ill prepared initial data where Q( u ε ) = 0. In this situation, the source terms are initially O(1/ε) in (1.1) and this introduces a microscopic O(ε) scale in time, absent for well prepared inital data. This microscopic time scale naturally lead to inital layers involving the fast time τ = t/ε. The structure and exponential decrease O exp(-δτ ) for some δ > 0 of such initial layers are then key issues for the relaxation dynamics. Denoting by η a mathematical entropy for (1.1), the quasinormal variable is defined by w = (u, q) t where u is assumed to be in the form u = (u, v) t and q = (∂ v η) t . The transformation u → w is a smooth diffeomorphism, the quasinormal variable w is well suited to the source term Q and the slow manifold E, and we define w ε = w (u ε ). An Hilbert formal asymptotic expansion of w ε is built as ε goes to zero by considering outer as well as inner expansions with initial-layer correctors [START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System[END_REF][START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System with Stiff Source Terms[END_REF][START_REF] Yong | Basic aspects of hyperbolic relaxation system, in advances in the theory of shock waves[END_REF]. Using this formal expansion, an existence proof for the nonlinear system (1.1) is established as well as rigorous error estimates for ill prepared inital data. Key tools are a convergence-stability lemma for hyperbolic-parabolic systems which guarantee that the maximum existence time interval for system (1.1) remains bounded away from zero as ε → 0, the existence and properties of formal expansions, the exponential decrease of initial layers correctors, and the existence of a standard normal form. The accuracy of the composite expansion is shown to be O(ε 2 ) and is improved from a former O(ε 3/2 ) estimate [START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System[END_REF]. Such an error estimate is obtained by using new modified expansions including outer second-order fast-component correctors and does not require extra regularity assumptions.

The nonequilibrium system (1.1) further admits a second-order accurate reduced or equilibrium system governing the slow variable u in the form (1.2)

∂ t u e + j∈D ∂ j f e j (u e ) -ε i,j∈D ∂ i B e ij (u e )∂ j u e = 0,
where u e is the equilibrium solution, f e j = f e j (u e ) the equilibrium convective fluxes and B e ij = B e ij (u e ) the equilibrium diffusion matrices. These fluxes and diffusion matrices are defined for u e ∈ O ue where O ue is a convex open domain. Such a reduced system (1.2) may either be derived by using a twoterm Chapman-Enskog expansion [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF] or equivalently the Maxwellian iteration method that generally yields similar results at second-order [START_REF] Truesdell | Fundamental of Maxwell's Kinetic Theory of a Simple Monoatomic Gas[END_REF]. The resulting diffusion coefficients at equilibrium B e ij have contributions either directly inherited from the diffusion terms B ij out of equilibrium or arising from perturbed convective terms. This yields in particular the viscous tensor in a one-temperature fluid with the shear viscosity inherited from the non-equilibrium model and the volume viscosity coefficient associated with the relaxation of internal energy [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF]. Denoting by u ε (x, t) the nonequilibrium solution of (1.1), u ε = Π t e u ε the corresponding slow variable out of equilibrium, u e (x, t) the local equilibrium solution of (1.2), it has been established that for well prepared initial data [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF][START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation: Local convergence results, Nonlinear Analysis[END_REF], when the relaxation parameter ε is small, the error estimate is in the form u ε (x, t) -u e (x, t) = O(ε 2 ). In this work, for the more general situation of ill prepared initial data, we establish that the out of equilibrium solution is such that (1.3) u ε (x, t) -u e (x, t) + εu il (x, t/ε) = O(ε 2 ), where u il (x, t/ε) is an initial-layer corrector that decays exponentially to zero as t/ε goes to infinity. This estimate is also improved from a former O(ε 3/2 ) accuracy by using the estimates out of equilibrium obtained with the new formal expansions. The assumptions on system (1.1) as well as the local equilibrium approximation (1.2) are presented in Section 2. The formal asymptotic construction is presented in Section 3 and the mathematical justifications are presented in Section 4.

A class of hyperbolic-parabolic systems

The mathematical assumptions on the system of partial differential equations (1.1) are presented in this section and the second-order reduced system (1.2) is investigated.

Quasinormal variable

We denote by η a mathematical entropy defined on O u , an open convex set of R n . The entropy η = η(u), the fluxes F j = F j (u), the diffusion matrices B ij = B ij (u), and the source term Q = Q(u) are assumed to be C ∞ over the open convex set O u and the Jacobian matrices of convective fluxes are denoted by A i = ∂ u F i . The unit sphere in d dimension is denoted by Σ d-1 and x, y is the Euclidean product between vectors x and y. The system (1.1) is assumed to satisfy the following properties.

(u 1 ) The Hessian matrix ∂ 2 u η is positive definite over O u . (u 2 ) The products A j (∂ 2 u η) -1 , j ∈ D, are symmetric over O u .

(u 3 ) We have B ij (∂ 2 u η) -1 t = B ji (∂ 2 u η) -1
, i, j ∈ D, and for any ξ ∈ Σ d-1 the matrix B(u, ξ) = i,j∈D B ij (∂ 2 u η) -1 ξ i ξ j is positive semi-definite over O u .

(u 4 ) The source term reads Q = -L (∂ u η) t where L is a matrix of size n and there exists a partition R n = R ne ×R nr with n = n e + n r and L and Q in the form

(2.1) L = - 0 0 0 S , Q = - 0 0 0 S (∂ u η) t ,
where S is a symmetric positive definite matrix of size n r and u → S(u) is C ∞ over O u .

The assumptions (u 1 )-(u 2 ) associated with convective transport have been adapted from [START_REF] Godunov | An interesting class of quasilinear system[END_REF][START_REF] Friedrichs | System of conservation laws with a convex extension[END_REF], the assumption (u 3 ) associated with diffusion terms from [START_REF] Kawashima | System of a Hyperbolic-Parabolic Composite Type, with applications to the equations of Magnetohydrodynamics[END_REF][START_REF] Kawashima | Large-time Behavior of Solutions to Hyperbolic-Parabolic System of Conservations Laws and Applications[END_REF][START_REF] Kawashima | On the Normal Form of the Symmetric Hyperbolic-parabolic Systems Associated with the Conservation Laws[END_REF] and the assumption (u 4 ) concerning the source term from [START_REF] Yong | An interesting class of partial differential equations[END_REF]. A more general form for the source term Q introduced in [START_REF] Yong | An interesting class of partial differential equations[END_REF] may also be rewritten as in (u 4 ) by using a linear transform. The conditions (u 1 )(u 2 )(u 4 ) have recently been used as a criterion to construct constitutive equations in non-equilibrium thermodynamics [START_REF] Zhu | Conservation-dissipation formalism of irreversible thermodynamics[END_REF]. The C ∞ smoothness assumptions on the system coefficients may also easily be weakened [START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF].

An entropy balance equation may further be obtained by multiplying the natural system (1.1) by the classical entropic variable v = (∂ u η) t . The corresponding time derivative term is such that ∂ u η ∂ t u = ∂ t η whereas the convective terms read ∂ u ηA i ∂ i u = ∂ i β i where β i are the entropy fluxes with ∂ u η A i = ∂ u β i . The existence of such entropy fluxes is classically obtained from (u 2 ) and Poincaré Lemma, using that the image O v of O u under the map u → v is simply connected. The resulting balance equation for η is in the form (2.2)

∂ t η + i∈D ∂ i β i -ε i,j∈D ∂ i ∂ u ηB ij ∂ j u + ε i,j∈D ∂ 2 u ηB ij ∂ j u, ∂ i u - 1 ε ∂ u η Q = 0,
where the diffusion terms have been integrated by part and we have used that

∂ i v = (∂ 2 u η)∂ i u.
In particular, the physical entropy production rate due to the source term reads -∂ u η Q/ε and that due to gradients is related in Fourier space to the positive semi-definite diffusion matrix B(u, ξ)

= i,j∈D B ij (∂ 2 u η) -1 ξ i ξ j .
From Property (u 4 ), the equilibrium manifold or slow manifold is given by E = R ne ×{0} and n r = n -n e is the dimension of the fast manifold E ⊥ . The superscript or subscript e will generally be associated with the equilibrium or slow manifold and r associated with the fast or rapid manifold. The manifold E is naturally termed the equilibrium manifold since Q, x = 0 when x ∈ E so that the variable u, x is governed by a standard time dynamics and E ⊥ = {0} × R nr is naturally termed the fast manifold. We denote the block structure of convective fluxes and diffusion matrices as

(2.3) F j (u) = F e j F r j , B ij (u) = B e,e ij B e,r ij B r,e ij B r,r ij ,
and Π e = [e 1 , . . . , e ne ] is the natural projection operator over the equilibrium space E with e 1 . . . , e n denoting the basis vectors of R n . We also denote by u ∈ R ne and v ∈ R nr the slow and fast components of the basic physical variable u = (u, v) t and by q = (∂ v η) t the column vector associated with the partial derivative of entropy with respect to the variable v.

Using the block structure induced by the partition R n = R ne × R nr in (u 4 ), we may introduce the quasinormal variable (2.4)

w = u q = u (∂ v η) t .
This variable w is termed quasinormal since it shares similarities with standard normal variables. The variable w is naturally suited to the source term Q and thus especially convenient for investigating the asymptotic equilibrium limit ε → 0. It is easily established that the map u → w is a C ∞ diffeomorphism from the open set O u onto an open set O w of R n . It is indeed one to one thanks to the strict convexity (u 1 ) of η and the jacobian matrix

(2.5) ∂ u w = I ne 0 ∂ 2 vu η ∂ 2 v η
, where I k denotes the unit matrix of order k, is invertible since the matrix ∂ 2 v η is positive definite. We will denote for short by f j = f j (w ) and g j = g j (w ) the slow and fast convective fluxes in the jth direction as functions of the quasinormal variable w

f j (w ) = F e j u(w ) , g j (w ) = F r j u(w ) .
These fluxes f j and g j are C ∞ function of w since F e , and F r , are C ∞ as well as w → u. Similarly, the coefficients of the Hessian matrix ∂ 2 u η as functions of w are denoted for short by

(2.6) ∂ 2 u η u(w ) = η uu (w ) η uv (w )
η vu (w ) η vv (w ) .

We will sometimes commit the abuse of notation of denoting by the same letter functions of u or w but the context will always be clear. A balance equation for w is now established from (1.1) by using a change of variable.

Proposition 2.1. The quasinormal variable satisfies the system of partial differential equations

(2.7) w t + j∈D A j (w )∂ j w -ε i,j∈D ∂ i B ij (w )∂ j w = 1 ε Q (w ) + εd (w , ∂ x w ),
where the matrices

A j = ∂ u w A j ∂ w u, B ij = ∂ u w B ij ∂ w u, the source term Q = ∂ u w Q, and the quadratic residual d = -i,j∈D ∂ i (∂ u w ) B ij ∂ w u ∂ j w have the following properties. (w 1 ) The matrix A 0 = (∂ w u) t ∂ 2 u η ∂ w u is symmetric positive-definite block-diagonal. (w 2 ) The matrices A j A -1 0 , j ∈ D, are symmetric over O w . (w 3 ) We have B ij A -1 0 t = B ji A -1 0 , i, j ∈ D, and for any ξ ∈ Σ d-1 the diffusion matrix B = i,j∈D B ij A -1 0 ξ i ξ j is positive semi-definite over O w . (w 4 ) The source term Q = (0, Q r ) t is in the form Q = -L(w )w where Q r = -η vv S q = -L r,r q.
(w 5 ) The slow components of the quadratic residual d vanishes: d e = 0 over O w .

The matrices A 0 , A j , B ij , L, and the source term Q have regularity C ∞ and the quadratic residual may be written

(2.8) d = i,j∈D M ij (w )∂ i w ∂ j w ,
where the third order tensors M ij (w ) also have regularity C ∞ . Moreover, using the partitioning induced by the decomposition R n = R ne × R nr , the matrices A 0 , A j , B ij and L are given by

A 0 = η uu -η uv η -1 vv η vu 0 0 η -1 vv , (2.9) A j = ∂ u f j ∂ q f j η vu ∂ u f j + η vv ∂ u g j η vu ∂ q f j + η vv ∂ q g j , (2.10) B ij = B e,e ij -B e,r ij η -1 vv η vu η vu B e,e ij -B e,r ij η -1 vv η vu + η vv B r,e ij -B r,r ij η -1 vv η vu • • • (2.11) • • • B e,r ij η -1 vv η vu B e,r ij η -1 vv + η vv B r,r ij η -1 vv , L = A -1 0 0 0 0 S = 0 0 0 η vv S . (2.12)
Proof. Letting u = u(w ) in (1.1) and multiplying on the left by ∂ u w first yields the system (2.7). From A 0 = (∂ w u) t ∂ 2 u η ∂ w u and (u 1 ) we deduce that A 0 is symmetric positive definite and the blockdiagonal expression (2.9) is obtained by a direct calculation using (2.5) and (2.6). This yields (w 1 ) and the blocks η -1 vv and η uu -η uv η -1 vv η vu are both positive definite. We have

A j A -1 0 = ∂ u w A j (∂ 2 u η) -1 (∂ u w ) t
and this matrix is symmetric from (u 2 ) and (w 2 ) holds.

Similarly we have

B ij A -1 0 = ∂ u w B ij (∂ 2 u η) -1 (∂ u w ) t ,

and arguing as above yields that (B

ij A -1 0 ) t = B ji A -1 0 . In addition B(w , ξ) = i,j∈D B ij A -1 0 ξ i ξ j = ∂ u w i,j∈D B ij (∂ 2 u η) -1 ξ i ξ j (∂ u w ) t ,
so that the matrix B is positive semi-definite from (u 3 ) and (w 3 ) is established. The expressions (2.10) and (2.11) for A j and B ij are also obtained by a direct calculation. Finally, using (u 4 ) the modified source term Q = ∂ u w Q may be written Q = -L(w ) w with L(w ) given by (2.12) and (w 4 ) is established. The regularity class of A 0 , A j , B ij , L, and the source term Q is a direct consequence of (2.9)-(2.12), the regularity class of A 0 , A j , B ij , L, and Q, and the regularity class of ∂ u η, ∂ u w and ∂ w u.

Since d =i,j∈D ∂ i (∂ u w )B ij ∂ w u∂ j w , the first n e components vanish from (2.5). The kth component d k is also in the form

d k = 1≤l,l ′ ≤n (M ij ) kll ′ ∂ i w l ∂ j w l ′ with coefficients given by (M ij ) kll ′ = -1≤r,s≤n ∂ w l (∂ u r w k )(B ij ) rs ∂ w l ′ u s where the coefficients of B ij ∈ R n,n are denoted by (B ij ) rs .
The quadratic residual may thus be written as (2.8) where M ij (w ) are third order tensors that are C ∞ over O w and this completes the proof.

We introduce a structural assumption about the open set O w that will be convenient for investigating asymptotic expansions in the quasinormal variable w .

(u 5 ) The open set O w is in the form O w = O u × R nr where O u is a convex open set of R ne .
This assumption is natural since the fast variable q = (∂ v η) t is the Legendre conjugate of v and its components are slopes of the convex function η. The assumption (u 5 ) thus means that the gradient of η is infinite at the boundaries of the domain O u which is natural from a thermodynamic point of view. This condition also corresponds to the notion of 'essentially smooth functions' in convex analysis [33]. In addition, the q component will converge rapidly to zero during the initiallayers of multitime expansions in such a way that it is a natural requirement that such initial-layer trajectories in the form (u, q(τ )) t -where u ∈ O u is fixed-lie in O w . Such a structure condition holds in particular in the situation of two temperature fluids [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF] as well as complex chemistry fluids [START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF]. With (u 5 ) the open set O w is also convex as O u and this will allow to derive various differential identities.

Equilibrium limit and reduced system

We study equilibrium states and then investigate formally the second-order reduced system as ε → 0. Proposition 2.2. For any u ∈ O u , letting w = w (u) = (u, q) t , the following properties are equivalent:

(i) Physical entropy production due to the source term vanishes: -∂ u η Q = 0.

(ii) The source term and the fast component vanish: Q = 0 and q = 0.

(iii) The entropic variable v and the quasinormal variable w belong to the equilibrium manifold E.

Any u ∈ O u that satisfies these equivalent properties is termed an equilibrium state as well as its image w under the map u → w .

Proof. These properties are directly obtained from the structure of the source term and using that E = R ne × {0}. Physical entropy production due to source terms is in the form

-∂ u η Q = -(∂ u η) t , Q = Sq, q so that ∂ u η Q = 0 is equivalent to q = 0 or Q = -(0, Sq) t = 0 since S is symmetric positive definite. Moreover, since v = (∂ u η) t ,
q) t , we have v ∈ E if and only if q = 0, and similarly w ∈ E if and only if q = 0 and this completes the proof. Such properties of the source term have been discussed in an abstract framework by Chen, Levermore and Liu [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF]. These properties hold in particular for gas mixtures with complex chemistry [START_REF] Giovangigli | Multicomponent flow modeling[END_REF][START_REF] Giovangigli | Mathematical modeling of supercritical multicomponent reactive fluids[END_REF][START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF] as well as two temperature fluids [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF]. Another interesting property of the source term is that the jacobian matrix at an equilibrium states is symmetric [START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System[END_REF][START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF]. The set of equilibrium points characterized by the condition q = 0 is naturally parametrized by the slow variable u. Indeed, for any u ∈ O u , there exists a unique equilibrium point w eq ∈ O w such that Π t e w eq = u simply given by w eq = (u, 0) t . Similarly, there exists a unique equilibrium point u eq ∈ O u such that Π t e u eq = u simply given by u eq = u(w eq ). Letting u eq = (u, v eq ) t , the map u → v eq (u) is then C ∞ over O u and by differentiating the equilibrium relation q u, v eq (u) = 0 it is obtained that ∂ u v eq = -(η vv ) -1 η vu where we have denoted for short η vu (u) = η vu (w eq ) = η vu (u, 0) and η vv (u) = η vv (w eq ) = η vv (u, 0). Lemma 2.3. Defining the entropy at equilibrium by η e (u) = η u eq (u) for u ∈ O u then u → η e (u) is C ∞ over O u , ∂ u η e = ∂ u η(u eq ) and

(2.13) ∂ 2 u η e = η uu -η uv (η vv ) -1 η vu = A e,e 0 (u, 0), where η uu (u) = η uu w eq (u) = η uu (u, 0) and η uv (u) = η uv w eq (u) = η uv (u, 0).

Proof. Differentiating the entropy η e first yields that ∂ u η e = ∂ u η u, v eq (u) since ∂ v η u, v eq (u) = 0. Differentiating again and using ∂ u v eq = -(η vv ) -1 η vu then yields (2.13) and ∂ 2 u η e = η uuη uv (η vv ) -1 η vu coincides with A e,e 0 (u, 0) from (2.9). Our aim is now to obtain a second-order accurate reduced system for the slow variable u as ε → 0 for the standard time dynamics associated with time t. The fast component of Equation (2.7) may first be rewritten as

q = -εS -1 η -1 vv ∂ t q + η vu j∈D ∂ j f j (u, q) + η vv j ∂ j g j (u, q) + O(ε) .
We thus obtain formally that q = O(ε) and thus that ∂ t q = O(ε) for the standard time dynamics, so that near equilibrium

q = -εS -1 η -1 vv j∈D η vu ∂ u f j (u, 0) + η vv ∂ u g j (u, 0) ∂ j u + O(ε 2 ),
where S(u) = S u eq (u) . On the other hand, the slow component of Equation (2.7) may be rewritten

(2.14) ∂ t u + j∈D ∂ j f j -ε i,j∈D ∂ i B e,e ij -B e,r ij η -1 vv η vu ∂ j u + B e,r ij η -1
vv ∂ j q = 0, using (2.10) and (2.11). We may further expand the coefficients near equilibrium B e,e ij (u) -B e,r ij (u)η -1 vv η vu = B e,e ij (u eq ) -B e,r ij (u eq )η -1 vv η vu + O(ε), since ww eq = (0, q) t = O(ε) and uu eq = O(ε) as well as expand the flux f j (u, q) in (2.14) into f j (u, q) = f j (u, 0) + ∂ q f j (u, 0)q + O(ε 2 ). Combining these expansions with (2.14) yields the following reduced system called the second-order equilibrium approximation of (1.1)

(2.15)

∂ t u + j∈D A e j (u)∂ j u -ε i,j∈D ∂ i B e ij (u)∂ j u = 0,
where f e j (u) = f j (u, 0) = F e u eq (u) , A e j (u) = ∂ u f e j (u) = A e,e j (u, 0), and

B e ij (u) = B e,e ij (u eq ) -B e,r ij (u eq )η -1 vv η vu + ∂ q f i (u, 0)S -1 η -1 vv η vu ∂ u f j (u, 0) + η vv ∂ u g j (u, 0) .
The matrices A e j , j ∈ D, and B e ij , i, j ∈ D, defined over O u have regularity C ∞ . Further using the expression (2.11) of B ij as well as the relation

∂ q f i (u, 0) t ∂ 2 u η e = η -1 vv η vu ∂ u f j (u, 0) + η vv ∂ u g j (u, 0
) deduced from the symmetry of A 0 A j , the equilibrium diffusion matrix B e ij can be written (2.16)

B e ij (u) = B e,e ij (u eq ) + ∂ q f i (u, 0) S -1 (∂ q f j (u, 0)) t ∂ 2 u η e .
This procedure to derive (2.15) corresponds to a second-order Chapman-Enskog expansion using the variable w [START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF] or equivalently a Maxwellian iteration that generally yields similar results at second-order [START_REF] Truesdell | Fundamental of Maxwell's Kinetic Theory of a Simple Monoatomic Gas[END_REF]. The first term B e,e ij (u eq ) in the expression of the equilibrium diffusion matrices (2.16) is inherited from the original diffusion term B ij whereas the second term

∂ q f j (u, 0) S -1 (∂ q f i (u, 0)) t ∂ 2
u η e comes from relaxation of the fast variable q. It leads in particular to the complete viscous tensor in one-temperature equilibrium fluid models with the shear viscosity inherited from the non-equilibrium two-temperature system and the volume viscosity coefficient associated with the relaxation of internal energy [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF][START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation: Local convergence results, Nonlinear Analysis[END_REF]. The structure of the reduced system at equilibrium is summarized in the following proposition. (u 2 ) The products A e j (∂ 2 u η e ) -1 , j ∈ D, are symmetric over u ∈ O u .

(u 3 ) We have B e ij (∂ 2 u η e ) -1 t = B e ji (∂ 2 u η e ) -1
, i, j ∈ D, and for any ξ ∈ Σ d-1 the diffusion matrix

B e = i,j∈D B e ij (∂ 2 u η e ) -1 ξ i ξ j is positive semi-definite over O u .
Proof. It has been established in Lemma 2.3 that ∂ 2 u η e = η uu -η uv (η vv ) -1 η vu so that it is positive definite as a Schur complement of the positive definite matrix η vv in the positive definite matrix ∂ 2 u η evaluated at u eq (u) and (u 1 ) is established.

From the symmetry of A j A -1 0 established in (w 2 ) and evaluated at the equilibrium point u eq , we deduce that the products A e j (∂

B e = i,j∈D B e,e ij (∂ 2 u η e ) -1 ξ i ξ j + i∈D ∂ q f i (u, 0)ξ i S -1 j∈D ∂ q f j (u, 0)ξ j t ,
and the positive semi-definiteness of B e is a consequence from that of B and this establishes (u 3 ).

The reduced equilibrium system (2.15) has been obtained by using a Chapman-Enskog expansion with the quasinormal variable w . It is established in Appendix A. that it naturally coincides with that obtained in previous work [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF]. In particular, the properties of the entropic symmetrized equilibrium system obtained in [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF] also applies to (2.15). We finally investigate the relative entropy with respect to equilibrium that is a natural distance to equilibrium. Relative entropies have been found to be key tools for investigating hyperbolic as well as hyperbolic-parabolic systems [START_REF] Dafermos | The second law of thermodynamics and stability[END_REF][START_REF] Diperna | Uniquenass of solutions to hyperbolic conservation laws[END_REF][START_REF] Kawashima | System of a Hyperbolic-Parabolic Composite Type, with applications to the equations of Magnetohydrodynamics[END_REF][START_REF] Dafermos | Hyperbolic Conservation Laws in Continuum Physics[END_REF][START_REF] Yong | Entropy and global existence for hyperbolic balance laws[END_REF][START_REF] Tzavaras | Relative entropy in hyperbolic relaxation[END_REF][START_REF] Feireisl | Relative entropies, suitable weak solutions, and weakstrong uniqueness for the compressible Navier?Stokes system[END_REF][START_REF] Lattanzio | Relative entropy in diffusive relaxation[END_REF][START_REF] Donatelli | Low Mach quantum number limit for the quantum-hydrodynamics system[END_REF][START_REF] Bresch | Relative entropy for compressible Navier-Stokes equations with density dependent viscosities and applications[END_REF][START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF][START_REF] Christoforou | Relative entropy for hyperbolic-parabolic systems and application to the constitutive equations of thermoelasticity[END_REF]. The relative entropy with respect to equilibrium is nonnegative by the convexity of η and locally behaves quadratically with respect to q.

Lemma 2.5. The relative entropy with respect to equilibrium (2.17) η il (u, q) = η u(u, q) -η u(u, 0) -∂ u η u(u, 0) u(u, q)u(u, 0) , is such that ∂ q η il = q t η -1 vv and may be written η il (u, q) = Hq, q where H = 1 0 η -1 vv u, αq)αdα is positive definite matrix.

Proof. Using u(u, q)u(u, 0) = 0, v(u, q) -v(u, 0) as well as that ∂ u η u(u, 0) ∈ E from Proposition 2.2 it is first obtained that ∂ u η u(u, 0) u(u, q)u(u, 0) = 0 so that

η il (u, q) = η u(u, q) -η u(u, 0) = η u, v(u, q) -η u, v(u, 0) .
Differentiating with respect to q we obtain ∂ q η il = ∂ v η∂ q v = q t η -1 vv and by calculus η il (u, q) = 1 0 (∂ v η∂ q v)(u, αq) dα q keeping in mind that all states (u, αq) t are in O w from (u 5 ). The identity η il (u, q) = Hq, q then follows since ∂ v η u, v(u, αq) = αq t and ∂ q v u, v(u, αq) = η -1 vv u, αq). The matrix H is finally positive definite as an average of positive definite matrices over a compact set.

Formal asymptotic expansion

The solution to the Cauchy problem for system (2.7) with initial data w ε is denoted by w ε = (u ε , q ε ) t . A formal approximation of w ε is built as ε goes to zero by using composite or multitime expansions [START_REF] O'malley | Introduction to Singular Perturbations[END_REF][START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System[END_REF]. The variable w is more convenient than u since it is adapted to the slow manifold E. Such a construction is similar to that for hyperbolic systems [START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System with Stiff Source Terms[END_REF] but now includes diffusion terms.

Composite expansion

We seek a composite expansion in the form

(3.1) w ε (x, t) = w 0 (x, t) + εw 1 (x, t) + w il 0 (x, t/ε) + εw il 1 (x, t/ε) + O(ε 2 ),
where w 0 and w 1 depend on the standard time t whereas w il 0 and w il 1 depend on the fast time τ = t/ε. The expansion (3.1) thus includes an outer expansion

(3.2) w 0 (x, t) + εw 1 (x, t) + O(ε 2 ),
associated with standard time dynamics and standard time t. This outer expansion (3.2), however, cannot generally satisfy the prescribed initial value w ε . It is then necessary to add a fast evolving initial-layer corrector in the form

w il 0 (x, t/ε) + εw il 1 (x, t/ε) + O(ε 2 ),
associated with fast time dynamics and fast time τ = t/ε. The initial-layer correctors w il 0 and w il 1 have to decay exponentially to zero as τ goes to infinity from general matching principles [START_REF] Eckhaus | Matching principles and composite expansions[END_REF]. The complete approximation of w ε is then a multitime expansion since it involves both the standard time t and the fast time τ = t/ε. We will notably use the second-order truncated approximation w a ε defined by

(3.3) w a ε (x, t) = w 0 (x, t) + εw 1 (x, t) + w il 0 (x, t/ε) + εw il 1 (x, t/ε).
This four term expansion w a ε will be shown to be well defined over some time interval independent of ε and used to establish an existence theorem for the out of equilibrium system (2.7). The expansion w a ε will also be shown to be an accurate approximation of the out of equilibrium solution w ε in the next section.

As a formal solution of (2.7) with initial value w ε , the composite expansion must also take the prescribed initial value

w 0 (x, 0) + εw 1 (x, 0) + w il 0 (x, 0) + εw il 1 (x, 0) + O(ε 2 ) = w ε (x).
Assuming naturally that w ε has an expansion

w ε = w 0 + ε w 1 + O(ε 2 ), we get (3.4) w 0 (x, 0) + w il 0 (x, 0) = w 0 (x), w 1 (x, 0) + w il 1 (x, 0) = w 1 (x).
The difficulty is then to find the proper initial values for w 0 and w 1 , based on the above initial value relations, as well as to determine the initial-layer correctors w il 0 and w il 1 by using the matched asymptotic expansion principle [START_REF] Eckhaus | Matching principles and composite expansions[END_REF][START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System[END_REF].

For the projections u ε and q ε the composite expansions are denoted by

u ε = u 0 (x, t) + εu 1 (x, t) + u il 0 (x, t/ε) + εu il 1 (x, t/ε) + O(ε 2 ), q ε = q 0 (x, t) + εq 1 (x, t) + q il 0 (x, t/ε) + εq il 1 (x, t/ε) + O(ε 2 ),
using the straightforward notation w ε = (u ε , q ε ) t , w 0 = (u 0 , q 0 ) t , w 1 = (u 1 , q 1 ) t , w il 0 = (u il 0 , q il 0 ) t , and w il 1 = (u il 1 , q il 1 ) t . From (3.4) we also have the relations

u 0 (x, 0) + u il 0 (x, 0) = u 0 (x), u 1 (x, 0) + u il 1 (x, 0) = u 1 (x), (3.5) q 0 (x, 0) + q il 0 (x, 0) = q 0 (x), q 1 (x, 0) + q il 1 (x, 0) = q 1 (x), (3.6)
where w 0 = ( u 0 , q 0 ) t and w 1 = ( u 1 , q 1 ) t .

Outer expansion

The equations governing the coefficients w 0 and w 1 are obtained by substituting the outer expansion (3.2) into system (2.7) and equating each power of ε. The governing equations obtained for w 0 and w 1 are then projected onto the slow and fast manifolds E and E ⊥ in order to get the governing equations for u 0 , q 0 , u 1 and q 1 . Equating the kth power of ε generally yields an equation for u k and q k+1 so that there is no clear separation of the asymptotic orders.

At the order ε -1 we first obtain from (2.7) that Q (w 0 ) = 0 and from the structure of Q = (0, Q r ) t with Q r = -η vv Sq this yields that (3.7) q 0 = 0.

The zeroth order solution w 0 = (u 0 , q 0 ) t = (u 0 , 0) t thus represents equilibrium states. At the next order ε 0 it is obtained that

(3.8) ∂ t w 0 + j∈D A j (w 0 )∂ j w 0 = -L(w 0 )w 1 .
Projecting on the slow manifold, using q 0 = 0 and the structure (2.10) of A j , we obtain (3.9)

∂ t u 0 + j∈D ∂ j f j (u 0 , 0) = 0.
From the structural properties of the out of equilibrium system and its symmetrizability, the system (3.9) is a hyperbolic system of conservation laws in such a way that there exists a unique local solution starting from an appropriate initial value u 0 (0). On the other hand, the fast component of (3.8) yields

(3.10) j∈D η vu (u 0 , 0)∂ u f j (u 0 , 0) + η vv (u 0 , 0)∂ u g j (u 0 , 0) ∂ j u 0 = -L r,r (w 0 )q 1 ,
where L r,r (w 0 ) = η vv (u 0 )S(u 0 ). This relation uniquely defines q 1 since L r,r is invertible keeping in mind that both η vv u 0 ) and S(u 0 ) are invertible. Finally, at the order ε 1 , the u 1 equation on the slow manifold is found in the form

(3.11) ∂ t u 1 + j∈D ∂ j ∂ u f j (u 0 , 0)u 1 + j∈D ∂ j ∂ q f j (u 0 , 0)q 1 - i,j∈D ∂ i B e i,j (u 0 , 0)∂ j u 0 = 0.
This is again a hyperbolic system that uniquely defines u 1 for an appropriate initial condition u 1 (0). In summary, from the symmetry properties of the equilibrium fluxes f e j established in Proposition 2.4, u 0 satisfies a symmetrizable hyperbolic system (3.9) and standard existence theory [START_REF] Godunov | An interesting class of quasilinear system[END_REF][START_REF] Friedrichs | System of conservation laws with a convex extension[END_REF][START_REF] Majda | Compressible Fluid Flow and System of Conservation Laws in Several Space Variables[END_REF][START_REF] Kawashima | System of a Hyperbolic-Parabolic Composite Type, with applications to the equations of Magnetohydrodynamics[END_REF][START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF][START_REF] Dafermos | Hyperbolic Conservation Laws in Continuum Physics[END_REF][START_REF] Kawashima | Dissipative structure and entropy for hyperbolic systems of conservation laws[END_REF] yields a local classical solution u 0 (x, t) once u 0 (x, 0) is given appropriately. The perturbed fast component q 1 (x, t) is then explicitely given in terms of u 0 and its first spatial derivatives ∂ x u 0 . Once u 0 (x, t) and q 1 (x, t) are determined, the system (3.11) is linear symmetrizable hyperbolic for u 1 and standard existence theory yields a local classical u 1 (x, t) once the initial value u 1 (x, 0) is given appropriately. A key issue is therefore to find the initial values u 0 (•, 0) and u 1 (•, 0) and this requires to further study the composite expansions (3.1).

Inner expansion

In order to derive the equations governing the initial-layer correctors w il 0 and w il 1 associated with the fast time τ = t/ε and to investigate composite expansions, it is first required to write the outer expansion in terms of τ by using t = ετ . To this aim we note that

w 0 (x, ετ ) = w 0 (x, 0) + ετ ∂ t w 0 (x, 0) + ε 2 τ 2 1 0 1 0 ∂ 2 t w 0 (x, εαβτ )αdαdβ, and letting R ′ 0 = τ 2 1 0 1 0 ∂ 2 t w 0 (x, εαβτ )αdαdβ and R 0 = τ ∂ t w 0 (x, 0) + εR ′ 0 yields w 0 = w 0 (x, 0) + εR 0 and w 0 = w 0 (x, 0) + ετ ∂ t w 0 (x, 0) + ε 2 R ′ 0 .
Similarly, we have

w 1 (x, ετ ) = w 1 (x, 0) + ετ 1 0 ∂ t w 1 (x, εατ )dα,
and letting

R 1 = τ 1 0 ∂ t w 1 (x, εατ )dα yields w 1 = w 1 (x, 0) + εR 1 .
By combining these expressions we obtain (3.12)

w 0 (x, t) + εw 1 (x, t) + O(ε 2 ) = w p 0 (x, τ ) + εw p 1 (x, τ ) + O(ε 2 ),
where

(3.13) w p 0 (x, τ ) = w 0 (x, 0), w p 1 (x, τ ) = w 1 (x, 0) + τ ∂ t w 0 (x, 0).
The superscript stands for 'polynomial' since w p i is generally found to be a polynomial in τ of the ith degree [START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System with Stiff Source Terms[END_REF]. The relations (3.13) may also be written component wise u p 0 = u 0 (x, 0), q p 0 = q 0 (x, 0) = 0, u p 1 = u 1 (x, 0) + τ ∂ t u 0 (x, 0), and q p 1 = q 1 (x, 0), where w p 0 = (u p 0 , q p 0 ) t and w p 1 = (u p 1 , q p 1 ) t . The resulting inner expansion in terms of τ is then

(3.14) w ε (x, τ ) = w p 0 (x, τ ) + w il 0 (x, τ ) + ε w p 1 (x, τ ) + w il 1 (x, τ ) + O(ε 2 ),
and should asymptotically satisfy system (2.7). Substituting the inner expansion (3.14) into (2.7) written in terms of the fast time τ , expanding into power series of ε, and equating the coefficients yield the governing equations for w p 0 + w il 0 and w p 1 + w il 1 . At the order ε -1 we obtain that

(3.15) ∂ τ (w p 0 + w il 0 ) = -L(w p 0 + w il 0 ) (w p 0 + w il 0 ) = Q (w p 0 + w il 0 ),
and projecting on the slow and fast manifolds yields

∂ τ (u p 0 + u il 0 ) = 0, (3.16) ∂ τ (q p 0 + q il 0 ) = Q r (w p 0 + w il 0 ). (3.17) At the order ε 0 it is next obtained (3.18) ∂ τ (w p 1 + w il 1 ) = - j∈D ∂ j A j (w p 0 + w il 0 )∂ j (w p 0 + w il 0 ) + ∂ w Q (w p 0 + w il 0 ) (w p 1 + w il 1 ).
Projecting for the slow and fast components yields that

∂ τ (u p 1 + u il 1 ) = - j∈D ∂ j f j (w p 0 + w il 0 ), (3.19) ∂ τ (q p 1 + q il 1 ) = ∂ u Q r (w p 0 + w il 0 ) (u p 1 + u il 1 ) (3.20) + ∂ q Q r (w p 0 + w il 0 ) (q p 1 + q il 1 ) + C(w p 0 + w il 0 ),
where

C(w p 0 + w il 0 ) = - j∈D η vu (w p 0 + w il 0 )∂ j f j (w p 0 + w il 0 ) + η vv (w p 0 + w il 0 )∂ j g j (w p 0 + w il 0 ) .
On the other hand, the outer expansion w 0 (x, t) + εw 1 (x, t) + O(ε 2 ) may also be rewritten in terms of the fast time w p 0 (x, τ ) + εw p 1 (x, τ ) + O(ε 2 ) and this expression naturally satisfies the system (2.7). This procedure yields the governing equations satisfied by u p 0 and u p 1 . Equivalently, we may let ε → 0 in the governing equations for u p 0 + u il 0 and u p 1 + u il 1 , keeping in mind that u il 0 and u il 1 converge exponentially to zero as τ → ∞. At order ε -1 this yields after some algebra

∂ τ w p 0 = -L(w p 0 ) w p 0 = Q (w p 0 ) = 0, (3.21) ∂ τ u p 0 = 0, ∂ τ q p 0 = Q r (w p 0 ). (3.22)
Similarly, at the order ε 0 , we obtain that

∂ τ w p 1 = - j∈D ∂ j A j (w p 0 )∂ j w p 0 + ∂ w Q (w p 0 ) w p 1 , (3.23)
∂ τ u p 1 = - j∈D ∂ j f j (w p 0 ), ∂ τ q p 1 = ∂ u Q r (w p 0 )u p 1 + ∂ q Q r (w p 0 )q p 1 + C(w p 0 ), (3.24)
where C(w p 0 ) = -j∈D η vu (w p 0 )∂ j f j (w p 0 ) + η vv (w p 0 )∂ j g j (w p 0 ) . The relations (3.21)-(3.24) may also be derived directly from the expressions of the outer expansion coefficients w p 0 and w p 1 and the governing equations for u 0 , q 0 , u 1 and q 1 .

Initial-layer correctors

We investigate the initial-layer correctors w il 0 and w il 1 by using the governing equations (3.15)-(3.17) and (3.18)-(3.20). Both w p 0 and w p 1 have been expressed in terms of the outer expansion and in particular, u p 0 (x, τ ) = u 0 (x, 0), q p 0 (x, τ ) = q 0 (x, 0) = 0, u p 1 (x, τ ) = u 1 (x, 0) + τ ∂ t u 0 (x, 0), and q p 1 (x, τ ) = q 1 (x, 0) keeping in mind that τ ∂ t q 0 (x, 0) = 0 since q 0 (x, t) = 0. At order ε -1 we obtain from (3.15) and (3.21) that

(3.25) ∂ τ w il 0 = Q (w p 0 + w il 0 ) -Q (w p 0 ) = Q (w p 0 + w il 0 ),
and projecting on the slow and fast manifold-or equivalently using (3.16), (3.17), and (3.22)-we obtain that ∂ τ u il 0 = 0 and ∂ τ q il 0 = Q r (u p 0 + u il 0 , q il 0 ) keeping in mind that q p 0 (x, τ ) = 0. Using ∂ τ u il 0 = 0 and that the initial-layer corrector w il 0 (x, τ ) goes to zero as τ goes to infinity, we deduce that

(3.26) u il 0 (x, τ ) = u il 0 (x, ∞) = 0.
Combining u il 0 (x, τ ) = 0, u p 0 (x, τ ) = u 0 (x, 0), q 0 (x, 0) = 0, with w 0 (x, 0) + w il 0 (x, 0) = w 0 (x) obtained from (3.4) we get (3.27) u 0 (x, 0) = u 0 (x), q il 0 (x, 0) = q 0 (x).

In particular w 0 (x, 0) = w p 0 (x, 0) = u 0 (x), 0 t , w il 0 (x, 0) = 0, q 0 (x) t , and for any τ ≥ 0

w p 0 (x, τ ) = u 0 (x), 0 t , w il 0 (x, τ ) = 0, q il 0 (x, τ ) t .
In particular, q il 0 (x, τ ) satisfies the following nonlinear ordinary differential equation involving solely the fast time τ

(3.28) ∂ τ q il 0 = Q r u 0 (x), q il 0 , q il 0 (x, 0) = q 0 (x).
The existence of a global in time classical solution q il 0 (x, τ ) that decays exponentially to zero as τ → ∞ will be established by using Q r ( u 0 , q il 0 ) = -η vv ( u 0 , q il 0 ) S u( u 0 , q il 0 ) q il 0 , the mathematical entropy η, as well as the symmetry and positive definiteness of η vv and S.

Considering next the first order initial-layer corrector w il 1 (x, τ ) we get from (3.18) and (3.23) that

∂ τ w il 1 = - j∈D A j (w p 0 + w il 0 )∂ j (w p 0 + w il 0 ) -A j (w p 0 )∂ j w p 0 (3.29) + ∂ w Q (w p 0 + w il 0 ) (w p 1 + w il 1 ) -∂ w Q (w p 0 ) w p 1 .
Projecting on the slow and fast manifolds-or equivalently using (3.19), (3.20), and (3.24)-we obtain the following differential equations involving solely the fast time τ

∂ τ u il 1 = - j∈D ∂ j f j (w p 0 + w il 0 ) -∂ j f j (w p 0 ) , (3.30) ∂ τ q il 1 = ∂ u Q r (w p 0 + w il 0 )(u p 1 + u il 1 ) -∂ u Q r (w p 0 )u p 1 (3.31) + ∂ q Q r (w p 0 + w il 0 )(q p 1 + q il 1 ) -∂ q Q r (w p 0 )q p 1 + C(w p 0 + w il 0 ) -C(w p 0 ).
The matching principle yields that

(3.32) u il 1 (x, ∞) = 0, q il 1 (x, 0) = q 1 (x) -q 1 (x, 0),
where q 1 (x, 0) can be obtained from (3.10)

q 1 (x, 0) = -S( u 0 ) -1 j∈D η -1 vv ( u 0 )η vu ( u 0 )∂ u f j ( u 0 , 0) + ∂ j g j ( u 0 , 0) ∂ j u 0 .
In particular, we may evaluate u il 1 (x, τ ) from (3.30) and (3.32) in the form

(3.33) u il 1 (x, τ ) = ∞ τ j∈D ∂ j f j w p 0 + w il 0 (τ ′ ) -∂ j f j (w p 0 ) dτ ′ ,
and u 1 (x, 0) is then finally given by u 1 (x, 0) = u 1 (x) -u il 1 (x, 0). The existence of a global in time classical solution q il (x, τ ) may be obtained by using (3.31), the corresponding initial condition (3.32) and the symmetry and positive definiteness of η vv and S. The exponential decay for w il 1 (x, τ ) is obtained by using the exponential decay of w il 0 (x, τ ) [START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System with Stiff Source Terms[END_REF]. In summary, u 0 and q il 0 are first obtained whereas q 0 = 0 and u il 0 = 0. Then q 1 is directly evaluated as well as u il 1 and this finally yields the required initial conditions u 1 (0) and q il 1 (0) in order to obtain u 1 and q il 1 .

Approximate equations

We investigate in this section the approximate equation satisfied by the second-order truncated approximation w a ε (x, t) = w 0 (x, t)+εw 1 (x, t)+w il 0 (x, t/ε)+εw il 1 (x, t/ε). Defining the residual operator R(f) for any function f(x, t) twice differentiable in space and once in time by

(3.34) R(f) = ∂ t f + j∈D A j (f)∂ j f -ε i,j∈D ∂ i B i,j (f)∂ j f + 1 ε L(f)f -εd (f, ∂ x f), then R a ε = R(w a ε
) is a measure of the accuracy of the approximated solution w a ε . In order to estimate the composite residual R a ε , we first estimate the outer residual R out ε = R w 0 (x, t)+εw 1 (x, t) associated with the outer expansion w 0 (x, t) + εw 1 (x, t), and next estimate the difference of residuals

δR ε = R a ε -R out ε associated with initial-layer corrector w il 0 + εw il 1 (3.35) δR ε = R a ε -R out ε = R(w 0 + εw 1 + w il 0 + εw il 1 ) -R(w 0 + εw 1 )
. The outer residual R out ε associated with the standard time dynamics is investigated in the following proposition.

Proposition 3.1. The outer residual R out ε = R w 0 (x, t) + εw 1 (x, t) can be written R out ε = εG ε , where G ε = ∂ t w 1 + j∈D 1 0 ∂ w A j (w 0 + εαw 1 ) dα w 1 ∂ j w 0 (3.36) + j∈D A j (w 0 + εw 1 )∂ j w 1 + 1 0 ∂ w L(w 0 + εαw 1 ) dα w 1 w 1 - i,j∈D ∂ i B ij (w 0 + εw 1 )∂ j (w 0 + εw 1 ) -d (w 0 + εw 1 , ∂ x w 0 + ε∂ x w 1 ).

Moreover we have εG

ε = εG 1 + ε 2 G 2ε where G 1 =∂ t w 1 + j∈D ∂ w A j (w 0 )w 1 ∂ j w 0 + j∈D A j (w 0 )∂ j w 1 (3.37) + ∂ w L(w 0 )w 1 w 1 - i,j∈D ∂ i B ij (w 0 )∂ j w 0 -d (w 0 , ∂ x w 0 ), G 2ε = j∈D 1 0 1 0 ∂ 2 w A j (w 0 + εαβw 1 ) αdαdβ w 1 w 1 ∂ j w 0 + j∈D 1 0 ∂ w A j (w 0 + εαw 1 ) dα w 1 ∂ j w 1 + 1 0 1 0 ∂ 2 w L(w 0 + εαβw 1 ) αdαdβ w 1 w 1 w 1 + i,j∈D ∂ i 1 0 ∂ w B ij (w 0 + εαw 1 ) dα w 1 ∂ j w 0 + i,j∈D ∂ i B j (w 0 + εw 1 )∂ j w 1 + i,j∈D 1 0 ∂ w M ij (w 0 + εαw 1 ) dα w 1 ∂ i w 0 ∂ j w 0 + i,j∈D M ij (w 0 + εw 1 ) ε∂ i w 1 ∂ j w 1 + ∂ i w 0 ∂ j w 1 + ∂ i w 1 ∂ j w 0 .
In addition, letting

G 1 = (G 1e , G 1r ) t we have G 1e = 0 and G 1 ∈ E ⊥ .
Proof. These relations are established after some algebra by using the governing equations for w 0 and w 1 and calculus identities in the form

A i (w 0 + εw 1 ) -A i (w 0 ) = ε 1 0 ∂ w A i (w 0 + εαw 1 ) dα w 1 .
Such identities are notably valid since the open set O w is convex so that any segment [w 0 , w 0 + εw 1 ] lies in O w whenever w 0 ∈ O w and w 0 + εw 1 ∈ O w . The expression of G ε is first obtained and consequently that of G 1 by letting ε → 0. The slow component of G 1 is then found to be zero by using the governing equation (3.11) of u 1 , keeping in mind that d e = 0.

We now investigate the difference residual δR ε that is naturally expanded around τ = 0 in order to use the governing equations of w il 0 and w il 1 . Proposition 3.2. The residual difference δR ε may be written in the form δR ε = εF ε where the expression of F ε is detailed in Appendix B. Moreover F ε decreases exponentially towards zero as τ → ∞ when w il 0 and w il 1 decrease exponentially to zero.

Proof. The proof is essentially a lengthy calculation summarized in Appendix B. The difference residual δR ε is naturally obtained in a form involving the initial-layer correctors w il 0 and w il 1 or their derivatives in each contributing term and this yields the exponential convergence towards zero.

It will be convenient in the following to consider R out ε as a function of t and δR ε as a function of τ with τ = t/ε. Combining the expressions of G ε and F ε we may now state the following result on the approximate system. Proposition 3.3. The second-order approximation w a ε satisfies the system of equations

∂ t w a ε + j∈D A j (w a ε )∂ j w a ε -ε i,j∈D ∂ i B i,j (w a ε )∂ j w a ε + 1 ε L(w a ε )w a ε (3.38) -εd (w a ε , ∂ x w a ε ) = ε(G ε + F ε ),
where G ε and F ε are detailed in Propositions 3.1 and 3.2 Proposition 3.3 shows that whenever the outer coefficients w 0 and w 1 are estimated as well as the initial-layer correctors w il 0 and w il 1 , then the residuals R out

ε and δR ε = R a ε -R out ε
are also estimated, in such a way that w a ε satisfies the approximated system (3.38) with a O(ε) residual.

The outer second-order fast-component corrector

In order to establish that the truncated expansion w a ε yields a second-order accurate approximation of w ε , we need to consider modified expansions. More specifically, due to the stiff sources, the outer residual G ε in the approximated equation (3.38) has a first order contribution εG 1 . Even though using the truncated expansion w a ε leads to a O(ε 3/2 ) accuracy estimate [START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System with Stiff Source Terms[END_REF], new modified expansions are required in order to obtain an O(ε 2 ) accuracy estimate without extra regularity assumptions.

Since G 1 belongs to the fast manifold G 1 ∈ E ⊥ from Proposition 3.1, a natural idea is to absorb this residual by using the source term. This precisely leads to add the second-order fast-component q 2 corrector term in the outer expansion as established in the following proposition. Proposition 3.4. Denote by w 2 = (u 2 , q 2 ) t the second-order outer corrector expansion coefficient and by w ′ 2 = (0, q 2 ) t its projection on the fast manifold. Then we have the identity

(3.39) L(w 0 )w 2 = L(w 0 )w ′ 2 = -G 1 .
Proof. The ε 1 relation deduced from (2.7) is in the form

∂ t w 1 + j∈D ∂ w A j (w 0 )w 1 ∂ j w 0 + j∈D A j (w 0 )∂ j w 1 + ∂ w L(w 0 )w 1 w 1 (3.40) - i,j∈D ∂ i B ij (w 0 )∂ j w 0 -d (w 0 , ∂ x w 0 ) + L(w 0 )w 2 = 0.
The identity (3.39) is then a consequence of (3.40

w ′a ε = w a ε + ε 2 w ′ 2 ,
that appears to be the simplest modified expansion with a second order accurate outer residual. In particular, only the fast part q 2 of the second-order outer correction w 2 is required in order to cancel the residual εG 1 and this will be sufficient in order to establish a second-order accuracy for w a ε as well as w ′a ε . Such a procedure appears to be natural since the orders are 'mixed' so that we first get q 0 = 0 at the order ε -1 , then we get u 0 and q 1 at the order ε 0 and next u 1 and q 2 at the order ε 1 simultaneously, so that the 'orders' are not clearly separated and we need to consider q 2 even if we only want to estimate the error arising from q 0 , u 0 , q 1 , and u 1 . Another point of view is that some of the natural norms arising in a priori estimates-once rewritten in terms of the quasinormal variable-are in the form |w -

w ⋆ | 2 k + 1 ε |q| 2 k-1
where w ⋆ denotes an equilibrium state and | • | k the norm in the Sobolev space H k (R d ) [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation: Local convergence results, Nonlinear Analysis[END_REF][START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF], in such a way that, loosely speaking, more accuracy is required on the fast component because of the 1/ε factor.

Improved approximate equations

We investigate in this section the residual associated with the modified truncated approximation

w ′a ε = w a ε + ε 2 w ′ 2 .
This residual is denoted by R ′a ε = R(w ′a ε ) and is a measure of the accuracy of the modified approximated solution w ′a ε .

Proposition 3.5. Keeping the notation of Propositions 3.1 and 3.2 the modified second-order approximation w ′a ε satisfies the system

∂ t w ′a ε + j∈D A j (w ′a ε )∂ j w ′a ε -ε i,j∈D ∂ i B i,j (w ′a ε )∂ j w ′a ε + 1 ε L(w ′a ε )w ′a ε (3.42) -εd (w ′a ε , ∂ x w ′a ε ) = ε 2 G ′ 2ε + ε 3 G ′ 3ε + εF ′ ε , with a residual R(w ′a ε ) = ε 2 G ′ 2ε + ε 3 G ′ 3ε + εF ′ ε split into G ′ 2ε = G 2ε + ∂ t w ′ 2 + i∈D A i (w ′a ε )∂ i w ′ 2 + 1 0 ∂ w L w 0 + w il 0 + αε(w 1 + w il 1 + εw ′ 2 ) dα (w 1 + w il 1 + εw ′ 2 ) w ′ 2 + i∈D 1 0 ∂ w A i (w a ε + αε 2 w ′ 2 ) dα w ′ 2 ∂ i w a ε + 1 0 ∂ w L(w a ε + αε 2 w ′ 2 ) dα w ′ 2 w 1 -ε i,j∈D 1 0 ∂ w M ij (w a ε + αε 2 w ′ 2 ) dα w ′ 2 ∂ i w a ε ∂ j w a ε -ε i,j∈D M ij (w ′a ε ) ε 2 ∂ i w ′ 2 ∂ j w ′ 2 + ∂ i w a ε ∂ j w ′ 2 + ∂ i w ′ 2 ∂ j w a ε -ε i,j∈D ∂ w B ij (w ′a ε )∂ i w ′a ε ∂ j w ′ 2 -ε i,j∈D ∂ i 1 0 ∂ w B ij (w a ε + αε 2 w ′ 2 ) dα w ′ 2 ∂ j w a ε , G ′ 3ε = - i,j∈D B ij (w ′a ε )∂ i ∂ j w ′ 2 , F ′ ε = F ε + 1 0 ∂ w L(w 0 + αw il 0 ) dαw il 0 w ′ 2 + 1 0 ∂ w L(w a ε + αε 2 w ′ 2 ) dα w ′ 2 (w il 0 + εw il 1 ).
Proof. The residual R(w ′a ε ) is first written in the form

R(w ′a ε ) = R(w a ε ) + ε 2 ∂ t w ′ 2 + i∈D A i (w ′a ε )∂ i w ′ 2 -ε i,j∈D ∂ i B ij (w ′a ε )∂ j w ′ 2 + ε L(w ′a ε ) -L(w 0 + w il 0 ) w ′ 2 + L(w 0 + w il 0 ) -L(w 0 ) w ′ 2 + L(w 0 )w ′ 2 + 1 ε L(w ′a ε ) -L(w a ε ) w a ε + i∈D A i (w ′a ε ) -A i (w a ε ) ∂ i w a ε -ε i,j∈D ∂ i B ij (w ′a ε ) -B ij (w a ε ) ∂ j w a ε -ε d (w ′a ε , ∂ x w ′a ε ) -d (w a ε , ∂ x w a ε ) .
The first order term εL(w 0 )w ′ 2 then compensates with εG 1 arising from εG ε in R(w a ε ) and this yields the decomposition of R(w ′a ε ).

The new outer residual has been decomposed into

ε 2 G ′ 2ε + ε 3 G ′ 3ε
where G ′ 3ε is less regular and need to be considered separately. We finally note that the initial condition associated with the new approximated solution is now in the form

(3.43) w ′a ε (0) = w 0 + ε w 1 + ε 2 w ′ 2 (0),
where w ′ 2 (0) = 0, q 2 (0) t and q 2 (0) is directly obtained from (3.39).

Mathematical analysis

We establish in this section the existence of strong solutions to the Cauchy problem for system (2.7).

We further rigorously establish the accuracy of the second-order truncated approximations w a ε and w ′a ε . To this aim, we first establish the existence of appropriate outer expansion coefficients w 0 and w 1 and initial-layer correctors w il 0 and w il 1 and discuss standard normal forms that uncouple hyperbolic and parabolic variables. We then use a convergence-stability lemma in order to establish existence of solutions and the accuracy of the approximate solutions.

Existence of approximate solutions

|∂ k φ| 2 L 2 = |α|=k k! α! R d |∂ α φ| 2 dx,
where k!/α! are the multinomial coefficients with a similar definition for vector functions. For any map φ : R d ×[0, t ] → R n where t > 0 is positive and for any t ∈ [0, t ], we denote by φ(t) or φ(•, t) the partial map x → φ(x, t) defined over R d with a similar convention for functions that depends on the fast time τ . We collect in the following Lemma various useful properties of Sobolev spaces [START_REF] Nirenberg | On Elliptic Partial Differential Equations[END_REF][START_REF] Kawashima | System of a Hyperbolic-Parabolic Composite Type, with applications to the equations of Magnetohydrodynamics[END_REF][START_REF] Hudjaev | On the Cauchy Problem for Composite Systems of Nonlinear Differential Equations[END_REF]. (ii) Assume that l ≥ l 0 + 1, f ∈ H l and g ∈ H l-1 . Then for all multi-indices |α| ≤ l we have [f, ∂ α ]g = f∂ α g -∂ α (fg) ∈ L 2 and there exists a constant c 0 depending on l and d such that

|f∂ α g -∂ α (fg)| 0 ≤ c 0 |f| l |g| |α|-1 .
(iii) The following nonlinear estimate holds for k ≥ 1 and

φ ∈ H k ∩ L ∞ (4.2) |f (φ) -f (0)| k ≤ c 0 |f | C k (O φ ) (1 + |φ| k-1 L ∞ ) |φ| k ,
where O φ denotes a ball of radius strictly larger than |φ| L ∞ and c 0 a constant depending on k and d.

We denote by w ⋆ an equilibrium state w ⋆ = (u ⋆ , 0) with q ⋆ = 0 and look for solutions w such that w

ε -w ⋆ ∈ C 0 ([0, t * ], H l+2 ) and similarly w a ε -w ⋆ ∈ C 0 ([0, t * ], H l+2 ).
There is a natural shift of two order derivatives, however, between the governing equations of u 0 and u 1 , since there are secondorder derivatives of u 0 in the governing equation of u 1 . We must thus have u

0 -u ⋆ ∈ C 0 ([0, t * ], H l+4 ) in order to obtain that u 1 ∈ C 0 ([0, t * ], H l+2 ) and w a ε -w ⋆ ∈ C 0 ([0, t * ], H l+2 ).
The initial data is naturally expanded in powers of ε

(4.3) w ε = w 0 + ε w 1 + O(ε 2 ),
and the coefficients are correspondingly assumed such that w 0w ⋆ ∈ H l+4 , q 0 = 0 and w 1 ∈ H l+2 . In order to establish the existence and properties of approximated solutions the condition (u 5 ) will be especially useful. Moreover, we assume that the relative entropy η il is level bounded in the following sense.

(L) The relative entropy η il is such that lim |q|→∞ η il (u, q) = ∞ uniformly with respect to u on compact sets of O u .

Heuristically, this condition means that the entropy η is infinite at finite boundaries of the equilibrium set O u ×{0}. Such an assumption is notably used in convex analysis [33] and holds for the two-temperature fluids presented in [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF][START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation: Local convergence results, Nonlinear Analysis[END_REF]. Assuming (u 1 )-(u 5 ) and (L) we now rigorously establish the existence of the multitime expansion coefficients investigated formally in Section 3. Proposition 4.2. Let l ≥ l 0 + 1 and l 0 = [d/2] + 1 be integers and b > 0 be given. Assume that the initial data is such that w 0w ⋆ ∈ H l+4 and

w 1 ∈ H l+2 with | w 0 -w ⋆ | l+4 + | w 1 | l+2 < b. Assume that (u 1 )-(u 5 ) and (L) hold, that w 0 ∈ O w = O u × (-q, q ) nr ⊂ O w ,
where O u is a bounded convex open set with its closure in O u and q > 0 is a given initial bound for fast variable components. Then there exist a time t * > 0, outer expansion coefficients w 0 (t) and w 1 (t) defined over [0, t * ], as well as initial-layer correctors w il 0 (τ ) and w il 1 (τ ) defined over [0, ∞) and decreasing exponentially as τ → ∞, with regularity

(4.4) w 0 -w ⋆ ∈ C 0 ([0, t * ], H l+4 ), w 1 ∈ C 0 ([0, t * ], H l+2 ), (4.5)
w il 0 ∈ C 0 ([0, +∞), H l+4 ), w il 1 ∈ C 0 ([0, +∞), H l+2 ),
such that the outer expansion equations (3.7), (3.8), and (3.11), the inner expansion equations (3.25) and (3.29), and the matching conditions (3.4) are satisfied. These coefficients are uniquely defined and there exists positive constants c and δ depending on O w , O u , q, b and the system coefficients A j , B ij , L and η such that

(4.6) |w 0 (t) -w ⋆ | l+4 + |w 1 (t) | l+2 ≤ c, 0 ≤ t ≤ t * , (4.7
)

|w il 0 (τ ) | l+4 + |w il 1 (τ ) | l+2 ≤ c exp(-δτ ), 0 ≤ τ < ∞.
Moreover, there exists O w 0 = O u0 × (-q 0 , q 0 ) nr ⊂ O w with O u ⊂ O u0 ⊂ O u0 ⊂ O u and 0 < q < q 0 , as well as ε ∈ (0, 1], only depending on O w , O u , q, b and the system coefficients A j , B ij , L and η, such that for any ε ∈ [0, ε], the functions w 0 , w 0 + εw 1 , and w 0 + εw 1 + ε 2 w ′ 2 , as well as w 0 + w il 0 , w 0 + w il 0 + ε(w 1 + w il 1 ), and w 0 + w il 0 + ε(w

1 + w il 1 ) + ε 2 w ′ 2 have their value in O w 0 ⊂ O w 0 ⊂ O w .
Corollary 4.3. Letting w 0 = (u 0 , q 0 ) t , w 1 = (u 1 , q 1 ) t , w il 0 = (u il 0 , q il 0 ) t , and w il 1 = (u il 1 , q il 1 ) t , there exist u 0 (t), q 0 (t), u 1 (t) and q 1 (t) with

u 0 -u ⋆ ∈ C 0 ([0, t * ], H l+4 ) ∩ C 1 ([0, t * ], H l+3 ), q 0 = 0, (4.8) u 1 ∈ C 0 ([0, t * ], H l+2 ) ∩ C 1 ([0, t * ], H l+1 ), (4.9) q 1 ∈ C 0 ([0, t * ], H l+3 ) ∩ C 1 ([0, t * ], H l+2 ), (4.10)
and initial-layer coefficients u il 0 (τ ), q il 0 (τ ), u il 1 (τ ), q il 1 (τ ) with Proof. In the following, a > 1 and δ < 1 denote respectively large and small generic constants only depending on O w , O u , q, and the system coefficients A j , B ij , L and η, whereas c > 1 denotes a generic constant that also depends on b. In the proof, we will successively obtain the existence and properties of w 0 = (u 0 , q 0 ) t , the existence and exponential decrease of w il 0 = (u il 0 , q il 0 ) t , the integrability properties of w il 0 , the existence and properties of w 1 = (u 1 , q 1 ) t , the existence, exponential decrease and integrability of w il 1 = (u il 1 , q il 1 ) t .

u il 0 = 0, q il 0 ∈ C 1 ([0, ∞), H l+4 ), (4.11) u il 1 ∈ C 1 ([0, ∞), H l+3 ), q il 1 ∈ C 1 ([0, ∞), H l+2 ), (4.
Step 1. The outer coefficient w 0 . From (3.7) we have q 0 = 0 and from (3.9) and Proposition 2.4, the zeroth order slow variable u 0 satisfies a symmetrizable hyperbolic system of conservation laws. More over u 0 ∈ O u , u 0 (•, 0) -u ⋆ = u 0 -u ⋆ ∈ H l+4 , and | u 0 -u ⋆ | l+4 < b. The initial value u 0 has its value in O u and we arbitrary select an open bounded set O u0 that contains the closure of O u and with O u0 ⊂ O u . According to the local-in-time existence theory of initial-value problems for symmetrizable hyperbolic system [START_REF] Godunov | An interesting class of quasilinear system[END_REF][START_REF] Friedrichs | System of conservation laws with a convex extension[END_REF][START_REF] Majda | Compressible Fluid Flow and System of Conservation Laws in Several Space Variables[END_REF][START_REF] Kawashima | System of a Hyperbolic-Parabolic Composite Type, with applications to the equations of Magnetohydrodynamics[END_REF][START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF][START_REF] Dafermos | Hyperbolic Conservation Laws in Continuum Physics[END_REF][START_REF] Kawashima | Dissipative structure and entropy for hyperbolic systems of conservation laws[END_REF], there exists a unique solution u 0 = u 0 (x, t) with initial value u 0 (x, 0) = u 0 (x) defined over [0, t * ] where t * > 0 with values in O u0 . This solution also satisfies the estimates

|u 0 (t) -u ⋆ | l+4 ≤ c loc | u 0 -u ⋆ | l+4 ≤ c, 0 ≤ t ≤ t * ,
where c loc only depend on O u0 , b and the system coefficients. This establishes the existence of w 0 with the regularity properties (4.4).

Step 2. Existence of the initial-layer corrector w il 0 . From (3.26) we have u il 0 = 0 and the ordinary differential equation (3.28) governing q il 0 yields, for any fixed x ∈ R d , a unique local solution q il 0 (x, τ ) with respect to the fast time τ starting from q 0 (x) ∈ (q, q) nr . In order to establish that this solution q 0 (x, τ ) is defined over the whole fast time interval τ ∈ [0, ∞) and decreases exponentially towards zero as τ → ∞, we use the relative entropy η il . This entropy η il is nonnegative and locally behaves quadratically η il = Hq, q where H = 1 0 η -1 vv u, αq)αdα as established in Lemma 2.5. Letting ϕ(τ ) = η il u 0 (x), q il 0 (x, τ ) and using ∂ τ ϕ = ∂ v η∂ q v∂ τ q il 0 , q = (∂ v η) t , and ∂ q v∂ τ q il 0 = -Sq, we obtain that ∂ τ ϕ + q, Sq = 0. The modified entropy ϕ is thus decreasing and is a Lyapunov function for the q il 0 differential equation. By integration, leaving implicit the dependence on x for the sake of notational simplicity, we get (4.13) η il ( u 0 , q il 0 (τ ) + τ 0 q il 0 (τ ′ ), S u( u 0 , q il 0 (τ ′ )) q il 0 (τ ′ ) dτ ′ = η il ( u 0 , q 0 ).

Letting M = sup{ η( u, q); u ∈ O u , q ∈ (q, q ) nr }, that is finite since O u , ×(q, q ) nr has compact closure in O w , and using that η il is level bounded (L), there exists q 0 such that q / ∈ (-q 0 , q 0 ) nr implies η( u 0 , q) > M . Therefore, (4.13) implies that whenever u 0 ∈ O u and q 0 ∈ (q, q ) nr , we must have q il 0 (τ ) ∈ (-q 0 , q 0 ) nr for any τ . All trajectories therefore remain uniformly bounded and thus exist for all positive τ . Letting then O ′ w 0 = O u0 ×(-q 0 , q 0 ) nr we have O ′ w 0 = O u0 ×[-q 0 , q 0 ] nr ⊂ O w and all trajectories are such that u 0 (t) ∈ O u0 and q il 0 (τ ) ∈ (-q 0 , q 0 ) nr and w 0 (t) + w il 0 (τ ) ∈ O ′ w 0 for any t ∈ [0, t * ] and τ ∈ [0, ∞). On the other hand, on the compact set O ′ w 0 of O w , that includes all trajectories starting from O w , there exists δ > 0 such that 2δϕ ≤ q, Sq and δ q, q ≤ ϕ ≤ (1/δ) q, q since ϕ, q, q , and Sq, q all behave quadratically with respect to q. From the inequality 2δϕ ≤ q, Sq we get that ∂ τ ϕ + 2δϕ ≤ 0 and by integration ϕ(τ ) ≤ η il ( u 0 , q 0 ) exp(-2δτ ). Further using δ q, q ≤ ϕ we obtain |q il 0 (x, τ )| ≤ a| q 0 (x)| exp(-δτ ) and (4.14) δτ ). This yields the exponential decrease towards zero of the initial-layer corrector uniformly in x ∈ R d and uniformly for the initial conditions u 0 ∈ O u and q 0 ∈ (q, q ) nr .

|q il 0 (τ )| L ∞ ≤ a exp(-
Step 3. Integrability of the initial-layer corrector w il 0 . We first note that by using mollifiers, we may assume that u 0 -u ⋆ and q 0 are smooth functions with compact support. In this situation, q il 0 is also smooth with compact support so that it is sufficient to establish a priori estimates of q il 0 in Sobolev spaces in this particular situation.

Since q 0 is in L 2 and |q il 0 (x, τ )| ≤ a| q 0 (x)| exp(-δτ ), we first obtain that q il 0 is in L 2 and that q il 0 (τ ) 0 ≤ a | q 0 | 0 exp(-δτ ) ≤ c exp(-δτ ). We now establish by induction that for any 0 ≤ k ≤ l + 4 (4.15)

|q il 0 | k ≤ c exp(-δ k τ ),
for some decreasing constants δ k < δ k-1 < • • • < δ depending on O w , O u , q, and the system coefficients A j , B ij , L and η. This estimate is already established for k = 0 and we now consider the situation where 1 ≤ k ≤ l + 4. Letting α be a multiindex with |α| = k, we may write α = e j + α ′ for some j ∈ D where e 1 , . . . , e d denotes the canonical basis of R d and |α ′ | = k -1. Applying ∂ j to the equation ∂ τ q il 0 = Q r ( u 0 , q il 0 ) we first get that ∂ τ ∂ j q il 0 = ∂ q Q r ∂ j q il 0 + ∂ u Q r ∂ j u 0 . Further applying the derivation operator ∂ α ′ we obtain (4.16)

∂ τ ∂ α q il 0 = ∂ q Q r ∂ α q il 0 + [∂ α ′ , ∂ q Q r ]∂ j q il 0 + ∂ α ′ (∂ u Q r ∂ j u 0 ).
The partial derivative ∂ q Q r may be decomposed into

∂ q Q r ( u 0 , q il 0 ) = -η vv S -∂ q (η vv S)q = -η vv S -(η vv S -η vv S) -∂ q (η vv S)q,
where η vv = η vv ( u 0 , 0) and S = S u( u 0 , 0) . Letting

J ( u 0 , q) = η -1 vv 1 0 ∂ q (η vv S)( u 0 , αq) dα + ∂ q (η vv S) ,
we have ∂ q Q r ( u 0 , q il 0 ) = -η vv S -η vv J q and ∂ u Q (u, q) = -∂ u (η vv S) q. Multiplying equation (4.16) by η -1 vv = η -1 vv ( u 0 ), that does not depend on τ , and next scalarly by ∂ α q il 0 , it is obtained that

∂ τ 1 2 ∂ α q il 0 ,η -1 vv ∂ α q il 0 + ∂ α q il 0 , S∂ α q il 0 + ∂ α q il 0 , J q il 0 ∂ α q il 0 = + ∂ α q il 0 , η vv [∂ α ′ , ∂ q Q r ]∂ j q il 0 -∂ α q il 0 , η vv ∂ α ′ ∂ u (η vv S) q il 0 ∂ j u 0 .
Letting ϕ = ∂ α q il 0 , η -1 vv ∂ α q il 0 and using the inequality 4δ ′ ϕ ≤ ∂ α q il 0 , S∂ α q il 0 as well as the inequality

| ∂ α q il 0 J q il 0 ∂ α q il 0 | ≤ a|q il 0 (τ )| L ∞ ϕ both valid over O ′
w 0 we obtain after some algebra

∂ τ 1 2 ϕ + 3δ ′ -a|q il 0 (τ )| L ∞ ϕ ≤ c [∂ α ′ , ∂ q Q r ]∂ j q il 0 2 0 + c ∂ u (η vv S) q il 0 ∂ j u 0 2 k-1 .
Using the induction assumption, Leibnitz formula, and the commutator estimates then yields

∂ τ 1 2 ϕ + 3δ ′ -a|q il 0 (τ )| L ∞ ϕ ≤ c exp(-2δ k-1 τ ).
Using (4.14) there exists τ 0 such that a|q il 0 (τ )| L ∞ ≤ δ ′ for τ ≥ τ 0 and then (4.17)

∂ τ 1 2 ϕ + 2δ ′ ϕ ≤ c exp(-2δ k-1 τ ).
Summing for |α| = k, integrating over [τ 0 , ∞), and selecting δ k ≤ δ ′ and δ k < δ k-1 we obtain that q il 0 (τ ) k ≤ c q il 0 (τ 0 ) k exp(-δ k τ ) + c exp(-δ k τ ) for τ ≥ τ 0 . On the other hand, we may also directly write that on the interval [0, τ 0 ] we have ∂ τ 1 2 ϕ ≤ aϕ + c, so that applying Gronwall lemma yields that q il 0 (τ ) k ≤ c over [0, τ 0 ]. Combining both estimates over [0, τ 0 ] and [τ 0 , +∞) then completes the proof of (4.15).

Step 4. Existence of the outer corrector w 1 . Once u 0 is known, q 1 may directly be evaluated by using (3.10). Taking into account the regularity properties of u 0 , and since the expression of q 1 only involves first order spatial derivatives of u 0 , we obtain that

q 1 ∈ C 0 ([0, t * ], H l+3 ) ∩ C 1 ([0, t * ], H l+2 ).
Similarly, once u 0 and q il 0 are known, the initial-layer corrector u il 1 is directly obtained from (3.33) in the form

u il 1 (x, τ ) = j∈D ∂ j ∞ τ 1 0 ∂ q f j u, αq il 0 (τ ′ ) dα q il 0 (τ ′ )dτ ′ .
Using the properties of q il 0 established in Step 2, in particular the uniform exponential decrease, it is obtained that

u il 1 ∈ C 1 ([0, ∞), H l+3 ) and |u il 1 (τ )| l+3 ≤ c exp(-δτ ). Once u il
1 is known, we may evaluate u il 1 (0) ∈ H l+3 , and thus u 1 (0) = u 1 -u il 1 (0) so that u 1 (0) ∈ H l+2 since u 1 ∈ H l+2 and u il 1 (0) ∈ H l+3 . The first order slow variable u 1 also satisfies a semilinear and symmetrizable hyperbolic system with nonhomogeneous terms that are functions of u 0 , ∂ x u 0 and ∂ 2

x u 0 and are C 0 ([0, t * ], H l+2 ). According to the local-in-time existence theory for symmetrizable hyperbolic system, there exists a unique solution u 1 = u 1 (x, t) with initial condition u 1 (0) = u 1 -u il 1 (0) for a positive time interval that may be chosen to be [0, t * ], eventually after a change of notation. This solution u 1 also satisfies the estimate

|u 1 (t)| l+2 ≤ c|u 1 (0)| l+2 ≤ c, 0 ≤ t ≤ t * ,
and the existence of w 0 and w 1 with regularity (4.4) is established.

Step 5. Existence and integrability of the initial-layer corrector w il 1 . The global existence of q il 1 (x, τ ) in C 1 ([0, ∞), R nr ) follows from the general theory of linear differential equations since q il 1 satisfies the linear equation with bounded coefficients

∂ τ q il 1 = ∂ q Q r (w p 0 + w il 0 )q il 1 + r il 1
, where

r il 1 = ∂ u Q r (w p 0 + w il 0 ) -∂ u Q r (w p 0 ) u p 1 + ∂ u Q r (w p 0 + w il 0 )u il 1 + ∂ q Q r (w p 0 + w il 0 ) -∂ q Q r (w p 0 ) q p 1 + C(w p 0 + w il 0 ) -C(w p 0 ).
The polynomial factor in τ of u p 1 is in particular compensated by the exponential decrease towards zero of u il 0 . Using the τ -independent metrics already used for q il 0 as well as the exponential decrease of r il 1 , the exponential decrease of q il 1 is next established. The integrability properties of q il 1 are consequences of that of r il 1 obtained from (4.15) and proceeding as in Step 3 yields that

q il 1 ∈ C 1 ([0, ∞), H l+2 ) and |q il 1 (τ )| l+2 ≤ c exp(-δτ ).
Finally, the values of w 0 and w il 0 remain in the compact set O ′ w 0 = O u0 ×[-q 0 , q 0 ] nr and we may find

O w 0 with O ′ w 0 ⊂ O w 0 ⊂ O w 0 ⊂ O w .
On the other hand, the values of w 1 , w il 1 , and

w ′ 2 remain in a compact set of R n since |w 1 | L ∞ + |w il 1 | L ∞ + |w ′ 2 | L ∞ ≤ c.
As a consequence, for ε small enough,

w 0 + εw 1 , w 0 + εw 1 + ε 2 w ′ 2 , w 0 + w il 0 + ε(w 1 + w il 1 )
, and

w 0 + w il 0 + ε(w 1 + w il 1 ) + ε 2 w ′ 2 all have their value in the open set O w 0 and O w 0 ⊂ O w 0 ⊂ O w for any ε ∈ [0, ε].
In the following proposition we investigate the approximated system satisfied by the second-order truncated expansions w a ε and w ′a ε . Proposition 4.4. Keeping the assumption and notation of Proposition 4.2, for any ε ∈ (0, ε], the approximated solution w a ε = w 0 + εw 1 + w il 0 + εw il 0 is such that w a εw ⋆ ∈ C 0 ([0, t * ], H l+2 ). The residuals G ε , G 1 , G 2ε , and F ε defined in Propositions 3.1 and 3.2 are such that

G ε , G 2ε ∈ C 0 ([0, t * ], H l ), G 1 ∈ C 0 ([0, t * ], H l+1 ), F ε ∈ C 0 ([0, t * /ε], H l )
, and F ε is decreasing exponentially with respect to the fast time. Moreover, there exist c and δ independent of ε with

(4.18) |G ε (t)| l + |G 2ε (t)| l + |G 1 (t)| l+1 ≤ c 0 ≤ t ≤ t * , (4.19) |F ε (τ ) | l ≤ c exp(-δτ ), 0 ≤ τ ≤ t * /ε.
Similarly, the second-order projected outer corrector w ′ 2 = (0, q 2 ) t is such that w ′ 2 ∈ C 0 ([0, t * ], H l+1 ) and the modified approximate solution

w ′a ε = w a ε + ε 2 w ′ 2 is such that w ′a ε -w ⋆ ∈ C 0 ([0, t * ], H l+1 ). Finally, the residuals G ′ ε , G ′ 2ε , G ′ 3ε , F ′ ε defined in Proposition 3.5 are such that G ′ 2ε ∈ C 0 ([0, t * ], H l ), G ′ 3ε ∈ C 0 ([0, t * ], H l-1 ), F ′ ε ∈ C 0 ([0, t * /ε], H l
) and there exist constants c and δ independent of ε such that Similarly, the regularity properties of w ′a ε and G ′ 2ε , G ′ 3ε , and F ′ ε and the corresponding estimates are also consequences of Proposition 4.2 and Proposition 3.5. Finally, the estimating constants only depend on O w , O u , q, b and the system coefficients A j , B ij , L and η.

(4.20) |G ′ 2ε (t)| l + |G ′ 3ε (t)| l-1 ≤ c, 0 ≤ t ≤ t * , (4.21) |F ′ ε (τ ) | l ≤ c exp(-δτ ), 0 ≤ τ ≤ t * /ε.

Normal variable

In order to establish the existence of solutions w ε to (2.7) as well as the accuracy of w a ε and w ′a ε , it is necessary to have a local existence theorem. To this aim we have to put the system into a normal form that distinguishes between hyperbolic and parabolic variables [START_REF] Kawashima | System of a Hyperbolic-Parabolic Composite Type, with applications to the equations of Magnetohydrodynamics[END_REF][START_REF] Kawashima | Large-time Behavior of Solutions to Hyperbolic-Parabolic System of Conservations Laws and Applications[END_REF][START_REF] Kawashima | On the Normal Form of the Symmetric Hyperbolic-parabolic Systems Associated with the Conservation Laws[END_REF][START_REF] Giovangigli | Asymptotic stability of equilibrium states for multicomponent reactive flows[END_REF][START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF][START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation: Local convergence results, Nonlinear Analysis[END_REF][START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF]. A natural sufficient condition for system (1.1) to be rewritten into a norm form is Condition N introduced by Kawashima and Shizuta [START_REF] Kawashima | On the Normal Form of the Symmetric Hyperbolic-parabolic Systems Associated with the Conservation Laws[END_REF] and strengthened in [START_REF] Giovangigli | Multicomponent flow modeling[END_REF].

(N) The nullspace N B(u, ξ) of B(u, ξ) = i,j∈D B ij (∂ 2 u η) -1 ξ i ξ j does not depend on u ∈ O u and ξ ∈ Σ d-1 , and B ij (u)(∂ 2 u η) -1 N B(u, ξ) = 0, for i, j = 1, . . . , d.
Heuristically, the independence of N ( B) on ξ is an isotropy condition of diffusion processes and the independence on x is an ubiquitous condition of diffusion processes. For fluid models, Condition (N) generally results from natural constraints on transport coefficients associated with mass, momentum, or energy conservation and derived from the kinetic theory [START_REF] Chapman | The Mathematical Theory of Non-Uniform Gases[END_REF][START_REF] Ferziger | Mathematical Theory of Transport Processes in Gases[END_REF][START_REF] Giovangigli | Multicomponent flow modeling[END_REF]. Under Condition (N) it is then possible to describe all possible normal variables [START_REF] Giovangigli | Asymptotic stability of equilibrium states for multicomponent reactive flows[END_REF][START_REF] Giovangigli | Multicomponent flow modeling[END_REF]. Letting n i = dim N ( B) and n ii = n -n i , let P denote an arbitrary constant nonsingular matrix of dimension n such that its first n i columns span the nullspace N ( B), and define the auxiliary variables by u [START_REF] Giovangigli | Asymptotic stability of equilibrium states for multicomponent reactive flows[END_REF][START_REF] Giovangigli | Multicomponent flow modeling[END_REF]. The diffusion matrices corresponding to these auxiliary variables have nonzero coefficients only in the lower right block of size n ii = n -n i and normal forms are equivalently-and more easily-obtained from the v ′ symmetric equation [START_REF] Giovangigli | Asymptotic stability of equilibrium states for multicomponent reactive flows[END_REF][START_REF] Giovangigli | Multicomponent flow modeling[END_REF]. Lemma 4.5. Assume that the nullspace invariance property (N) is satisfied and denote by u ′ = P t u and v ′ = P -1 v the auxiliary variable. All normal variables of system (1.1) are in the form w

′ = P t u = (u ′ i , u ′ ii ) t and v ′ = P -1 v = (v ′ i , v ′ ii ) t where v = (∂ u η) t [
= φ i (u ′ i ), φ ii (v ′ ii )
t where φ i and φ ii are two C ∞ diffeomorphisms of R ni and R nii , respectively. The

component w i = φ i (u ′ i ) is then the hyperbolic variable and w ii = φ ii (v ′ ii )
is the parabolic variable. Assuming that Condition (N) holds, the structure of the normal form is investigated in the following proposition using the partition

R n = R ni × R nii .
Proposition 4.6. Let w = (w i , w ii ) t be a normal variable with a C ∞ diffeomorphism u → w. Letting u = u(w) in (1.1) and multiplying on the left by

(∂ w u) t ∂ 2 u η = (∂ w v) t yields (4.22) A 0 (w)w t + j∈D A j (w)∂ j w -ε i,j∈D ∂ i B ij (w)∂ j w = Q(w) ε + ε d(w, ∂ x w),
where the matrices

A 0 = (∂ w u) t ∂ 2 u η ∂ w u, A j = (∂ w u) t ∂ 2 u η A j ∂ w u, B ij = (∂ w u) t ∂ 2 u η B ij ∂ w u, the residual d(w, ∂ x w) = -i,j∈D ∂ i (∂ w u) t ∂ 2 u η B ij ∂ j u,
B t ij = B ji , i, j ∈ D, over O w . For any ξ ∈ Σ d-1 the diffusion matrix B ii,ii = i,j∈D B ii,ii ij ξ i ξ j is positive definite over O w .
The matrices A 0 , A j , B ij , and the source term Q have regularity C ∞ and the quadratic residual d may be written

(4.23) d(w, ∂ x w) = i,j∈D m ij (w)∂ i w∂ j w,
where m ij are third order tensors depending on w, have regularity C ∞ and only involve parabolic

components d(w, ∂ x w) = 0, i,j∈D m ii,ii,ii ij (w)∂ i w ii ∂ j w ii t .
We also assume in the following that the source term Q is in quasilinear form.

(w 4 ) The source term Q may be written Q(w) = -L(w) w where L is a symmetric positive semidefinite matrix of size n with a fixed nullspace N (L) = E and w → L(w) has regularity C ∞ .

We use the terminology quasilinear for such source terms Q since the matrix L has several invariant properties, being symmetric positive semi-definite with fixed range and nullspace. In particular, the linear subspace E is the slow manifold with respect to the normal variable w, E ⊥ is the corresponding fast manifold, and we denote by π the orthogonal projector onto E ⊥ . Source terms Q that are in quasilinear form with respect to the entropic variable v = (∂ u η) t are often encountered in mathematical physics [START_REF] Yong | An interesting class of partial differential equations[END_REF]. The assumption (w 4 ) concerning the source term of the normal form Q is thus natural in terms of properly chosen normal variables w [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF]. This is notably the case when the normal variable is such that πw = π(∂ v w)v as investigated in [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF] as well as for multitemperature nonequilibrium gas models [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF] and complex chemistry fluids [START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF]. The assumption (w 4 ) further implies extra properties already observed for the quasinormal variables like for instance Q(w), w = 0 if and only if Q(w) = 0 and if and only if w ∈ E. Incidentally, it is often the case that the source terms Q belongs to the parabolic subspace {0} × R nii but such a property will not be needed in the analysis. Finally, we introduce a commutation relation that simplifies a priori estimates [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF][START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation: Local convergence results, Nonlinear Analysis[END_REF][START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF].

(C) The commutation relation A 0 (w)π = πA 0 (w) holds where π is the orthogonal projector onto the fast manifold E ⊥ .

This property is equivalent to

A 0 E ⊂ E or A 0 E ⊥ ⊂ E ⊥ and naturally implies A -1 0 (w)L(w) = π t A -1
0 (w)L(w)π since L = Lπ = πL and π t = π. These relations simplify some commutator estimates for nonequilibrium variables [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF][START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF].

Since the change of variable u → w is a C ∞ diffeomorphism, we may use the induced C ∞ diffeomorphism w → w from O w onto O w and define naturally the second order approximations w a ε and w ′a ε in the normal variable as (4.24)

w a ε = w(w a ε ), w ′a ε = w(w ′a ε ).
The approximate equation satisfied by w ′a ε is investigated in the following proposition and an analogous equation for w a ε is similarly obtained. Proposition 4.7. Keeping the assumptions and notation of Proposition 4.2, for any ε ∈ (0, ε], the approximation w ′a ε is such that w ′a εw ⋆ ∈ C 0 ([0, t * ], H l+1 ) and satisfies the system

A 0 (w ′a ε )∂ t w ′a ε + j∈D A j (w ′a ε )∂ j w ′a ε -ε i,j∈D ∂ i B ij (w ′a ε )∂ j w ′a ε = Q(w ′a ε ) ε (4.25) + ε d(w ′a ε , ∂ x w ′a ε ) + ε 2 G ′ 2ε + ε 3 G ′ 3ε + εF ′ ε ,
where

G ′ 2ε = (∂ w v) t ∂ w u G ′ 2ε , G ′ 3ε = (∂ w v) t ∂ w u G ′ 3ε , and
F ′ ε = (∂ w v) t ∂ w u F ′ ε .
Moreover, there exists positive constants c and δ independent of ε such that

|G ′ 2ε (t)| l + |G ′ 3ε (t)| l-1 ≤ c, G ′ 3ε (t) ∈ {0}×R nii , 0 ≤ t ≤ t * , (4.26) |F ′ ε (τ )| l ≤ c exp(-δτ ), 0 ≤ τ ≤ t * /ε. (4.27)
Proof. The approximate equation in normal form (4.25) is obtained from the approximate equation in the quasinormal variable (3.42) after the change of variable w = w (w) and multiplication by (∂ w v) t ∂ w u = A 0 ∂ w w. The second-order derivative terms are also integrated by part in order to include the multiplying factor (∂ w v) t ∂ w u inside the derivation operators in such a way that the quadratic residual reads

d = (∂ w v) t ∂ w ud + i,j∈D (∂ w v) t ∂ w uB ij ∂ j w ′a ε . Proposition 4.4 then yields that w ′a ε -w ⋆ ∈ C 0 ([0, t * ], H l+1 ) and w a ε -w ⋆ ∈ C 0 ([0, t * ], H l+2 ) using that w → w is of regularity class C ∞ .
The integrability of the residuals also directly results from the expressions

G ′ 2ε = (∂ w v) t ∂ w u G ′ 2ε , G ′ 3ε = (∂ w v) t ∂ w u G ′ 3ε , and
F ′ ε = (∂ w v) t ∂ w u F ′ ε and Proposition 4.4.
Finally, the residual G ′ 3ε may be rewritten by using the definition of G ′ 3ε from Proposition 3.5 in the form

G ′ 3ε = -i,j∈D (∂ w v) t ∂ w uB ij ∂ i ∂ j w ′ 2 where all matrices (∂ w v) t , ∂ w u and B ij are evaluated at w ′a ε . Using then B ij = (∂ w u) -1 B ij ∂ w u we obtain that G ′ 3ε = - i,j∈D (∂ w v) t B ij ∂ w u∂ i ∂ j w ′ 2 , = - i,j∈D B ij ∂ w w∂ i ∂ j w ′ 2 ,
in such a way that G ′ 3ε ∈ {0}×R nii from the properties of the diffusion matrices B ij , i, j ∈ D. The initial conditions w ε , w 0 , and w 1 , are assumed to be such that

| w ε -w 0 -ε w 1 | l ≤ cε 2 in such a way that |w a ε (0) -w ε | l ≤ cε 2 and |w ′a ε (0) -w ε | l ≤ cε 2 .
Letting naturally w ε = w( w ε ) and w ′a ε (0) = w w ′a ε (0) , we thus have (4.28)

|w ′a ε (0) -w ε | l ≤ cε 2 .
Finally, we define the zeroth order outer coefficient w 0 and the limit at infinity w ⋆ in the normal variable as (4.29)

w 0 = w(w 0 ), w ⋆ = w(w ⋆ ).
Since

w ′a ε = w a ε + ε 2 w ′ 2
, where w a ε = w 0 + εw 1 + w il 0 + εw il 1 , we also obtain from Proposition 4.2 that there exists constants c and δ with δτ ) . Note that one may define similarly approximate solutions in the conservative variable u a ε = u(w a ε ), and u ′a ε = u(w ′a ε ) that also satisfy approximate equations analogous to (1.1) but only w a ε and w ′a ε are used in the following.

|w ′a ε -w 0 | l ≤ c|w ′a ε -w 0 | l ≤ c ε + exp(-

Existence of solutions

Let w be a normal variable and consider the Cauchy problem for the system of partial differential equations (4.22) with initial condition w ε where ε ∈ (0, ε] and ε is fixed with 0 < ε ≤ 1. The initial value w ε is assumed such that w εw ⋆ ∈ H l (R d ) where l ≥ l 0 + 1 and

w ε (x) ∈ O 0 , x ∈ R d , where O 0 is a bounded open set with O 0 ⊂ O 0 ⊂ O w . We further introduce a d 1 -neighborhood of O 0 with 0 < d 1 < d(O 0 , ∂O w ) defined by (4.30) O 1 = { w ∈ O w ; d(w, O 0 ) < d 1 }, in such a way that O 0 ⊂ O 1 ⊂ O 1 ⊂ O w and O 1 is bounded.
According to the local existence theory for hyperbolic-parabolic systems [START_REF] Hudjaev | On the Cauchy Problem for Composite Systems of Nonlinear Differential Equations[END_REF][START_REF] Kawashima | System of a Hyperbolic-Parabolic Composite Type, with applications to the equations of Magnetohydrodynamics[END_REF][START_REF] Serre | The structure of dissipative viscous system of conservation laws[END_REF], for any ε ∈ (0, ε], the system (4.22) has a unique classical solution w ε taking its value in the open set O 1

w ε -w ⋆ ∈ C 0 ([0, t ], H l ) ∩ C 1 ([0, t ], H l-2 ), w εii -w ⋆ ii ∈ L 2 ([0, t ], H l+1 ), (4.31)
w ε (x, t) ∈ O 1 , (x, t) ∈ R d ×[0, t ]. (4.32)
For each fixed ε ∈ (0, ε] we may thus define a maximal time tε by tε = sup{ t > 0 : w ε satisfies (4.31) and (4.32) }.

This maximum time tε depends on the relaxation parameter ε and may tend to zero as ε → 0. In order to establish that tε remains at finite distance form zero, we may use the existence of approximate solutions defined on a finite interval that satisfy the following properties (L 1 )-(L 3 ).

(L 1 ) There exists a finite time t * > 0 independent of ε and a family of approximate solutions

w ′a ε ∈ C 0 ([0, t * ], H l ) for ε ∈ (0, ε] such that w ′a ε -w ⋆ ∈ C 0 ([0, t * ], H l ) ∩ C 1 ([0, t * ], H l-2 ), (4.33) w ′a εii -w ⋆ ii ∈ L 2 ([0, t * ], H l+1 ), (4.34)
w ′a ε (x, t) ∈ O 0 , (x, t, ε) ∈ R d × [0, t * ]× ∈ (0, ε]. (4.35) (L 2 ) Letting t m = min(t * , tε ), we have (4.36) lim ε→0 sup t∈[0,tm) |w ε (t) -w ′a ε (t)| L ∞ = 0.
(L 3 ) There exists a constant c independent of ε ∈ (0, ε] such that

(4.37) sup t∈[0,tm) |w ε (t) -w ′a ε (t)| 2 l + ε tm 0 |∂ x w εii (t) -∂ x w ′a εii (t)| 2 l dt ≤ c.
A key tool in order to establish existence results by using approximate solutions with (L 1 )-(L 3 ) is the following convergence-stability lemma [START_REF] Yong | Basic aspects of hyperbolic relaxation system, in advances in the theory of shock waves[END_REF].

Lemma 4.8. Assume that there exists a family of approximated solutions w ′a ε for ε ∈ (0, ε] such that (4.33)-(4.37) hold. Let tε be the maximal existence time for a classical solution with H l regularity and values in O 1 . Then for ε small enough we have tε > t * .

Proof. Arguing by contradiction we assume that there exists a sequence (ε i ) i∈N with lim i→∞ ε i = 0 and tεi ≤ t * . From (4.37) and the properties of w ′a ε we obtain that sup t∈[0, tε i ) |w εi (t)w ⋆ | 2 l ≤ c and from the governing equations in normal form that ∂ t w εi i (t) ∈ L 2 (0, tεi ), H l-1 . Similarly, since

tε i 0 |∂ x w εiii (t)| 2
l dt ≤ c/ε i from (4.37) we also get that ∂ t w εi ii (t) ∈ L 2 (0, tεi ), H l-1 , so that finally ∂ t w εi (t) ∈ L 2 (0, tεi ), H l-1 . The existence of lim t→ tε i w εi (t) in H l-1 then follows from the estimate

|w εi (t ′ ) -w εi (t)| l-1 ≤ c √ t ′ -t valid for 0 ≤ t ≤ t ′ < tεi .
The limit lim t→ tε w ε (t) then exists in H l-1 and also as a continuous function since l -1 ≥ [d/2] + 1. Moreover, from (4.37) and standard functional analysis arguments, the limit w εi ( tεi ) also belongs to H l . Since the trajectories w ′a ε (t) for 0 ≤ t ≤ tε ≤ t * are included in the compact set O 0 ⊂ O 1 , we obtain from (4.36) that the values of w ε (x, t) are necessarily in the open set O 1 and remain at finite distance from its boundary for ε small. In particular, the limits w εi (x, tεi ) also belongs to the open set O 1 for ε i small enough. We may then use a local solution starting from w εi ( tεi ) showing that the solution can be extended beyond tεi and contradicting the definition of tεi .

Another method of proof consists in noting that the local existence time for parabolic systems with smooth uniformly bounded coefficients is uniform with respect to initial solutions that are a priori bounded in H l . Classical solutions with a fixed and small enough ε can thus be extended beyond tεi from the uniform H l a priori bound [START_REF] Yong | Singular Perturbations of First-Order Hyperbolic System[END_REF].

In order to establish (L 1 )-(L 3 ) and the properties (4.33)-(4.37), we may now use the approximate solutions w ′a ε constructed in Section 4.1 that satisfy (4.33) and (4.35) when w (O w 0 ) ⊂ O 0 . We further establish in the next section a stronger version of (4.36) and (4.37) by proving that (4.38) sup

t∈[0,tm) |w ε (t) -w ′a ε (t)| 2 l + ε tm 0 |∂ x w εii (t) -∂ x w ′a εii (t)| 2 l dt ≤ cε 4 , 0 < ε ≤ ε.
By combining these estimates with the convergence-stability lemma, we obtain the following existence and approximation theorem.

Proof of the error estimates

We establish in this section the error estimates (4.38) for w ε -w ′a ε used in Theorem 4.9. To this aim, we consider a solution w ε of the equations in normal form with the regularity class w ε -w ⋆ ∈ C [0, t ], H l ) and initial condition w ε . We also consider the approximate solution w ′a ε = w(w ′a ε ) that satisfies the approximate system (4.25) with initial condition (3.43). The solution w ε is defined over an interval [0, tε ) and the approximate solution over [0, t * ] and we have to compare these solutions over [0, t m ) where t m = min(t * , tε ). Note that by assumption we have w ε ∈ O 1 and w ′a ε ∈ O 0 ⊂ O 1 in such a way that both |w ε | L ∞ and |w ′a ε | L ∞ are already a priori bounded. Moreover, it is sufficient to establish the estimates (4.38) for smooth solutions since we may use mollifiers and convolution operators.

Step 0. Preliminaries. In the following δ ≤ 1 denotes a generic small constant and c ≥ 1 a generic large constant, both depending on the system coefficients and initial data but, independent of ε. The various occurrences of these constants may be distinguished and the minimum of all δ and the maxima of all c may be taken at the end of the proof so that only single constants ultimately remain. For k ≥ 0 and φ ∈ H k we introduce (4.40)

E 2 k (φ) = 0≤|α|≤k |α|! α! R d A 0 (w ε )∂ α φ, ∂ α φ dx,
we define the error e ε by e ε = w εw ′a ε and the error for the parabolic variable by e εii = w εiiw ′a εii . The error is such that e ε ∈ C [0, t m ), H l ∩ L 2 [0, t m ), H l+1 and ∂ t e ε ∈ C [0, t m ), H l-2 and we may form the difference between equations (4.22) and (4.25). This yields the following equation for e ε written in nonconservative form for convenience (4.41)

A 0 (w ε )∂ t e ε + j∈D A j (w ε )∂ j e ε -ε i,j∈D B ij (w ε )∂ i ∂ j e ε + 1 ε L(w ε )e ε = h,
where h = h cv + h ds + h sr + h pt and

h cv = - i∈D A 0 (w ε ) A -1 0 (w ε )A i (w ε ) -A -1 0 (w ′a ε )A i (w ′a ε ) ∂ j w ′a ε , h ds = i,j∈D εA 0 (w ε ) A -1 0 (w ε )B ij (w ε ) -A -1 0 (w ′a ε )B ij (w ′a ε ) ∂ i ∂ j w ′a ε + εA 0 (w ε ) A -1 0 (w ε ) d′ (w ε , ∂ x w ε ) -A -1 0 (w ′a ε ) d′ (w ′a ε , ∂ x w ′a ε ) , h sr = - 1 ε A 0 (w ε ) A -1 0 (w ε )L(w ε ) -A -1 0 (w ′a ε )L(w ′a ε ) w ′a ε , h pt = -A 0 (w ε )A -1 0 (w ′a ε ) ε 2 G ′ 2ε + ε 3 G ′ 3ε + εF ′ ε .
These residuals h cv , h ds , h sr , and h pt are respectively due to convection, diffusion, sources and perturbations and the quadratic residual d′ is defined by d′ (w,

∂ x w) = d(w, ∂ x w)+ i,j∈D ∂ w B ij (w)∂ i w∂ j w.
Applying ∂ α to equation (4.41) where |α| ≤ l we get that (4.42)

A 0 (w ε )∂ t ∂ α e ε + j∈D A j (w ε )∂ j ∂ α e ε -ε i,j∈D B ij (w ε )∂ i ∂ j ∂ α e ε + 1 ε L(w ε )∂ α e ε = h α + g α ,
where the residual h α is associated with the derivatives of h whereas g α is due to commutators:

h α = A 0 (w ε )∂ α A -1 0 (w ε )h , (4.43) g α = i,j∈D εA 0 (w ε )[A -1 0 (w ε )B ij (w ε ), ∂ α ]∂ i ∂ j e ε (4.44) - i∈D A 0 (w ε )[A -1 0 (w ε )A i (w ε ), ∂ α ]∂ i e ε - 1 ε A 0 (w ε )[A -1 0 (w ε )L ij (w ε ), ∂ α ]e ε .
Step 1. The lth order estimate. Multiplying equation (4.42) by ∂ α e ε and |α|!/α!, integrating over R d , summing over 0 ≤ |α| ≤ l, using the inequality |φ| L ∞ ≤ c 0 |φ| l-1 , the Gårding inequality

δ|∂ x φ ii | 2 0 ≤ i,j∈D R d B ij (w ε )∂ i φ, ∂ j φ dx + c 0 |∂ x w ε | 2 L ∞ |φ ii | 2 0 , keeping in mind that B ij (w ε )∂ i φ, ∂ j φ = B ii,ii ij (w ε )∂ i φ ii , ∂ j φ
ii where φ = (φ i , φ ii ) t , and integrating by part the second-order terms

- i,j∈D R d B ij (w ε )∂ i ∂ j ∂ α e ε , ∂ α e ε dx = i,j∈D R d ∂ w B ij (w ε )∂ i w ε ∂ j ∂ α e ε , ∂ α e ε dx + i,j∈D R d B ij (w ε )∂ j ∂ α e ε , ∂ i ∂ α e ε dx, it is obtained that ∂ t E 2 l (e ε ) + εδ|∂ x e εii | 2 l + δ ε |πe ε | 2 l ≤ c |∂ t A 0 (w ε )| l-1 + j∈D |∂ j A j (w ε )| l-1 |e ε | 2 l (4.45) + cε|∂ x w ε | 2 l-1 |e εii | 2 l + c 0≤|α|≤l R d h α + g α , ∂ α e ε dx .
The integration by part of second-order terms and the Gårding inequality both yield extra similar contributions in the form cε|∂

x w ε | 2 l-1 |e εii | 2 l ≤ cε|∂ x w ε | 2 l-1 |e ε | 2
l on the right hand side. The terms arising from |∂ j A j (w ε )| l-1 are majorized by c(1+|∂ x e ε | l-1 ) upon writing w ε = w ′a ε +e ε and using the a priori bound of w ′a ε in H l . For the time derivative |∂ t A 0 (w ε )| l-1 we may use the governing equation (4.22) for w ε and πw 0 = 0, to obtain that

|∂ t A 0 (w ε )| l-1 ≤ c j∈D |∂ j w ε | l-1 + j∈D ε|∂ j w εii | 2 l-1 (4.46) + i,j∈D ε|∂ i ∂ j w εii | l-1 + |π(w ε -w 0 )| l-1 ε .
Combining (4.45) and (4.46), using the identity w ε = w 0 + (w ′a εw 0 ) + e ε as well as |w ′a εw 0 | l-1 ≤ c ε + exp(-δτ ) , we obtain the lth order estimate

∂ t E 2 l (e ε ) + εδ|∂ x e εii | 2 l + δ ε |πe ε | 2 l ≤ c(1 + |e ε | l + ε|e ε | 2 l )|e ε | 2 l (4.47) + c ε exp(-δτ )|e ε | 2 l + c ε |πe ε | l-1 |e ε | 2 l + 0≤|α|≤l R d h α + g α , ∂ α e ε dx .
More specifically, the terms l /ε. We next have to investigate in the right hand side of (4.47) the terms associated with h α = A 0 ∂ α (A -1 0 h) where h = h cv + h ds + h sr + h pt as well as the terms associated with commutators g α .

|∂ j w ε | l-1 |e ε |
Step 2. Convective terms h cv . In order to estimate A 0 (w ε )∂ α A -1 0 (w ε )h cv , we first note that A 0 (w ε ) is bounded in L ∞ and that A -1 0 (w ε )h cv is in the form

A -1 0 (w ε )h cv = - i∈D A -1 0 (w ε )A i (w ε ) -A -1 0 (w ′a ε )A i (w ′a ε ) ∂ j w ′a ε .
Using the nonlinear estimate (4.2) with the variable w εw ′a ε , the L ∞ bounds of w εw ′a ε , and the estimates of

w ′a ε in in C [0, t * ], H l+1 yields 0≤|α|≤l ∂ α e ε , A 0 (w ε ) ∂ α A -1 0 (w ε ) h cv dx ≤ c|e ε | 2 l .
Step 3. Diffusion terms h ds . We first note that

A -1 0 (w ε )h ds = ε i,j∈D A -1 0 (w ε )B ij (w ε ) -A -1 0 (w ′a ε )B ij (w ′a ε ) ∂ i ∂ j w ′a ε + ε ij∈D A -1 0 (w ε )m ij (w ε ) -A -1 0 (w ′a ε )m ij (w ′a ε ) ∂ i w ′a ε ∂ j w ′a ε + ε ij∈D A -1 0 (w ε )m ij (w ε ) ∂ i e ε ∂ j e ε + ∂ i e ε ∂ j w ′a ε + ∂ i w ′a ε ∂ j e ε ,
where only parabolic components are involved for h ds , ∂ i w ′a ε , and ∂ i e ε , with the blocks B ii,ii ij , A ii,ii 0 , and m ii,ii,ii ij , due to the block structure of B ij , A 0 , and m ij . The norm A -1 0 (w ε )h ds l-1 is then directly majorized in the form

A -1 0 (w ε )h ds l-1 ≤ cε|e ε | l-1 + cε(1 + |e ε | 2 l-1 )(|e εii | l + |e εii | 2 l ),
using the L ∞ bounds of e ε as well as the H l+1 bounds of w ′a ε . This inequality may then be simplified into

A -1 0 (w ε )h ds l-1 ≤ cε(|e ε | l + |e ε | 4 l ), using |e εii | l ≤ |e ε | l and inequalities in the form x β ≤ c(x + x 4 ) valid for x ≥ 0 when 1 ≤ β ≤ 4.
We first consider the terms A 0 (w ε )∂ α A -1 0 (w ε )h ds such that |α| ≤ l -1. From the previous estimates, these terms can directly be estimated as

0≤|α|≤l-1 ∂ α e ε , A 0 (w ε ) ∂ α A -1 0 (w ε ) h ds dx ≤ cε(|e ε | 2 l + |e ε | 5 l ).
On the other hand, for the terms such that |α| = l, an integration by part in performed by writing α = e i0 + α ′ for some i 0 ∈ D where e 1 , . . . , e d denotes the canonical basis of R d and |α ′ | = l -1. We may then write

∂ α e ε , A 0 (w ε ) ∂ α A -1 0 (w ε ) h ds dx = ∂ i0 A 0 (w ε )∂ α e ε , ∂ α ′ A -1 0 (w ε ) h ds dx = ∂ w A 0 (w ε )∂ i0 w ε ∂ α e ε + A 0 (w ε )∂ i0 ∂ α e ε , ∂ α ′ A -1 0 (w ε ) h ds dx .
The term involving the higher order derivative ∂ i0 ∂ α e ε yields a majorizing term in the form

cε|∂ x e εii | l (|e ε | l + |e ε | 4 l ),
keeping in mind that only parabolic components are involved, whereas the remaining terms are easily majorized by cε(|e

ε | 2 l + |e ε | 6 l ). We have finally established that 0≤|α|≤l ∂ α e ε , A 0 (w ε ) ∂ α A -1 0 (w ε ) h ds dx ≤ cε|∂ x e εii | l (|e ε | l + |e ε | 4 l ) + cε(|e ε | 2 l + |e ε | 6 l ).
Step 4. Source terms h sr . We first note that A -1 0 (w ε )h sr may be rewritten in the form

A -1 0 (w ε )h sr = - 1 ε A -1 0 (w ε )L(w ε ) -A -1 0 (w ′a ε )L(w ′a ε ) (w ′a ε -w 0 ),
since L(w ε )w 0 = 0 and L(w ′a ε )w 0 = 0. We may then write

0≤|α|≤l ∂ α e ε , A 0 (w ε ) ∂ α A -1 0 (w ε ) h sr dx ≤ c 1 + 1 ε exp(-δτ ) |e ε | 2 l , by using |w ′a ε -w 0 | l ≤ c ε + exp(-δτ ) .
Step 5. Perturbation terms h pt . The terms arising from h pt may be split by considering separately

ε 2 G ′ 2ε , ε 3 G ′ 3ε , and εF ′ ε . The term associated with ε 2 G ′ 2ε is directly majorized in the form 0≤|α|≤l ∂ α e ε , A 0 (w ε ) ∂ α ε 2 A -1 0 (w ′a ε ) G ′ 2ε dx ≤ cε 2 |e ε | l ,
using the a priori bounds of G ′ 2ε in H l from Proposition 4.7. Similarly the term associated with εF ′ ε is majorized as

0≤|α|≤l ∂ α e ε , A 0 (w ε ) ∂ α εA -1 0 (w ′a ε ) F ′ ε dx ≤ cε exp(-δτ )|e ε | l ,
using the a priori bound of F ′ ε in H l from Proposition 4.7. We next consider the terms arising from G ′ 3ε such that |α| ≤ l -1. These terms may directly be estimated by using the a priori bound of G ′ 3ε in H l-1 and yields that 0≤|α|≤l-1

∂ α e ε , A 0 (w ε ) ∂ α ε 3 A -1 0 (w ′a ε ) G ′ 3ε dx ≤ cε 3 |e εii | l .
On the other hand for the terms arising from G ′ 3ε such that |α| = l, an integration by part in performed by writing α = e i0 + α ′ where i 0 ∈ D and |α ′ | = l -1. We may then write

∂ α e ε , A 0 (w ε ) ∂ α A -1 0 (w ε ) G ′ 3ε dx = ∂ i0 A 0 (w ε )∂ α e ε , ∂ α ′ A -1 0 (w ε ) G ′ 3ε dx = ∂ w A 0 (w ε )∂ i0 w ε ∂ α e ε + A 0 (w ε )∂ i0 ∂ α e ε , ∂ α ′ A -1 0 (w ε ) G ′ 3ε dx .
The higher order derivative contribution ∂ i0 ∂ α e ε yields a majorizing term in the form

|α|=l A 0 (w ε )∂ i0 ∂ α e ε , ∂ α ′ ε 3 A -1 0 (w ε ) G ′ 3ε dx ≤ cε 3 |∂ x e εii | l .
On the other hand, the remaining terms are directly majorized in the form

|α|=l ∂ w A 0 (w ε )∂ i0 w ε ∂ α e ε , ∂ α ′ ε 3 A -1 0 (w ε ) G ′ 3ε dx ≤ cε 3 (1 + |e ε | l )|e εii | l ,
and we have established that

0≤|α|≤l ∂ α e ε , A 0 (w ε ) ∂ α A -1 0 (w ε ) h pt dx ≤ cε 2 |e ε | l + cε exp(-δτ )|e ε | l + cε 3 (|e ε | l + |e ε | 2 l ) + cε 3 |∂ x e εii | l .
Step 6. Combined estimates. From the estimates associated with h cv , h ds , h sr , and h pt we deduce after some algebra, keeping only the dominant terms and using that ε ≤ 1, that 0≤|α|≤l

∂ α e ε , h α dx ≤cε|∂ x e εii | l (ε 2 + |e ε | l + |e ε | 4 l ) + cε(ε|e ε | l + |e ε | 6 l ) + c 1 + 1 ε exp(-δτ ) |e ε | 2 l + cε exp(-δτ )|e ε | l .
This estimate for the residuals associated with h α may then be combined with the lth order inequality to get that

∂ t E 2 l (e ε ) + εδ|∂ x e εii | 2 l + δ ε |πe ε | 2 l ≤ c(1 + |e ε | l + ε|e ε | 2 l )|e ε | 2 l (4.48) + c ε exp(-δτ )|e ε | 2 l + c ε |πe ε | l-1 |e ε | 2 l + 0≤|α|≤l R d g α , ∂ α e ε dx + cε|∂ x e εii | l (ε 2 + |e ε | l + |e ε | 4 l ) + cε(ε|e ε | l + |e ε | 6 l ) + c 1 + 1 ε exp(-δτ ) |e ε | 2 l + cε exp(-δτ )|e ε | l .
We may then only keep dominant terms in the right hand side, using ab ≤ 1 2 (a 2 + b 2 ) and ε ≤ 1. The term cε 2 |e ε | l yields in particular a majorizing factor in the form c|e ε | 2 l + cε 4 and the term ε exp(-δτ )|e ε | l yields the majorizing factors (c/ε) exp(-δt/ε)|e ε | 2 l + cε 3 exp(-δt/ε). In addition, we may absorb all ε 1/2 |∂ x e εii | l terms with the help of the parabolic regularity term ε|∂ x e εii | 2 l of the left hand side. This yields the extra terms cε|e ε | 2 l that is dominated by c|e ε | 2 l the term cε 5 that is dominated by cε 4 , as well as the term cε|e ε | 8 l that is majorized by c|e ε | 8 l . Using finally ε ≤ 1 and inequalities in the form x β ≤ c(x 2 + x 8 ) valid for x ≤ 0 when 2 ≤ β ≤ 8, it is established that

∂ t E 2 l (e ε ) + εδ|∂ x e εii | 2 l + δ ε |πe ε | 2 l ≤ c(|e ε | 2 l + |e ε | 8 l ) + c ε exp(-δτ )|e ε | 2 l (4.49) + c ε |πe ε | l-1 |e ε | 2 l + c ε 4 + ε 3 exp(-δτ ) + 0≤|α|≤l R d g α , ∂ α e ε dx .
It will be useful in the following to have an inequality for E 2 0 sufficient for our purposes. In this particular situation α = 0 there are no commutators so that h 0 = h and g 0 = 0. Proceeding as for the lth estimates it is first obtained that

∂ t E 2 0 (e ε ) + εδ|∂ x e εii | 2 0 + δ ε |πe ε | 2 0 ≤ c(1 + |e ε | l + ε|e ε | 2 l )|e ε | 2 0 (4.50) + c ε exp(-δτ )|e ε | 2 0 + cε|∂ x e εii | l |e ε | 2 0 + c ε |πe ε | l-1 |e ε | 2 0 + R d h, e ε dx .
Using indeed (4.45) and (4.46), the terms arising from

|∂ j w ε | l-1 |e ε | 2 0 are majorized by c(1+|e ε | l )|e ε | 2 0 , the terms ε|∂ j w εii | 2 l-1 |e ε | 2 0 are majorized by cε(1 + |e ε | 2 l )|e ε | 2 0 , the terms ε|∂ i ∂ j w εii | l-1 |e ε | 2 0 are ma- jorized by cε(1 + |∂ x e εii | l )|e ε | 2 0 .
We may use for the source terms the majorizing factors c |πe ε | l-1 + ε + exp(-δτ ) |e ε | 2 0 /ε. We next have to investigate the terms associated with h = h cv + h ds + h sr + h pt . From the expression of these residuals, the L ∞ bounds of w εw ′a ε , and the estimates of

w ′a ε in in C [0, t * ], H l+1 , it is directly obtained that e ε , h cv dx ≤ c|e ε | 2 0 , e ε , h ds dx ≤ cε|e ε | 2 0 + cε|e ε | l-1 |∂ x e εii | 2 0 + cε|e ε | 0 |∂ x e εii | 0 , e ε , h sr dx ≤ c 1 + 1 ε exp(-δτ ) |e ε | 2 0 .
Furthermore, for the terms arising from h pt , we split the contributions by considering separately

ε 2 G ′ 2ε , ε 3 G ′ 3ε
, and εF ′ ε and we obtain that

e ε , G ′ 2ε dx ≤ cε 2 |e ε | 0 ≤ c|e ε | 2 0 + cε 4 , e ε , F ′ ε dx ≤ cε exp(-δt/ε)|e ε | 0 ≤ c ε exp(-δt/ε)|e ε | 2 0 + cε 3 exp(-δt/ε), e ε , G ′ 3ε dx ≤ cε 3 |e εii | 0 .
Combining finally the previous estimates for α = 0 we have established that

∂ t E 2 0 (e ε ) + εδ|∂ x e εii | 2 0 + δ ε |πe ε | 2 0 ≤ c(1 + |e ε | l + ε|e ε | 2 l )|e ε | 2 1 + c ε exp(-δτ )|e ε | 2 0 (4.51) + cε|∂ x e εii | l |e ε | 2 0 + c ε |πe ε | l-1 |e ε | 2 0 + c ε 4 + ε 3 exp(-δτ ) .
Step

[φ(w ε ), ∂ α ]ψ 0 = 1≤|α|≤l [φ(w ε ), ∂ α ]ψ 0 ≤ c |φ(w ε ) -φ(w ⋆ )| l |ψ| l-1 .
Using the nonlinear estimate (4.2) and since both w ε and w ′a ε are a priori bounded in L ∞ , this inequality is simplified into 0≤|α|≤l [φ(w ε ), ∂ α ]ψ 0 ≤ c|w εw ⋆ | l |ψ| l-1 . Finally using w εw ⋆ = w εw ′a ε + w ′a εw ⋆ yields the useful estimate where we have used ψ 3/2 ≤ ψ 4 + ψ. Moreover ψ(0) = (1/ε 2 ) E 2 l w εw ′a ε (0) since w ε (0) = w ε and using | w εw ′a ε (0)| l ≤ cε 2 from (4.28) yields ψ(0) ≤ cε 2 . Assuming then that ψ(0) = cε 2 < 1, it remains true in an interval near zero in such a way that ψ 4 ≤ ψ on this interval. for the nonequilibrium system, then Π -1 e N ( B) is independent of u and ξ, so that a sufficient condition in order to obtain (N e ) is therefore that the nullspace of F = i∈D ∂ q f i ξ i contains Π -1 e N ( B). In this situation, we then have N ( B e ) = Π -1 e N ( B) that is independent of u and ξ and this situation notably arises for two-temperature fluids [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF] and multicomponent fluids [START_REF] Giovangigli | Multicomponent flow modeling[END_REF][START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF].

We also naturally assume that the initial value of u e is in H l and that u e (0) -u 0 (0) + εu 1 (0) l = u e (0) -u out ε (0) l = O(ε 2 ).

We then have a second-order difference estimate between the Hilbert and the Enskog expansions [START_REF] Yang | Validity of the Chapman-Enskog expansion for a class of hyperbolic relaxation system[END_REF].

Proposition 4.12. Let l ≥ l 0 + 1 and let t m > 0 such that both the second-order Hilbert expansion u 0 + εu 1 = u out Then there exists a constant c independent of ε such that u e (t) -u 0 (t) + εu 1 (t) l ≤ cε 2 .

Proof. This is a direct consequence of the estimate for the residual r out ε that is O(ε 2 ) and of the stability theorem presented in Appendix B of Reference [START_REF] Giovangigli | Asymptotic stability and relaxation for fast chemistry fluids[END_REF].

Thanks to the above result on u e = u ce 1 and u out ε = u 0 + εu 1 and Theorem 4.9, we obtain the following result about the initial-layer of Chapman-Enskog expansion. Proof. We may indeed write that u ε (t) -u e (t) + εu il 1 (t/ε) l ≤ u ε (t)-u 0 (t) + εu 1 (t) + εu il 1 (t/ε) l + u 0 (t) + εu 1 (t) -u e (t) l .

The estimate (4.60) then yields |u 0 (t) + εu 1 (t) -u e (t)| l ≤ cε 2 whereas we obtain from Theorem 4.9 that |u ε (•, t) -(u 0 (t) + εu 1 (t) + εu il 1 (x, t/ε))| l ≤ cε 2 and this completes the proof.

We have thus obtained the initial-layer for a two term Chapman-Enskog expansion, extending the results of [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF]. It is remarkable that, even though there exists a zeroth order initial-layer term q il 0 for the nonequilibrium system, the reduced equilibrium system is nevertheless a good approximated system from (4.61).

APPENDICES A Equivalence of equilibrium systems

Proposition A.1. The reduced equilibrium system (2.15) coincides with that obtained in the framework of entropic symmetrizable systems [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF].

Proof. The convective matrices have been obtained previously [START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF] in the form Π t e A i ∂ u u eq . Since and ∂ u v eq = (η vv ) -1 η vu at equilibrium, a calculation yields that Π t e A i ∂ u u eq = A e,e i -A e,r i (η vv ) -1 η vu which is easily shown to coincide with ∂ u F e u, v eq (u, 0) = ∂ u f i (u, 0) that is the convective matrix of (2.15).

The same calculation may be applied to the first part of the diffusion matrices and this yields the identity Π t e B ij ∂ u u eq = B e,e i -B e,r i (η vv ) -1 η vu so that these terms in both formulations also coincide. The remaining part of the diffusion matrices obtained previously is in the form Π e A i Λ A j ∂ u v eq = Π e A i Λ A j Π e ∂ 2 u η e where A i = A i ∂ u v = A i (∂ 2 u η) -1 and Λ is the generalized inverse of L = -∂ v Q with prescribed nullspace A 0 E and range ( A 0 E) ⊥ = A -1 0 E ⊥ where A 0 = ∂ v u = (∂ 2 u η) -1 . After some algebra, it is obtained that A i = ∂ w u A i ∂ v w and using ∂ v w = ∂ u w (∂ u v) -1 that ∂ v w = (η uu -η uv η -1 vv η vu ) -1 -(η uu -η uv η -1 vv η vu ) -1 η uv η -1 vv 0 I .

Using the special structure of E = R ne ×{0} it is then established that

Λ = 0 η uv η -1 vv 0 I 0 0 0 S -1 0 0 η -1
vv η vu I .

A direct calculation then yields that Π t e ∂ w u A i = [∂ u f i , ∂ q f i ] and ∂ v w Λ (∂ v w ) t simplifies into

∂ v w Λ (∂ v w ) t = 0 0 0 S -1 .

Finally

Π t e A i Λ A j Π e = Π t e ∂ w u A i ∂ v w Λ(∂ v w ) t (A j ) t (∂ w u) t Π e = ∂ q f i S -1 (∂ q f j ) t , and this completes the proof.

B The inner layer residual

Proposition 2 . 4 .(u 1 )

 241 Under Conditions (u 1 )-(u 4 ) the following properties hold. The Hessian matrix ∂ 2 u η e is positive definite over O u .

  ) and the expression (3.37) of G 1 established in Proposition 3.1. On the other hand, the slow projection of (3.40) yields the u 1 governing equation (3.11) already derived in Section 3.2.It is then natural to introduce the new expansion(3.41)

  We denote by | • | l the norm in the Sobolev space H l = H l (R d ) and otherwise | • | A in the functional space A. If α = (α 1 , . . . , α d ) ∈ N d is a multiindex, we denote as usual by ∂ α the differential operator ∂ α1 1 • • • ∂ α d d and by |α| its order |α| = α 1 + • • • + α d . For any scalar function φ the L 2 norm of k th derivatives is defined by (4.1)

Lemma 4 . 1 .

 41 We denote by l 0 the integer l 0 = [d/2] + 1. (i) Let k, l be integers with k ≤ l and l 0 ≤ l. Then for any f ∈ H k and g ∈ H l , we have f g ∈ H k and there exists a constant c 0 depending on l and d such that |f g| k ≤ c 0 |f| k |g| l .

12 )

 12 such that (3.7), (3.9)-(3.11), (3.26), (3.28), (3.30), (3.31), (3.5), and (3.6) are satisfied.

Proof.

  The regularity properties of w a ε and of the residuals G ε , G 1 , G 2ε , and F ε are direct consequences of the regularity properties established in Proposition 4.2, and of the expression of the residuals obtained in Proposition 3.1, Proposition 3.2 and Appendix B. The corresponding estimates are also consequences of the estimates established in Proposition 4.2.

(w 2 )(w 3 )

 23 and the source term Q = (∂ w u) t ∂ 2 u η Q have the following properties. (w 1 ) The matrix A 0 has a block structure A i,i 0 0 0 A ii,ii 0 and is symmetric positive definite over O w . The matrices A j , j ∈ D, are symmetric over O w . The matrix B ij has a block structure 0 0 0 B ii,ii ij and satisfy

2 l

 2 are majorized by c(1+|e ε | l )|e ε | 2 l , the terms ε|∂ j w εii | 2 l-1 |e ε | 2 l are majorized by (1+ε|e ε | 2 l )|e ε | 2 l , the terms ε|∂ i ∂ j w εii | l-1 |e ε | 2 l are partly compensated by the parabolic regularity term δε|∂ x e εii | 2 l of the left hand side of (4.45) and yield on the right hand side εc(1 + |e ε | 2 l )|e ε | 2 l , and the source term yields the majorizing quantity c |πe ε | l-1 + ε + exp(-δτ ) |e ε | 2

7 .

 7 Terms arising from commutators. The commutators terms in (4.44) are in the form [φ(w ε ), ∂ α ] where φ(w ε ) converges to φ(w ⋆ ) at infinity. Upon decomposing φ(w ε ) = φ(w ⋆ ) + φ(w ε )φ(w ⋆ ) and since [φ(w ⋆ ), ∂ α ] = 0 and [φ(w ε ), ∂ 0 ] = 0 we have the inequality 0≤|α|≤l

0≤|α|≤l[ 1 .( 1 + 2 l( 4 ( 1 +Step 8 . 2 l( 1 ++ δ ε |πe ε | 2 l ≤ c|e ε | 8 l + c ε |e ε | 3 l

 1124182123 φ(w ε ), ∂ α ]ψ 0 ≤ c(1 + |e ε | l )|ψ| l-Using this inequality with φ = A -1 0 A i and ψ = ∂ i e ε first yields that0≤|α|≤l i∈D ∂ α e ε , A 0 (w ε )[A -1 0 (w ε )A i (w ε ), ∂ α ]∂ i e ε dx ≤ c(1 + |e ε | l )|e ε | 2 l .Similarly with φ = A -1 0 B ij and ψ = ∂ i ∂ j e ε we have0≤|α|≤l i,j∈D ∂ α e ε , εA 0 (w ε )[A -1 0 (w ε )B ij (w ε ), ∂ α ]∂ i ∂ j e ε dx ≤ cε(1 + |e ε | l )|e εii | l |∂ x e εii | l .Finally, we also obtain for φ = A -1 0 L and ψ = e ε that 0≤|α|≤l∂ α e ε , 1 ε A 0 (w ε )[A -1 0 (w ε )L(w ε ), ∂ α ]e ε dx ≤ c ε |e ε | l )|πe ε | l-1 |πe ε | l .We have thus established that 0≤|α|≤l∂ α e ε ,g α dx ≤ c(1 + |e ε | l )|e ε | .52) + cε(1 + |e ε | l )|e εii | l |∂ x e εii | l + c ε |e ε | l )|πe ε | l-1 |πe ε | l .Final estimates. Combining (4.49) and (4.52) and using 2 ≤ l we obtain that∂ t E 2 l (e ε ) + εδ|∂ x e εii | 2 l + δ ε |πe ε | 2 l ≤ c(|e ε | 2 l + |e ε | 8 l ) + c ε exp(-δτ )|e ε | ε | l-1 |e ε | 2 l + cε 4 + cε 3 exp(-δτ ) + cε(1 + |e ε | l )|e εii | l |∂ x e εii | l + c ε |e ε | l )|πe ε | l-1 |πe ε | l .The ε 1/2 |∂ x e εii | l factor can then be absorbed using the term due to parabolic regularity ε |∂ x e εii | 2 l . This yields on the right hand side the quantity cε(1+ |e ε | 2 l )|e εii | 2 l that is dominated by c(|e ε | 2 l + |e ε | 8 l ). Moreover, the product (c/ε)|πe ε | l-1 |πe ε | l be eliminated from (4.53) by using the left hand side term (δ/ε)|πe ε | 2 l of the lth estimate, the left hand side term (δ/ε)|πe ε | 2 0 of the zeroth order estimate, multiplying the zeroth order estimate by a suitable large constant a, adding the result to the lth order estimate, and using the interpolation inequality |φ| l-1 ≤ c 0 |φ| θ l |φ| 1-θ 0 where θ = (l -1)/l. The extra term cε|∂ x e εii | l |e ε | 2 0 arising from the addition of the zeroth order estimate is then controled by using parabolic regularity and yields on the right hand side ε|e ε | 4 0 that is dominated by c(|e ε | 2 0 + |e ε | 8 0 ), whareas all other added terms from (4.51) are already present on the right hand side using 1 ≤ l. Letting E 2 l (e ε ) = E 2 l (e ε ) + aE 2 0 (e ε ), it has been established that ∂ t E 2 l (e ε ) + εδ|∂ x e εii | 2 l δτ ) |e ε | 2 l + c ε 4 + ε 3 exp(-δτ ) . 53) by ε 2 , and using |e ε | 2 l ≤ c E 2 l (w ε ) and ε ≤ 1, it is obtained that ∂ t ψ + δ ε |∂ x e εii | 2 l ≤ cψ 4 + c 1 + 1 ε exp(-δτ ) ψ + c ε 2 + ε exp(-δτ ) ,

ε and the equilibrium solution u ce 1 =

 1 u e exist over [0, t m ] with regularity(4.59) sup t∈[0,tm] |u e (t) -u ⋆ | 2 l+1 + ε tm 0 |u eii (t ′ ) -u ⋆ ii | 2 l+2 dt ′ ≤ c|u e (0) -u ⋆ | l+1 .

Proposition 4 . 13 .

 413 Keep the assumptions of Theorem 4.9, assume that the equilibrium solution also exists over [0, t * ] with regularity (4.59). Then there exist positive constants t * and c such that the systems (1.1) and (2.15) have unique solutions in C 0 ([0, t * ], H l ) for ε ∈ [0, ε] with

  |u ε (•, t) -(u e (•, t) + εu il 1 (•, t/ε))| l ≤ cε 2 .

  Π e = I 0 , A i = A e,e i A e,r i A r,e i A r,r i , ∂ u u eq = I ∂ u v eq ,

Proposition B. 1 . 1 0∂ 1 0∂ 2 w∂jw p 0 + j∈D 1 0∂w

 111201 The residual difference δR ε may be written δR ε = εF ε with F ε = F εcv + F εsr + F εds + F εqd where Fεcv = j∈D w Aj w0 + εw1 + α(w il 0 + εw il 1 ) dα (w il 0 + εw il 1 )∂j(w1 + R0) + j∈D w Aj w0 + εw1 + α(w il 0 + εw il 1 ) dα w il 1 ∂j w p 0 Aj w p 0 + αw il 0 + βε(w1 + αw il 1 + R0) (w1 + αw il 1 + R0) dα dβ w il 0 Aj w p 0 + w il 0 + αε(w1 + w il 1 + R0) dα (w1 + w il 1 + R0)∂jw il 0 + j∈D Aj w0 + εw1 + w il 0 + εw il 1 ∂jw il 1 ,

  Proposition 4.9. Keep the assumptions and notation of Proposition 4.2 and further assume that (w 1 )-(w 4 ), (N) and (C) hold. Assume that w(O w 0 ) ⊂ O 0 ⊂ O 0 ⊂ O w , and that w ε ∈ O 0 . Then there exist positive constants t * and c such that for ε sufficiently small, system (4.22) has an unique solution w ε over [0, t * ] with regularity class (4.31) and (4.32) such that Proof. The existence of w ε on a finite time interval is a consequence of the existence of the approximated solution w ′a ε , of Lemma 4.8, and of the error estimates (4.38) which is established in the next section.Corollary 4.10. The family of approximate solutions w a ε constructed in Section 4.1 also satisfies sup t∈[0,t * ] |w ε (t)w a ε (t)| l ≤ cε 2 . Using the diffeomorphism u → w and letting u a ε = u(w a ε ) and u ε = u(w ε ) we also have a second-order estimate in the formsup t∈[0,t * ] |u ε (t)u a ε (t)| l ≤ cε 2 .The estimate (4.39) is an improved version of a former O(ε 3/2 ) estimate and is established without extra regularity assumptions. The extra accuracy is obtained by first using w ′a ε instead of w a ε in the proof. Ultimately, since w ′a εw a ε = O(ε 2 ), it is also obtained that w ε (t)w a ε (t) = O(ε 2 ) but the proof requires to first use w ′a ε because of the entanglement of the asymptotic orders.

	(4.39)	sup t∈[0,t * ]	|w ε (t) -w ′a ε (t)| l ≤ cε 2 .

  /ε) ψ(t ′ ) dt ′ + cε 2 , /ε) dt ′ ≤ cε 2 exp c(t * + 1) , using t 0 exp(-δt ′ /ε) dt ′ = ε 1 -exp(-δt/ε) . Assuming finally that cε 2 exp c(t * + 1) ≤ 1,it is then guarantee that ψ remains lower than unity over [0, t m ) with t m ≤ t * , so that finally the estimate (4.56) is established over [0, t m ). Moreover, from ψ(t) ≤ cε 2 exp c(t * + 1) over [0, t m ), we further obtain that εδ tm 0 |∂ x e εii (t ′ )| 2 l dt ′ ≤ cε 4 and this completes the proof of (4.38).

							Then on this
	interval we get that
	(4.55) exp(-δt ′ and from Gronwall Lemma we obtain that ψ(t) ≤ c t 0 1 + 1 ε
	(4.56)	ψ(t) ≤ cε 2 exp c	0	t	1 +	1 ε	exp(-δt ′

Chapman-Enskog expansion

The Chapman-Enskog method allows to derive governing equations for reduced systems [START_REF] Chapman | The Mathematical Theory of Non-Uniform Gases[END_REF][START_REF] Ferziger | Mathematical Theory of Transport Processes in Gases[END_REF][START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Giovangigli | Multicomponent flow modeling[END_REF][START_REF] Kawashima | Dissipative structure and entropy for hyperbolic systems of conservation laws[END_REF][START_REF] Giovangigli | Volume viscosity and fast internal energy relaxation : Symmetrization and Chapman-Enskog Expansion[END_REF]. Considering in particular the system in quasinormal form (2.7), this asymptotic method yields an expansion in the form w ε = w ce 0 + εw ce 1 + O(ε 2 ) involving a zeroth order variable w ce 0 = (u ce , q ce 0 ) t with a slow component u ce that is a priori infinitely accurate u ε = u ce and higher order coefficients w ce i for i ≥ 1 in the form w ce i = (0, q ce i ) t . Its effective resulting accuracy then depends on the truncation order of the expansion. In other words, with Hilbert method, the slow variable u ε itself is expanded u ε = u 0 + εu 1 + O(ε 2 ) and each coefficient u i has its own system of equation whereas with the Chapman-Enskog method the slow variable equation is expanded and the corresponding solutions resulting after truncation, that may be denoted by u ce i , i ≥ 0, are ith order approximation of u ε = u ce . We first investigate in this section the difference between the zeroth and first order Hilbert and Chapman-Enskog expansions. We then discuss the inital layer associated with a two-term Chapman-Enskog expansion.

With the hyperbolic-parabolic system presented in Section 2, it is obtained after some algebra that at zeroth order the Hilbert and Chapman-Enskog expansions coincide with w ce 0 = w 0 = (u 0 , 0) t where u 0 = u ce 0 is the solution of the first order system (3.9). Such zeroth order approximations would typically yield Euler equations for a compressible gas. On the other hand, after some algebra, at first order of Enskog expansion, the coefficient w ce 1 is found in the form w ce 1 = (0, q 1 ) t where q 1 = q ce 1 is given by (3.10) and thus differ from the Hilbert first order term w 1 = (u 1 , q 1 ) t . Moreover the first order approximation u ce 1 of u ce is found to be the solution of the more accurate reduced system (2.14) that includes all O(ε) terms. This Chapman-Enskog first order approximation would typically yields the Navier-Stokes-Fourier system for a compressible gas. A first task in therefore to compare the first order Hilbert outer expansion (4.57)

to the solution u ce 1 = u e of (2.14). Since u e = u ce 1 is the solution of the reduced system (2.14), the natural method in order to compare u ce 1 = u e and u out ε = u 0 + εu 1 is to evaluate the outer residual r out ε defined by (4.58)

The difference between u ce 1 and u 0 + εu 1 will then be estimated by using a stability theorem for hyperbolic-parabolic systems.

Proof. This result from a direct calculation.

In order to compare the Chapman-Enskog first approximation u ce 1 = u e that is the solution of the reduced system (2.15) to the Hilbert approximate solution u 0 + εu 1 , we need some assumptions on the structure of the reduced equilibrium system. We assume in the following that Condition N e holds for the reduced system.

(N e ) The nullspace N ( B e (u, ξ)) of B e (u, ξ) = i,j∈D B e ij (∂ 2 u η e ) -1 ξ i ξ j does not depend on u ∈ O u and ξ ∈ Σ d-1 and B e ij (∂ 2 u η e ) -1 N ( B e ) = 0 for i, j ∈ D.

Note that from the expression (2.16) of the diffusion matrices B e ij (u) we have B e = Π t e B Π e + F S -1 F t where F = i∈D ξ i ∂ q f i in such a way that N ( B e ) = Π -1 e N ( B) ∩N (F ). With condition (N)

Proof. It is first established directly that the difference δR ε may be written

The full expression of δR ε = F ε is then obtained by using the governing equations of w il 0 and w il 1 , in order to eliminate the time derivatives ∂ t w il 0 and ∂ t w il 1 , and after lengthy algebra. The residual δR ε = εF 1ε may also be expanded in the form εF ε = εF 1 +ε 2 F 2ε where F 1 is directly associated with the second-order initial-layer corrector w il 2 . In particular, by adding both w ′ 2 and w il 2 contributions to the truncated expansion w a ε , one may cancel both first order residuals εG 1 and εF 1 . Proposition B.2. Denote by w il 2 = (u il 2 , q il 2 ) t the second-order initial-layer corrector. Then we have the ordinary differential equation Proof. The expression of F 1 is directly obtained by taking the limit ε → 0 in F ε . The ordinary differential equation for w il 2 is then obtained from the ε 1 inner expansion relation deduced from (2.7) after lengthy calculations.