A theory of generalised solutions for ideal gas mixtures with Maxwell-Stefan diffusion

After the pioneering work by Giovangigli on mathematics of multicomponent flows, several attempts were made to introduce global weak solutions for the PDEs describing the dynamics of fluid mixtures. While the incompressible case with constant density was enlighted well enough due to results by Chen and J\"ungel (isothermal case), or Marion and Temam, some open questions remain for the weak solution theory of gas mixtures with their corresponding equations of mixed parabolic-hyperbolic type. For instance, Mucha, Pokorny and Zatorska showed the possibility to stabilise the hyperbolic component by means of the Bresch-Desjardins technique and a regularisation of pressure preventing vacuum. The result by Dreyer, Druet, Gajewski and Guhlke avoids ex machina stabilisations, but the mathematical assumption that the Onsager matrix is uniformly positive on certain subspaces leads, in the dilute limit, to infinite diffusion velocities which are not compatible with the Maxwell-Stefan form of diffusion fluxes. In this paper, we prove the existence of global weak solutions for isothermal and ideal compressible mixtures with natural diffusion. The main new tool is an asymptotic condition imposed at low pressure on the binary Maxwell-Stefan diffusivities, which compensates possibly extreme behaviour of weak solutions in the rarefied regime.


Introduction: The Maxwell-Stefan Navier-Stokes equations for an isothermal ideal fluid mixture
The fluid mixture here under consideration consists of N ≥ 2 chemical species A 1 , . . . , A N subject to the following equations of mass and momentum balance ∂ t ρ i + div(ρ i v + J i ) = 0 for i = 1, . . . , N (1) The equations (1) are the partial mass balances for the partial mass densities ρ 1 , . . . , ρ N of the species. We shall use the abbreviation ̺ := fields b 1 , . . . , b N model the external body forces acting on each species per unit mass. The diffusion fluxes J 1 , . . . , J N are defined to be the non-convective part of the mass fluxes, in view of which they obey the necessary side-condition N i=1 J i = 0. One among the possibilities to close the model in thermodynamic consistent way, respecting also the constraint that the net diffusion flux vanishes, are the Maxwell-Stefan closure equations. In this approach, described a. o. by [Gio99], [BD15], [BDc], the diffusion fluxes J 1 , . . . , J N are defined implicitly as the solutions to a linear algebraic system with N × 3 rows. For i = 1, . . . , N the vector identities are required to be valid in R 3 . For 1 ≤ i < k ≤ N, the functions f i,k are given friction coefficients representing binary friction forces between the species; We define f k,i := f i,k whenever i > k (symmetric interaction). The friction coefficients f i,k are the inverse Maxwell-Stefan diffusivities − D ik , but we shall avoid this heavier notation. We moreover assume that these coefficients are strictly positive functions of the state variables. In (3), we write ∇µ k and b k instead of the correct ∇(µ k /θ) and b k /θ for the sake of visual convenience. This is justified as the temperature θ is overall assumed constant. We moreover denote y =ŷ(ρ) the vector of the mass fractions: y = (y 1 , . . . , y N ) with y i := ρ i / N j=1 ρ j =:ŷ i (ρ) for i = 1, . . . , N. The quantities d i := ρ i N k=1 (δ i,k − y k ) (∇µ k − b k ) on the right-hand side of (3) are called the generalised driving forces for diffusion. With the help of the matrix the equations (3) can also be rephrased in the more compact way where R = diag(ρ 1 , . . . , ρ N ), P (y) is the projector I − 1 N ⊗ y and 1 N = (1, . . . , 1) ∈ R N .
In order to close the model (1), (2) with (3), we must specify the mass based chemical potentials µ 1 , . . . , µ N in their relationship to the state variables and the pressure. In this paper we restrict our considerations to so-called ideal mixtures. In this case the chemical potentials obey The function g i (p) is given, it is the Gibbs-energy of constituent i at temperature θ. In other words, denotingρ i (θ, p) the mass density of the ith constituent as function of temperature and pressure, then the derivative g ′ i (p) is nothing else but 1/ρ i (θ, p). Concerning the second contribution in (6), the quantities x 1 , . . . , x N denote the number fractions. They are related to the mass densities via x i =:= ρ i /(m i N j=1 (ρ j /m j )) =:x i (ρ).
Here and throughout the paper, the positive constants m 1 , . . . , m N denote the elementary masses of the species.
In order to integrate the definition (6) into a thermodynamically consistent framework, we moreover have to provide the definition of µ as variational derivative of the Helmholtz free energy of the system. Here we shall rely on the assumption that the Helmholtz free energy possesses only a bulk contribution with density ̺ψ. Moreover this function possesses the special form ̺ψ = h(ρ 1 , . . . , ρ N ), where h : R N + → R is convex and sufficiently smooth on the entire range of admissible mass densities. Then µ 1 , . . . , µ N are related to the mass densities ρ 1 , . . . , ρ N via To make the link between, on the one hand, the free energy h as function of the mass densities and, on the other hand, the prescribed chemical potentials (6), we use the Euler equation that connects all state variables where for simplicity θ does not occur, as it is a constant. We solve (8) with the integration method proposed in the paper [BDa]. We first define a function of pressure and mass densities via V (p, ρ) := N i=1 g ′ i (p) ρ i . The dimensionless function V (p, ρ) acts as a 'volume potential' and it describes the contribution of volume extension to the free energy density. Then the resulting pressure state equation p =p(ρ) is implicitly obtained from the identity Under the assumption -to be made more precise below -that the functions p → g i (p) are sufficiently smooth, concave and non-decreasing, the implicit function ρ →p(ρ) is uniquely determined by (9) and sufficiently smooth. The Helmholtz free energy density h compatible with (6) and (7) assumes the form Up to the integration constants C(θ) and c 1 (θ), . . . , c N (θ), the function h is likewise uniquely determined by the definition (7).
As an illustration of these concepts and definitions, consider the case of all g i are logarithmic, more precisely in the form where g 0 ∈ R N are reference values,v 0 1 , . . . ,v 0 N > 0 are reference volumes per unit mass, and p 0 > 0 is a reference pressure. Then, inversion of (9) yields a state equationp(ρ) := p 0 N i=1v 0 i ρ i , and the associated form of the free energy is easily gained from (10).
Let us remark in passing that we call a mixture ideal if it is volume additive -see for instance [Brd81], section 4.1.2. If the chemical potentials are subject to (6), the specific mixture volume isv = N i=1 g ′ i (p) y i , and it is nothing else but the sum of the specific volumes of the constituents weighted by their relative contributions to the total mass. Reversely, it can be shown in rigorous mathematical framework (cf. [BDa]) that the assumption of volume additivity generates the splitting (6) of the chemical potentials, and that the only compatible form of the Helmholtz potential is given by (10).
The initial-boundary-value problem. The universal equations (1), (2) supplemented by the Maxwell-Stefan closure equations (3), in which µ and p are defined via (6) and (9), now form a closed system of partial differential equations for the densities ρ 1 , . . . , ρ N and the velocities v 1 , v 2 , v 3 . We call it the Maxwell-Stefan Navier-Stokes system for an isothermal ideal compressible fluid mixture.
In this paper, we shall investigate the resolvability of typical initial-boundary-value problems for this PDE-system in a cylindrical domain Q T := Ω×]0, T [ with a bounded, open, connected cross-section Ω ⊂ R 3 and T > 0 a finite time. We consider the initial conditions For simplicity, we consider only linear homogeneous boundary conditions on the lateral surface S T : Definition 1. The initial-boundary-value problem consisting of the Maxwell-Stefan Navier-Stokes system (1), (2), (3) on Q T with the initial conditions (12), (13) on {0} × Ω and (14), (15) on S T is denoted (P ).
Some notations. In this paper, vectors of mass densities ρ live in the state space R N + := {ρ ∈ R N : ρ i > 0 for i = 1, . . . , N} . Vectors of mass fractions y and number fractions x live on the hyper-surface Given a state ρ ∈ R N + , the associated mass and number fractions are given as We shall encounter several tensors of size N × 3 like for instance the diffusion fluxes J or the body forces b. For such rectangular matrices, the line index is the upward one, and the column index the downward one.
Scalar products are denoted · in R N and R 3 indifferently. For two m × n matrices A, B, we denote A : For a differentiable function f of the state variables, we denote ∂ ρ i f , or alternatively f ρ i the partial derivatives, and ∇ ρ f = (∂ ρ 1 f, . . . , ∂ ρ N f ) = f ρ is the gradient.
We denote 1 N = (1, . . . , 1) ∈ R N , and {1 N } ⊥ ⊂ R N is the orthogonal complement of the span of 1 N . The orthogonal projection In the identity (10) we shall choose C(θ) = 1 and c i (θ) = 0, which is of no consequence for the mathematical analysis.

The main result
2.1. State of the art. While the modelling of mass, momentum and energy transport in multicomponent reactive fluids is probably as old as thermodynamics, the mathematical investigation of well posedness issues concerning these models is rather recent. In the very comprehensive book [Gio99], Giovangigli was the first to develop a systematic view on modelling and analysis. In the chapters 8 and 9, the fundamental underlying hyperbolic-parabolic structure of multicomponent flows with Maxwell-Stefan diffusion is first exhibited. To reveal this structure he resorts to equivalent so-called entropic variables, a tool which had been introduced and developed before in the context of hyperbolic conservation laws: See [Gio99], Ch. 8 for details. In the case of smooth coefficients, these techniques allowed in [Gio99], Th. 9.4.1 to establish global-in-time existence results for strong solutions to the Cauchy problem, if the initial conditions are close enough to a constant equilibrium state. These breakthrough results in principle gave a positive answer to the theoretical questions of well-posedness for multicomponent flow models, but they do not exhaust the field. Indeed, they rely on several abstract assumptions concerning the transformation to entropic variables, and also on the assumption that there is a constant equilibrium. In particular, the weak solution analysis is entirely left open. Concerning the theory of strong solutions, several authors looked later at different aspects of the multicomponent flow models. Maxwell-Stefan diffusion in a much more reduced setting, without mechanics and energy balance, was for instance investigated in [HMPW17], while the incompressible Maxwell-Stefan-Navier-Stokes case with ̺ = const was studied in [BP17]. In both cases the hyperbolic component is avoided due to the assumptions of mechanical balance or incompressibility. In these papers, the PDE system is investigated in the original variables ρ 1 , . . . , ρ N , and local-in-time well-posedness is proved in classes of maximal parabolic regularity. In the investigation [BDb] on strong solutions, we investigate the well-posedness in classes of strong solutions (L p −setting) models similar to the one of the book [Gio99] (for the isothermal case, but more general free energy descriptions are included). By means of a similar passage to entropic variables in order to isolate the hyperbolic component, we prove the local-in-time well-posedness for initial-boundary-value problems. We do recover also the global existence result if the initial condition is sufficiently near to a stationary non-constant solution, but only for a finite time-horizon. Our estimates yield stability of the hyperbolic variable ρ only for smooth equilibrium solution and initial data. Another recent investigation worth mentioning on strong solutions is [PSZ18].
On the weak solution side, not only the mixed parabolic-hyperbolic character of the PDE system is the source of tough problems. The question of positivity and/or species vanishing is also a major issue. In the approach of strong or classical solutions, this problem is usually handled using the short-time vicinity to the initial data, the global vicinity of the equilibrium solution (cf. [Gio99], Ch. 9 and [BDb]) or using the maximum principle for smooth solutions (cf. [Gio99], 7.3.3, [HMPW17], [BDc], section 6). Moreover, relaxing the regularity of the velocity field in the Navier-Stokes equations implies that the solution to the continuity equation is not any more bounded away from zero and infinity: complete vacuum and blow-up of the total mass density are phenomena that we must expect. Thus, the main challenges for a theory of weak solutions are • The mixed parabolic-hyperbolic character; • The dilute limit (vanishing constituents); • Extreme behaviour of the total mass density (vacuum and blow-up).
To get rid of a part of these difficulties, the incompressibility assumption was used in [CJ15] to establish the global-in-time existence of some weak solutions. Interestingly, the passage to entropic variables was used again in this investigation, showing that this reformulation provides a natural framework both for the weak and the strong solution analysis. Very similar results are obtained in [MT15] published the same year. While the assumption of incompressibility or mechanical equilibrium oblige in [CJ15] and [MT15] to some readjustments of the thermodynamics, the full models exposed in [Gio99] were investigated from the viewpoint of weak solutions by [MPZ15]. In this investigation, the PDE system is handled with the help of several steps of stabilisation/regularisation: A pressure blow-up at density zero is admitted to exclude vacuum, while the hyperbolicity of the total mass density is tempered by means of the so-called Bresch-Desjardins device. A special structure of the viscous stress tensor yields some bounds on the density gradient, so that all partial mass densities, and with them the pressure, are weakly differentiable in space.
Afterwards, coming from the different context of electrochemistry, the team of authors of [DDGG17a], [DDGG17b], [DDGG17c] proposed a weak existence theory for equations describing the multicomponent flow of charged constituents. Accordingly, the underlying thermodynamic model relies on the structure (6) which is slightly more general than the models in [Gio99]. In this context it is more natural to admit that the pressure of the constituents is growing super linearly, so that it was possible to prove the global existence of certain weak solutions employing the Lions compactness method -thus, with almost the same methods as for the single component fluids. In [DDGG17a], we also employed a change of variables similar to [Gio99] to isolate the hyperbolic density component. In one point however, the analysis remains unsatisfactory. We assume there Fick-Onsager diffusion fluxes of the form J = −M(ρ) (∇µ − b) -which in general is equivalent to Maxwell-Stefan diffusion -but have to assume moreover that the Onsager operator is uniformly strictly positive on the orthogonal complement of 1 N . This provides a uniform control on the gradients of the entropic variables. A drawback of this assumption is that the Maxwell-Stefan diffusion case is not any more covered in the dilute limit (vanishing constituents), as infinite diffusion velocities must occur.
In order to obtain a result including Maxwell-Stefan diffusion, the author relaxed in the subsequent paper [Dru16] the necessary assumptions for the Onsager operator. The correct dilute behaviour is recovered in the range of finite, positive mass densities, and it degenerates only on the vacuum. Unfortunately, this result is very technical, it relies on the Fick-Onsager formulation and, due to the additional features coming from the electrochemistry model, the preprint does not convey a clear and direct message.
From this brief review of the literature, we conclude that there is no published result on global weak solutions to the equations of multicomponent Maxwell-Stefan diffusion in a compressible fluid, as far as it avoids both stabilisations and/or the assumption that the equivalent Fick-Onsager form is uniformly positive on {1 N } ⊥ . The aim of this paper is to provide such a result in the simplest possible setting where only Maxwell-Stefan diffusion and the isothermal multicomponent compressible Navier-Stokes equations are taken into account, with the ideal form (6) of the chemical potentials. Compared with the previous investigations, the existence of global solutions shall be enforced here by means of two types of modelling assumptions: (a) At high densities, the pressure of the single constituents is a sufficiently fast super linear function of their density; (b) Maxwell-Stefan frictions coefficients f i,k are regular and strictly positive at finite pressure, but they might exhibit at low pressure a singular dependence on the state (cp. the formula (18)).
2.2. The growth conditions and the main result. As temperature is a constant, we reinterpret the functions g i occurring in (6) as functions of pressure only. For i = 1, . . . , N the function g i : ]0, +∞[→ R is subject to the following assumptions (A): (A1) g i ∈ C 2 (]0, +∞[); (A2) g ′ i (p) > 0 and g ′′ i (p) < 0 for all p > 0; (A3) lim p→0 g ′ i (p) = +∞ and lim p→+∞ g ′ i (p) = 0; The conditions (A2) express for each constituent the requirement that: first its densitŷ ρ i (p) = 1/g ′ i (p) is positive (first condition), and second that its volume g ′ i (p) is a strictly decreasing function on pressure (second condition). The first of the conditions in (A3) means that the lower-pressure threshold corresponds to infinite volume, and the second one that infinite pressure leads to volume zero. Up to now the assumptions on g i are obvious. For the mathematical analysis, we shall moreover require specific asymptotic growth assumptions for small and large arguments: (A4) For small and moderate arguments, we assume that g i behaves likec i ln p, more precisely that there are a positive constants c 1 < c 2 and M 0 such that (A5) For large arguments, we assume that g i behaves likec i p 1 α i with α i > 1, more precisely that there are α 1 , . . . , α N ≥ β > 1 and M 1 > 0 such that Concerning the friction coefficients f i,k occurring in the Maxwell-Stefan equations (3) we assume that they are given as functions of the state variables expressed by the pressure p and the composition vector x (number fractions). The domain of definition of composition vectors is the positive unit-sphere of the one-norm S 1 + , and we denote S 1 + the closure in R N . For the matrix of coefficients {f i,k }, we make for all i < k the following natural modelling assumptions (B): In order to formulate the specific asymptotic conditions, we distinguish between a lowpressure range p < p 1 and a normal pressure range p > p 1 for some p 1 > 0. We assume that (B3) There are constants 0 < f 0 ≤ f 1 < +∞ such that for the singular functions σ i,k defined just here below in (16) Recalling that m 1 , . . . , m N are the elementary masses of the species and p 0 the reference pressure, the general form of the coefficients σ i,k in (B3) is given by the formula in which the positive factorχ(p, x) is defined implicitly as the root of the equation In order to interpret (16), we choose identical masses m 1 = m 2 = . . . = m N =: m. This allows to eliminate the inexplicit factorχ, and we obtain that Insert the special choice g i (p) =c i ln(p/p 0 ) to find that the functions define the asymptotic behaviour for p → 0.
The conditions (B3), (B4) can therefore be interpreted as follows. For finite pressures, or equivalently, for finite total mass density, the friction coefficients f i,k (the inverse Maxwell-Stefan diffusivities) are bounded and positive. For extreme behaviour of the density (rarefied regime), the friction coefficients are singular and depend on the composition and the reference values for mass m, pressure p 0 and partial volumesc/p 0 in a complex manner. In this paper, we do not attempt to justify the condition (B3) from physical considerations. Let us merely notice that the possibilities to measure diffusivities for multicomponent systems are generally limited, not to mention approaching extreme states like vacuum. As to molecular dynamical simulations, they for the most rely on isobaric models, and are for this reason more suited to capture the compositional dependence of the f i,k than the dependence on parameters like pressure or total density. We refer to [BDc] for details. Hence, keeping of course all cautions that are necessary in making such statements, we point at the fact that a mathematically motivated condition like (B3) might also serve as a guideline to meaningfully extrapolate known models for friction coefficients.
2.3. Main result. Assume that Ω ⊂ R 3 is bounded and Lipschitz. For T > 0, denote Q T = Ω×]0, T [ and S T = ∂Ω×]0, T [. We shall use standard parabolic Lebesgue spaces L p,q (Q T ) with 1 ≤ p, q ≤ +∞. To vary a little bit, we use the original notations of the monograph [LSU68], putting the space integration index first. Recall that -v. problem (P ) consists of (i) A vector of mass densities ρ ∈ L γ,∞ (Q T ; R N ) assuming nonnegative values almost everywhere in Q T . Here γ > 1 is some constant fixed by the data; , these fields are subject to the following integral relations In order to solve the Maxwell-Stefan equations (3), (5), we must make sense of the generalised driving forces. In particular, we must explain by possibly vanishing mass densities why we are allowed to build the spatial gradient of µ :=μ(p(ρ),x(ρ)). To this aim we introduce for large parameters k ∈ N and vectors ξ ∈ R N a cut-off operator [ξ] k := ([ξ 1 ] k , . . . , [ξ N ] k ) and for s ∈ R, we define [s] k := sign(s) min{|s|, k}. For a vector ρ subject to (i), we then require in addition to (i), (ii), (iii) the following additional technical condition: Obviously, the vector fieldμ(p(ρ),x(ρ)) is finite almost everywhere in the set Due to the condition (iv), the matrix ∇P {1 N } ⊥μ(p(ρ),x(ρ)) is therefore identical with an element of L 2 (Q T ; R N ×3 ) at almost every point of A + .
With B(ρ) defined via (4), with R := diag(ρ 1 , . . . , ρ N ) and P = I − 1 N ⊗ŷ(ρ), solving the Maxwell-Stefan equations in the weak sense shall consist of two pointwise identities: Remark 1. For a weak solution satisfying (22) and (23), the flux J remains non identified only on the pathological vacuum set. For positive total mass density -or equivalently positive pressures -J satisfies the Maxwell-Stefan equations (cf. (22)) and exhibits the expected behaviour in the dilute limit (cf. (23)).
Theorem 2. We assume that Ω ⊂ R 3 is bounded and of Lipschitz class. Let T > 0. For the data, we assume that (a) The functions g i are subject to conditions (A), and the numbers α 1 , . . . , α N occurring Then, the problem (P ) possesses a global weak solution (ρ, v, J). In particular, it is subject Remark 3. In analogy to the results known for single components Navier-Stokes equations ( [Lio98], [FNP01], a. o.) and for multicomponent systems ([DDGG17a]), we expect that a similar result holds for max i=1,...,N α i < 3. To verify this claim would however considerably increase the technicality.
Structure of the paper. Due to several preliminary investigations, the steps of regularising the PDE system and constructing sufficiently regular approximate solution are in principle clear enough. We focus in this paper on the new aspects in the estimates and the convergence proof. In the next section, the general strategy is first exposed: Change of variables, new estimates, compactness argument. The remainder of the paper is the technical part: We present separately in two separate sections the detailed proofs concerning the change of variables and the resulting new robust estimates. The final section of the paper discusses the convergence of a typical approximation scheme.

Main ideas
As explained in the previous section, the first main question is how to make sense of the Maxwell-Stefan equations (3) even for extreme states of the system, in particular dilute or rarefied ones. Both phenomena cannot be avoided for weak solutions. From the technical viewpoint of analysis, we introduce a new parametrisation of the variable ρ in the state space R N + , in such a way as to obtain weak parabolic estimates for auxiliary variables w 1 , . . . , w N independently of vanishing or blow-up of the total mass density. We show that a refined analysis of the dissipation due to diffusion, when combined with the growth assumptions (B3) for the friction coefficients in the low-pressure regime, shall provide the robust parabolic estimates we need. The variables w 1 , . . . , w N shall have the meaning of some normalisation of the mass densities to reference pressure.
Thanks to these estimates, we can rephrase the pressure as a function of the total mass density and the parabolic variable w via p =p(ρ) = P(̺, w). The function P is increasing in the first component, and bounded differentiable in the remaining parabolic components. This structure allows to prove that the total mass density associated with a weak solution is in a compact set of L 1 (Q T ). We now proceed to explaining these main ideas in more details.
3.1. A fundamental direction vector. Consider a solution ρ, v to the problem (P ). In the absence of sources, external fluxes and forces, the energy (in)equality for the system reads as follows: As is well known -this shall be recalled below in more details -inversion of the Maxwell- Hence (24) provides, among others, a control on the quantity ζ : , which represents the dissipation (or entropy production) of diffusion.
To further characterise this ζ, we use (7) and the chain rule compute the spatial gradient ∇µ = D 2 h(ρ) ∇ρ with the Hessian D 2 h of the free energy. We obtain that Due to the properties of M(ρ), the matrix D 2 h(ρ) M(ρ) D 2 h(ρ) is symmetric, positive semi-definite, and possesses rank N − 1. Its kernel is the vector u = u(ρ) solution to D 2 h(ρ) u = 1 N . We next want to determine u. In order to compute the Hessian, we exploit the particular representation (10), or directly the formula (6). If all g i (p) are strictly concave functions of class C 2 (]0, +∞[) -as required in assumptions (A) -the implicit equation (9) forp can be used to compute that Thus, setting k B θ = 1 for the sake of easier notation, we find for the Hessian of the free energy the expression Recall in this point thatn(ρ) = j (ρ j /m j ) is the total number density. As will be shown below, this matrix is explicitly invertible (see the inversion Lemma 6), so that we can now compute directly that the solution to D 2 h(ρ) u = 1 N is in which ̺ = j ρ j .

3.2.
Re-parametrisation of the state space. We next introduce an parametrisation of the state space R N + using alternative variables. The aim is a description ρ = X(s, w 1 , . . . , w N ), where w corresponds to mass densities normalised at some fixed reference pressure p 0 , and the parameter s stands forp(ρ) or the pressure. For fixed w, the curve s → X(s) shall be the characteristics following the vector u of (27).
To describe this procedure, we let S 0 := {ρ ∈ R N + : Due to the definition of the pressure state equation (9), the identity S 0 = {ρ ∈ R N + :p(ρ) = p 0 } is valid. For w ∈ S 0 , consider on an interval I =]p 0 − δ, p 0 + δ[ the system of ordinary differential equationsẊ Here the vector field u is given by (27), andẊ := d ds X. Assume that X is a solution to (28) such that X assumes values in R N + . Then we easily compute that In order to justify the two latter equations, recall that (8) guaranties that ∇ ρp (ρ) = D 2 h(ρ)ρ, while the definition of u implies that D 2 h(ρ)u(ρ) = 1 N . Thus, (29) implies thatp(X(s)) = s + c. For s = p 0 , we use that X(p 0 , w) = w. Since w ∈ S 0 , we havê p(X(s, w)) =p(w) = p 0 showing that c = 0 and thatp(X(s, w)) = s for all s > 0 and all w ∈ S 0 . Thus, the curves s → X(s, w), s ∈ I described by (28) are parametrised by the pressure.
Recall at this pointμ is a map defined on ]0, +∞[×S 1 + viaμ(p, x) = g(p) + 1 m ln x (cf. (6)). Due to (31), the dissipation ζ := N i,j=1 M i,j (ρ) ∇µ i · ∇µ j due to diffusion possesses the equivalent representation ζ = , t)) and a unique normalisation w = w(x, t) ∈ S 0 . Our idea -to be proved in details in the subsequent sections -is that the growth conditions (B3) for the friction coefficients of the Maxwell-Stefan system allow to prove for the matrix M(ρ) an inequality from below: where d 0 > 0 is a constant and y =ŷ(ρ) = ρ/( i ρ i ) are the mass fractions.
This means that the control for of the integral of ζ = M(ρ) ∇µ w · ∇µ w next implies an estimate in in L 1 (Q T ) for with the mass fractionsŷ i (w) := w i / j w j associated with the normalised state w 1 , . . . , w N . From the identity (33), we now proceed in two steps. Using first the convexity of | · | 2 , and the fact thatŷ(w) is a vector of fractions, we obtain that Next, we expand in (33) via Thus, combining with (33), we get Now, a brief calculation using the definitions of the mass fractionsŷ(w) and the number fractionsx(w) shows thatŷ i (w) = m ixi (w)/ j m jxj (w). Therefore, due to the fact that x attains values in S 1 + , we show that Thus, the inequality (35) yields This also implies a bound for the gradient of w. Recall indeed that w ·v 0 = 1, due to the fact that w ∈ S 0 . This allows to show the identity w i =x i (w)/(m i jx j (w) (v 0 j /m j )). For a i := x i (w)/m i , we thus get w i = a 2 i /(a 2 ·v 0 ). We clearly obtain that |∇w| 2 ≤ c 1 |∇a| 2 , where c 1 depends on the numbers min m, max m, minv 0 , maxv 0 . Since also w·1 N ≥ minv 0 we finally obtain that 3.4. Proving compactness. Imagine that we have approximated the problem (P ) with a sequence of problems (P σ ) easier to solve. For example, we discuss below the case of replacing the Onsager operator with a stabilisation M σ (ρ) := M(ρ) + σ I of full rank in order to obtain a parabolic problem. The question is now whether the vector of mass densities ρ σ solution to (P σ ) converges, at least for an appropriate subsequence, to a weak solution for (P ). In view of the nonlinearity of the system, it is clear that obtaining strong convergence is a necessary step to prove existence for the limit problem.
In a series of studies (in particular in [DDGG17c], [Dru19]) devoted to this type of problems, we have shown that the compactness problem is reduced to obtaining the strong convergence of the hyperbolic components, that is, the sequence of total mass densities ̺ σ . Indeed, once the latter property is known to hold, the weak parabolic estimates (37) allow to show the compactness of the entire vector ρ σ 1 , . . . , ρ σ N by means of adapted Aubin-Lions (J.-L.) techniques. Now, the compactness of the total mass densities is essentially obtained by adapting the method of Lions (P.-L.) for single component Navier-Stokes to our case. We show here in essence how to derive this adaptation from the estimate (37).
Assume first that the growth of the pressure ρ →p(ρ) is sufficiently strong. This means that there are positive c 0 , c 1 such thatp(ρ) ≥ c 0 |ρ| γ − c 1 with γ > 3/2. To simplify the technical discussion, we will in fact assume for the proof of our main result that γ > 3, which corresponds to choosing max i=1,...,N α i < 3/2 in (A5). In this case we elementarily obtain from the mathematical theory of Navier-Stokes equations an a priori estimate in L 1+1/γ (Q T ) for the pressure, and in L 1+γ (Q T ) for the density, so that the product p ̺ is at least in L 1 (Q T ).
Then, due to many investigations (work of Lions, Feireisl and, for the multicomponent case, also our study [DDGG17c] a.o.), the validity of the Lions argument for the compactness of the total mass density is reduced to showing lim σ→0 Qt where p is the weak limit of p σ in L 1+ 1 γ (Q T ) and ̺ the weak limit of ̺ σ in L 1+γ (Q T ). To prove the property (38), we essentially employ the monotonicity of the curves s → X(s, w) of (28). The strict monotonicity allows to introduce for r > 0 the coordinate transformation (s, w) → (r, w) via r = N j=1 X j (s, w), which defines an implicit function s = P(r, w). In particular, the inverse map r → P(r, w) is increasing on ]0, +∞[. Thus, for some sequence {σ n } n∈N , let us change variables and represent p σn (x, t) = s n (x, t) and ̺ σn (x, t) = X(s n (x, t), w n (x, t)) · 1 N . Thanks to the procedure just described, we find s n = P(̺ σn , w n ). Then, the Lions compactness argument is reduced to obtaining lim n→∞ Qt This property has been proved among others in [DDGG17c]. We recall it briefly, abbreviating for simplicity ̺ n = ̺ σn etc. We start from the inequality (P(̺ n , w n )−P(̺, w n )) (̺ n − ̺) ≥ 0 valid for arbitrary positive̺. Integrated against any smooth nonnegative testfunction φ, this yields For m ∈ N, we let̺ m be a smooth positive function approximating ̺ and such that we can find for all m ∈ N some positive constants a m , b m for which 0 < a m ≤ ̺ m ≤ b m < +∞. Suppose now that for m fixed, we can establish for all n the property Then, (40) can be combined with the uniform estimate (37) for ∇w n in L 2 . The chain rule for Sobolev functions yields a uniform bound in L 2 (Q T ; R 3 ) for the sequence {∇P(̺ m , w n )}. Since ∂ t ̺ n is uniformly bounded in L 2 (0, T ; (W 1,6γ/(5γ−6) (Ω)) * ) (see the uniform estimates in [DDGG17b]), a div-curl type argument then shows that Here β m is a weak limit of P(̺ m , w n ). We can show that the functions β m are bounded in L 1+1/γ (Q T ) independently of m. Indeed, owing to the pressure growth, there are some constants c 0 , c 1 such that c 0 (r γ − 1) ≤ P(r, w) ≤ c 1 (r γ + 1), and (39) is valid. Thus Letting m tend to ∞, we get (38). Therefore we can expect that, for typical approximation schemes, the total mass densities converges strongly in L 1 (Q T ). The remainder of the paper is devoted to providing the complete proofs of the ideas sketched in this section. At this level, a somewhat more technical presentation cannot be completely avoided.
For the fractionx j (X(s, w)), the latter yields a representation We sum up over j = 1, . . . , N and, sincex maps into S 1 + , obtain that With a j = a j (s, w) :=x j (w) exp(m j (g j (p 0 ) − g j (s))), the identity (44) is an equation of the form N j=1 a j χ m j = 1 for the variable χ := (x k (X(s, w))/x k (w)) 1/m k exp(−(g k (p 0 ) − g k (s))). A brief analysis recalled in Lemma 5 shows that this equation possesses a uniquely determined and strictly positive root χ. Moreover, we have χ =χ(a 1 , . . . , a n ) with a certain regular function of the entries of a. In particular, we can rely on the following estimates Plugging this into (43), we obtain the equation x j (X(s, w)) =a j (s, w)χ m j (a 1 , . . . , a n ) with a j (s, w) :=x j (w) exp(m j (g j (p 0 ) − g j (s))) .

(46)
With the help of the latter relation, we can now derive a full representation of the solution to (41). We make use of the identity N i=1 g ′ i (s) X i (s, w) = 1 valid by definition of the pressure (cf. (9)), and further of X i (s, w) =n(X(s, w)) m ixi (X(s, w)) (the definition of number fractions and number densities). Forn, we obtain the equivalent definition n(X(s, w)) = 1/ N i=1 g ′ i (s) m ixi (X(s, w)). Hence, invoking (46) yieldŝ n(X(s, w)) = 1 We combine (46) and (47) with X k (s, w) = m kn (X(s, w))x k (X(s, w)) to obtain a representation andχ(a) =: χ defined to be the unique solution to the algebraic equation N j=1 χ m j a j = 1. Now, we claim having found in (48) the global solution to (41), which can be verified briefly. Abbreviating X = X(s, w), we indeed first notice that (48) implies, for all j, k = 1, . . . , N, that from which it is readily verified that (e j − e k ) · (μ(s,x(X)) −μ(p 0 ,x(w))) = 0. Since k, j are arbitrary, this yields that the X of (48) satisfies η · ∇ ρ h(X(s, w)) = η · ∇ ρ h(w) for all η ∈ {1 N } ⊥ . By these means,Ẋ(s, w) is orthogonal to D 2 h(X(s, w)) η for all η ∈ {1 N } ⊥ , and we conclude thatẊ(s, w) must be a vector parallel to u(X(s, w)). Next, the representation (48) also implies that k g ′ k (s) X k (s, w) = 1 and therefore s =p(X(s, w)) by the implicit definition ofp. Computing the derivative, this allows to show that 1 = X(s, w) D 2 h(X(s, w)) X(s, w). As we already know thatẊ is parallel to u = (D 2 h) −1 1 N , we find thatẊ(s, w) = u(X(s, w))/(X(s, w) · 1 N ). Thus (48) provides the solution, and it can clearly be extended to a map defined on ]0, +∞[×S 0 with values in R N + . Sinceχ is differentiable, the formula (48) moreover shows that the map w → X(s, w) is differentiable on S 0 . We next prove the bijectivity of X, for which we first compute the inverse. The starting point is (43), which we rephrase aŝ We sum up, and we obtain the equation (X(s, w)) .
This is an equation of the form N j=1 b j φ m j = 1. It follows that φ =χ(b 1 , . . . , b N ) with b j = b j (X) given as b j (X) :=x j (X) exp(−m j (g j (p 0 ) − g j (p(X)))) and that x k (w) =x k (X) φ m k exp(−m k (g k (p 0 ) − g k (p(X)))) .
In order to obtain w as a function of X, we recall that w ∈ S 0 so that w ·v 0 = 1. Therefore, we can represent w k as a function ofx(w) via w k = m kxk (w)/( N ℓ=1v 0 ℓ m ℓxℓ (w)). Thus andχ(b) =: φ defined to be the unique solution to the algebraic equation N j=1 φ m j b j = 1. For ρ ∈ R N + , one therefore has ρ = X(s, w) if and only if s =p(ρ) and w = w(X) is given by (50) at X = ρ. This proves that X is bijective.
Next we prove the estimates. Ad (a). From (48), we directly deduce that the bound (a) is valid.

(51)
Owing to N j=1 (w j /m j ) ≥ w · 1 N /|m| ∞ ≥ 1/(|m| ∞ |v 0 | ∞ ), and to the definition of χ, it follows that a(s, w)) , and now (45) and the fact that |a| 1 ≥ min ℓ exp(m ℓ (g ℓ (p 0 ) − g ℓ (s))) yield Similarly, we obtain the lower bound We next estimate the derivatives. Ad (c). Recall (27) to see that Since N j=1 g ′ j (p(ρ)) ρ j = 1, we obtain the bound With the help of (a), we therefore verify for X satisfying (41) that Ad (d). Due to the formula (48), the map w → X(s, w) is differentiable on S 0 . To more easily compute the derivatives, we start again from the identity (∇ ρ h(X(s, w)) − ∇ ρ h(w)) · η = 0 for all η ∈ {1 N } ⊥ (see (30)). Consider τ i := e i − ν i ν a tangent on S 0 , and differentiate in the direction of τ i the latter identity. Then, with δ w i := τ i · ∇ w denoting the tangential derivatives on S 0 , we obtain for all η ∈ {1 N } ⊥ that η · (D 2 h(X(s, w)) δ w i X(s, w) − D 2 h(w)τ i ) = 0. Therefore, there must exist some real number r i such that We multiply from left with X(s, w). Recall that D 2 h(ρ) ρ = ∇ ρp (ρ). Since The inversion formula for D 2 h in Lemma 6, and use of δ w i X(s, w) ·p ρ (X(s, w)) = 0 yield the result Here we also made use of the formula (25) for the derivative of the pressure. In order to obtain estimates, we first use (26). Together with the fact that w ∈ S 0 , it implies that Therefore, for any k, ℓ, the expressions |D 2 k,ℓ h(w) X ℓ (s, w)| are bounded by some uniform constant C times the quotient F k (s, w) := X k (s, w)/w k . For the functions q introduced in (54), we obtain that Now, by means of (55), we find that and, recalling (a), it follows that | δ w X(s, w)| ≤C sup k F k (s, w) (1+max g ′ (s)/ min g ′ (s)). We invoke (52) to see that | δ w X(s, w)| is bounded by a continuous function depending only on s.
We next provide the two auxiliary lemmas used in the preceding proof.
Lemma 6. Let ρ ∈ R N + , and D 2 h(ρ) is given by (26). Then with Λ : In these formula, V (p, ρ) := i g ′ i (p) ρ i and its derivatives V ρ , V p are evaluated at (p(ρ), ρ). Proof. Assume that v, w ∈ R N are two vectors subject to D 2 h(ρ) v = w. We multiply from left with ρ. Then, using that Here V ρ = g ′ (p(ρ)) and V p = j g ′′ (p(ρ))ρ j . We introduce the positively homogeneous function With the volume V of (9), we alternatively have We multiply with m i ρ i V ρ i and sum up over i = 1, . . . , N to find that Since we know that V ρ · v = −V p w · ρ, we can eliminate the quantity 1 The claim follows.
In order to complete the picture about the change of variables, it remains to introduce the equivalent representation of the pressure needed to prove the compactness (see the Section 3.4). In particular, we must verify (40).
Lemma 7. There is a function P ∈ C 1 (]0, +∞[×S 0 ) such that for ρ ∈ R N + , s > 0 and w ∈ S 0 connected via ρ = X(s, w), the identities s =p(ρ) = P( i ρ i , w) are valid. The function ̺ → P(̺, w) is strictly increasing. The norm of the tangential derivatives δ w P(̺, w) is bounded above by a continuous function of ̺ only.
We define P(̺, w) to be the implicit solution to ̺ = X(P(̺, w), w) · 1 N . Then, the derivatives obey We invoke (59) to see that ∂ ̺ P > 0. Due to (59) and Lemma 4, (d), we first obtain that | δ w P(̺, w)| is bounded above by a function of s =p(ρ) = P. Then, we invoke Lemma 4, relation (a), and the fact that the monotonously decreasing functions g ′ are invertible, deducing that Thus, | δ w P(̺, w)| is bounded above by a continuous function of X · 1 N = ̺.

Analysis of the Maxwell-Stefan equations and the dissipation bound
For this section we assume that the friction coefficients obey the assumptions (B). We consider the Maxwell-Stefan algebraic system We define a matrix B(ρ) = {b i,k (ρ)} i,k=1,...,N (cp. (4)) via so as to reformulate (60) as Here we recall that R := diag(ρ 1 , . . . , ρ N ) and that P (y) := I − 1 N ⊗ y. In all these equations we have used the natural transformation y =ŷ(ρ) for the mass fractions. In (60) and (61), we have by assumption f i,k = f i,k (p(ρ),x(ρ)).
It is well known and it can be easily checked that the matrix B(ρ) is a singular M -matrix with left kernel {1 N } and right kernel {y} (see among others [Gio99], par. 7.7, [BDc]). The use of generalised inverse to solve the Maxwell-Stefan equations is an original product by the first reference. Nevertheless we quote for convenience from the latter paper, recalling an invertibility result for (62).
We now prove the second essential result for this paper in order to obtain the robust estimates, and introduce first the fundamental function to formulate our asymptotic conditions on the friction coefficients in the Maxwell-Stefan system. For s ∈]0, +∞[, w ∈ S 0 we recall the definitions of quotients F i (s, w) = X i (s, w)/w i (cf. Lemma 4 and (51)). Define ρ := X(s, w). Making use of the identity (50) for w, we obtain for F i an equivalent form as function of the pressure s =p(ρ) and the number densities x =x(ρ) via Making use of (63), (64) and (16), we can verify that the singular functions σ i,k (s, x) occurring in the condition (B3) are expressed via Making use of (64), we have at ρ = X(s, w) the equivalent expression (X(s, w)) .

(66)
We now prove the robust estimates.
We next consider ρ = X(s, w) where X is the solution map for (41), s =p(ρ) and w ∈ S 0 . We denote W = diag(w 1 , . . . , w N ), and consider a matrix K := W 1 2 R −1 B α W 1 2 which is positive definite, symmetric, and possesses the entries (51). We now make use of the structural assumption (B3). For s ∈]0, p 1 [ and x ∈ S 1 + , we might , with f reg uniformly bounded from below and above by positive constants f 0 < f 1 , and the σ i,j (s, x) are by definition (cp. (65), (66)) (X(s, w)) .
For i = j, we use F i √ w i = √ F i ρ i and F ·ŷ ̺ = F · ρ, and we therefore obtain that Hence |K i,j | ≤ f 1 + α |w|∞ ̺ . For j = i, we use (B3) again and we have We therefore can bound max i,j |K i,j | ≤ f 1 + α |w| ∞ /̺. Thus the spectral radius r(K) is bounded by the same quantity, and In other words, the matrix K −1 − d(α) I is positive semi-definite, which means that We multiply from the left with the matrix P T W 1 2 and from the right with W 1 2 P , preserving the inequality, and thus P T B −1 α R P ≥ d(α) P T W P . Recall that B −1 α = B D + 1 α y ⊗ 1 N and P T y = 0, to see that P T B D R P ≥ d(α) P T W P . Finally, B D R P = B D R and Thus, we can let α → 0, obtaining the claim (1). We next prove (2). Due to the fact that B(ρ)ŷ(ρ) = 0, we can restate (62) in the form B(ρ)J = −d withJ := J − J·F y·F y. Notice thatJ is orthogonal to the vector F = F (s, w) related to ρ via (51) whenever ρ = X(s, w).
With τ := B(ρ) R, we then have For fixed ℓ ∈ {1, 2, 3}, consider the row forJ ℓ ∈ R N in the latter system. We multiply in (68) from the right with R −1J ℓ . Recalling that P (y) T y = 0, note that Thus, invoking (67) and (B3), The numbers F k y k F ·y are positive and they sum up to one. Using convexity and the fact that F ·J = 0 by construction, we thus obtain that Due to Hölder's inequality, |J ℓ | 1 ≤ |w| Using the uniform bound for w, we get It remains to recall thatJ = J − J·F y·F y. Since J · 1 N = 0, we must have J =J − yJ · 1 N , so that J satisfies the bound claimed in (2).
Lemma 9 provides the needed estimates at small pressure. Of course, for normal pressure, we have also the expected behaviour, as stated in the next remark.
Remark 10. For ρ = X(s, w) with s =p(ρ) > p 1 and w ∈ S 0 , the assumption (B4) for the Maxwell-Stefan coefficients yield min i =k f i,k (ρ) ≥ f 2 > 0 and max i =k f i,k ≤ f 3 . Under these conditions we obtain instead of the estimates (1), (2) of Lemma 9 the normal Maxwell-Stefan bounds, and we deduce that )). Proof. In order to prove (3), we replace the matrix K in the proof of Lemma 9 with K := R 1 2 R −1 B α R 1 2 , and show for the spectral radius that r(K) ≤ f 3 + α/̺. Arguing analogously, we obtain first that and letting α → 0, it follows that B D R ≥ 1 f 3 P T R P . Recall that ρ i = F i (s, w) w i . Thus R = diag (F 1 (s, w), . . . , F N (s, w)) W . Since F i (s, w) ≥ F i (s) according to Lemma 4, (b), we obtain for the range s =p(ρ) ≥ p 1 that R ≥ inf i=1,...,N, s≥p 1 F i (s) W , proving (3).

Proof of the main result
Due to several prior investigations, it is not necessary to proof the result of Th. 2 in every detail. We shall present the general sketch putting some additional emphasis only on the main steps and the new ideas. Approximation. We construct approximate solutions using entropic variables and a stabilisation of the diffusion fluxes. For σ > 0, consider a regularised kinetic matrix where B D (ρ) R, R := diag(ρ) is the matrix resulting from inversion of the Maxwell-Stefan matrix (cf. Lemma 8).
We further denote h the free energy function (10). As shown in [DDGG17a], this function is of Legendre type on R N + and the convex conjugate h * is of Legendre type on R N . Our approximation scheme consists in solving the system for variables µ σ = (µ σ 1 , . . . , µ σ N ) and subject to the identities We consider the initial conditions (13) for v, and for µ µ σ i (x, 0) = ∂ ρ i h(ρ 0 (x)) for x ∈ Ω, i = 1, . . . , N , with ρ 0 from (12). On S T we assume (14) and ν ·  σ = 0.
Due to the strict convexity of h * , the approximation scheme (71), (72) constitutes a nonlinear parabolic system for (µ, v). There are no algebraic constraints on the unknowns. Existence for the regularised problem can be proved as in [DDGG16], [DDGG17b] by means of a Galerkin approximation scheme. The weak solution possesses the regularity and it satisfies for all t ≤ T the energy dissipation inequality After estimating the right-hand side by means of Hölder's and Young's inequalities, we attain the inequality Recall that M σ = B D R + σ I with B D R resulting from inversion of the Maxwell-Stefan equations. The Remark 10, (4) shows that |B D (ρ) R| ≤ C |ρ| ifp(ρ) is sufficiently large. Thus, we can obtain the linear-growth estimate |M σ (ρ)| ≤ |B D (ρ) R| + σ ≤ C (1 + |ρ|). Showing next that h growth at least linearly, we can apply the Gronwall Lemma to (77), and we see that the left-hand is uniformly bounded.
Notice that, while estimating ∇µ σ in L 2 (Q T ) is a straightforward consequence of (77) and of the choice of M σ ≥ σ I, the control on µ σ L 2 (Q T ) results from more subtle arguments. Indeed, all contributions of µ σ to the dissipation are due to the gradient. To obtain a bound, we must in particular exploit the fact that Ω ρ σ i (x, t) dx is a conserved quantity for each i. We refer to [DDGG17b] for details.
From (76) and Lemma 11, we can motivate the regularity ρ σ ∈ L γ,∞ (Q T ; R N ) and √ ̺ σ v σ ∈ L 2,∞ (Q T ; R 3 ). Using (73), the Gibbs-Duhem equation . With analogous arguments:  σ ∈ L 2γ/(2+γ),2 (Q T ). Uniform estimates. Several uniform estimates are standard consequences of the energy inequality (76), so we will state them without comment. We again refer to [DDGG17b] for details. Notice moreover that among these estimates, several are in principle the same as for single component Navier-Stokes equations: Restricting for simplicity to the case γ > 3 as announced, we moreover can uniformly estimate the approximate pressures Only the uniform estimate for  σ is missing.
where we abbreviated R σ := diag(ρ σ 1 , . . . , ρ σ N ). By the definition of B D (see Lemma 8), the first contribution J σ solves the Maxwell-Stefan equations To obtain estimates, we next decompose the domain Q T between points for which p σ (x, t) ≤ p 1 and for which p σ (x, t) > p 1 , where p 1 is the constant occurring in the condition (B3).
If p σ (x, t) ≥ p 1 , we obtain the standard Maxwell-Stefan bound (see Remark 10), which means that |J σ (x, t)| 2 ≤ C |ρ σ | (−J σ : (∇µ σ − b)). We combine both uniform bounds, and they yield Passing to variables ρ σ = X(p σ , w σ ) as described in Lemma 4, we obtain via Lemma 9 and Remark 10 the additional bound This was demonstrated in the section 3. It is readily seen that the second contribution in Limit passage. Now we possess enough informations to extract weak limits (ρ, v, p, J) for . These are well known structural consequences for the Navier-Stokes operator (cp. [Lio98]).
Passing to variables ρ σ = X(p σ , w σ ), we have p σ = P(̺ σ , w σ ), and the Lions method provides the strong convergence of the mass densities ̺ σ → ̺ in L 1 (Q T ) -as demonstrated in section 3.4.
Once this point is clear, we can rephrase ρ σ = R(̺ σ , w σ ) with a certain mapping R, and we obtain the compactness of the whole vector of mass densities as for instance in [DDGG17c]. Due to the bijectivity of the variable change, the auxiliary variables w σ also converge strongly, for instance in L 1 (Q T ; R N ). These properties suffice for passing to the limit in the integral relations resulting from the distributional formulations of (71) and (72).
We next establish for our limit element ρ the property (iv) of weak solutions. For µ = µ(p(ρ),x(ρ)) = ∇ ρ h(ρ), the trick (30) implies again that P {1 N } ⊥ µ = P {1 N } ⊥ ∇ ρ h(w). Now, the vector ∇ ρ h(w) is, up to inessential constants, equal to lnx 1 m (w), and there is always at least one among the fractionsx 1 (w), . . . ,x N (w) being larger than N −1 . Thus, boundedness of η · lnx 1 m (w) for all η ∈ {1 N } ⊥ implies the boundedness of all entries of w from below. With this argument, we can verify that for all k > 1, the truncations [P {1 N } ⊥ ∇ ρ h(w)] k are identical with P {1 N } ⊥ ∇ ρ h([w] ℓ(k) ), where [·] ℓ(k) is some cut-off operator at small positive level. Such truncations are, with w itself, elements of W 1,0 2 (Q T ; R N ). Moreover, the map r → ∇ ρ h(r) is Lipschitz continuous in the subset of R N + such that r i ≥ ℓ(k) > 0 for all i. Thus, we can verify [P {1 N } ⊥ µ] k ∈ W 1,0 2 (Q T ; R N ) for all k, proving the property (iv). This allows to transform the relation (90). If (x, t) ∈ A + , then we find k large enough for which µ(x, t) = [µ(x, t)] k and P {1 N } ⊥ µ = [P {1 N } ⊥ µ] k . We therefore can reinterpret the term N j=1 (δ i,j −ŷ j (ρ)) e j D 2 h(w)∇w in the right-hand of (90) and this yields This proves the first condition (22) for a weak solution of Maxwell-Stefan.
Finally we want to prove, as required in the second condition (23), that J i = 0 if ρ i = 0 and p(x, t) > 0. Here we employ the intermediate result (69) in the proof of Lemma 9. It shows for our approximations, that the quantities 3 are bounded above by the functions ζ σ := −(∇µ σ − b) : J σ . We multiply with the indicator functions φ σ,ǫ . The sequenceJ σ φ σ,ǫ is again bounded in L 2γ 1+γ ,2 (Q T ). For its weak limitJ, we obtain in the set S 0 ǫ the identitỹ As (J σ,i ) 2 is bounded above by ζ σ times w σ i , we find thatJ = 0 for almost all (x, t) in S 0 ǫ such that w i = 0. But in the set S 0 ǫ , we know that w i = 0 is equivalent to ρ i = 0. Hence, we can infer from (92) that 0 = J i − J · F (p, w) F (p, w) ·ŷ(ρ)ŷ i (ρ) = J i for almost all (x, t) such that ǫ ≤ p(x, t) ≤ ǫ −1 and ρ i (x, t) = 0. This completes the proof that the weak solution (ρ, v, J) also complies with the second condition (23).