Well-posedness analysis of multicomponent incompressible flow models

In this paper, we extend our study of mass transport in multicomponent isothermal fluids to the incompressible case. For a mixture, incompressibility is defined as the independence of average volume on pressure, and a weighted sum of the partial mass densities stays constant. In this type of models, the velocity field in the Navier–Stokes equations is not solenoidal and, due to different specific volumes of the species, the pressure remains connected to the densities by algebraic formula. By means of a change of variables in the transport problem, we equivalently reformulate the PDE system as to eliminate positivity and incompressibility constraints affecting the density, and prove two type of results: the local-in-time well-posedness in classes of strong solutions, and the global-in-time existence of solutions for initial data sufficiently close to a smooth equilibrium solution.


Multicomponent diffusion in an incompressible fluid
In this paper we study the well-posedness analysis in classes of strong solutions of class-one models 1 of mass transport in isothermal, incompressible multicomponent fluids. This investigation is a direct continuation of results obtained recently concerning the compressible case in [BDb], and the weak solvability of the incompressible model in [Dru19]. Performing the incompressible limit (the low-Mach number limit) in models for fluid mixtures and for multicomponent fluids is desirable both from the practical and the theoretical viewpoint. On the one hand, fluid mixtures occurring in applications are often incompressible, and the limit passage reduces the stiffness of the models by eliminating the parameter which is practically infinite. On the other hand, the low-Mach number limit leads to a type of incompressibility condition which has not yet been studied in the context of mathematical analysis for fluid dynamical equations. We are interested in the second type of issue, that is, the theoretical issues of unique solvability and continuous dependence in classes of strong solutions for the underlying PDEs.
The model class for multicomponent transport in fluids here under study is the one proposed in [BD15], also applied to mixtures with charged constituents in [DGM13], [DGM18]. Concerning the fundamentals of thermodynamics for fluid mixtures, the reader is referred to these papers, or to the book [Gio99]. The model for Mach-number zero (incompressibility constraint) is based on I. Müller's definition of incompressibility as invariance of the volume under pressure variations [Mü85], [GMR12]. More directly, we follow the recent example of [DGM13] (formal limit), and the more general road map proposed in the Section 16 of [BD15]. In [BDa] we propose a derivation of the incompressible limit starting from a few postulates of mathematical nature about the structure of the Helmholtz free energy. Similar concepts have been exposed and discussed in a few research papers like [Mil66,JHH96,DNB + 15]. Incompressible mixtures are also conceptualised in the book [PS14]. The corner stone of these works is that incompressibility for a multicomponent system means the invariance of average volume under pressure variations. For a fluid mixture of N ≥ 2 chemical species A 1 , . . . , A N , it assumes the form of a volume constraint whereV 1 , . . . ,V N > 0 are partial specific volumes of the molecules at reference temperature and pressure. The relation generalises the assumption of a constant mass density considered in other analytical investigations, a. o. [CJ15], [MT15], [HMPW17], [BP17]. In the present paper we are interested only in the general case that at least two indices exist such that V i 1 =V i 2 or, in vectorial notation, thatV = λ 1 N for all λ ∈ R, whereV = (V 1 ,V 2 , . . . ,V N ) and 1 N = (1, 1, . . . , 1) ∈ R N . Bulk. The convective and diffusive mass transport of these species and the momentum balance are described by the partial differential equations ∂ t ρ i + div(ρ i v + J i ) = r i for i = 1, . . . , N , (3) The physical system is assumed isothermal with absolute temperature θ > 0. The partial mass densities of the species are denoted ρ 1 , . . . , ρ N . Throughout the paper we shall use the abbreviation ̺ := N i=1 ρ i for the total mass density. The barycentric velocity of the fluid is called v and the thermodynamic pressure p. In the Navier-Stokes equations, S(∇v) denotes the viscous stress tensor, which we assume for simplicity of Newtonian form. The vector fields b 1 , . . . , b N are the external body forces. The diffusions fluxes J 1 , . . . , J N , that are defined to be the non-convective part of the mass fluxes, must satisfy by definition the necessary side-condition N i=1 J i = 0. A thermodynamic consistent Fick-Onsager closure respecting this constraint is assumed. This approach is described in great generality among others by [BD15], [DGM18] following older ideas by [MR59], [dM63]. The diffusions fluxes J 1 , . . . , J N obey M i,j (ρ 1 , . . . , ρ N ) (∇µ j − b j ) for i = 1, . . . , N . (4) The Onsager matrix M(ρ 1 , . . . , ρ N ) is a symmetric, positive semi-definite N × N matrix for every (ρ 1 , . . . , ρ N ) ∈ R N + . In all known linear closure approaches, this matrix satisfies N i=1 M i,j (ρ 1 , . . . , ρ N ) = 0 for all (ρ 1 , . . . , ρ N ) ∈ R N + .
One possibility to compute the special form of M is for instance to invert the Maxwell-Stefan balance equations. For the mathematical treatment of this algebraic system, the reader can consult [Gio99], [Bot11], [JS13], [HMPW17], [MT15]. Or M is constructed directly in the form P T M 0 P , where M 0 is a given matrix of full rank, and P is a projector guaranteeing that (5) is valid. The paper [BDc] establishes equivalence relations between the Fick-Onsager and the Maxwell-Stefan constitutive approaches, proposing moreover a novel unifying approach to close the diffusion model. The quantities µ 1 , . . . , µ N are the chemical potentials from which the thermodynamic driving forces for the diffusion phenomena are inferred. For an incompressible system, they are related to the mass densities ρ 1 , . . . , ρ N and to the pressure via Here the function k denotes the the positively homogeneous part of the free energy, which is independent of thermodynamical pressure. A typical choice discussed in [BD15] is where n i := ρ i /m i are the number densities with the molecular masses m 1 , . . . , m N > 0, y i = n i / N j=1 n j are the number fractions, and µ ref i are reference values of the chemical potentials. For the mathematical theory in this paper, more general structures in (7) will however be admitted. The isothermal Gibbs-Duhem equation: dp = N i=1 ρ i dµ i defines the intrinsic relationship between (1), (6), and the pressure field. The paper [BDa] shows that the relation (6) indeed occurs in the limit case when the bulk free energy density of the system adopts the singular form The relation (6) is an equivalent expression of µ ∈ ∂h ∞ (θ, ρ), where ∂ denote the subdifferential of the convex function h ∞ (θ, ·), and the function p = −h ∞ (θ, ρ) + N i=1 ρ i µ i can be understood as a 'Lagrange multiplier' associated with the constraint (1).
We notice that, multiplying the equations (2) with the constantsV i and summing up, the local change of volume is described by the equation Effects like diffusion and chemical reactions will induce a local change in the molecular composition, implying a net local change of the volume, independent of a mechanical compression or expansion. Concerning the presence of reaction terms in (2), we have to mention in respect with the compressible systems considered in [BDb] a subtle difference of the incompressible models. For the compressible case, the reactions densities r i in (2) are allowed to be general functions r i = r i (ρ 1 , . . . , ρ N ), without influencing qualitatively the well-posedness results or the mathematical methods. This is different in the incompressible case. At first, the restriction (1) implies that µ does not depend on ρ only, so that the structure r = r(ρ) does not comply with standard thermodynamically consistent reaction terms. At second, the 'elliptic equation' (9) defines a differential operator acting on a certain relative chemical potential (variable ζ, details below). This elliptic operator is linear for the pure diffusion case, but turns to non-linear in the presence of reactions of the general form r = r(µ). In this paper, we treat incompressible multicomponent diffusion in itself. We shall address the specific problems raised by chemical reactions in further research. Thus, allowing -as we shall do -for certain source terms r = r(ρ) in (2) means a bit more mathematical generality, but it remains clear that realistic models of chemical reactions require non trivial modifications of the methods used here. As to the stress tensor S we shall restrict for simplicity to the standard Newtonian form with constant coefficients. We, however, present methods which are sufficient to extend the results to the case of density and composition dependent viscosity coefficients. Boundary and initial conditions. We investigate the problem (2), (3) in a cylindrical domain Q T := Ω×]0, T [ where T is a finite time and Ω ⊂ R 3 a bounded domain. It is possible to treat the case Ω ⊂ R d for general d ≥ 2 with similar methods. We consider initial conditions v j (x, 0) = v 0 j (x) for x ∈ Ω, j = 1, 2, 3 .
As a matter of fact, these simplifying choices oblige us to make a further restriction. To see this, we recall the relation (9), that we integrate over Ω. If there is no mass flux through the boundary, we see that´Ω N i=1V i r i (x, t) dx = 0. This condition cannot be enforced for a general r = r(ρ), unless we assume that r takes values in {V } ⊥ . Recalling that realistic models for chemical reactions are to be treated in an upcoming paper, we here restrict to the case that r(ρ) ·V = 0 for all ρ.
2. State of the art and our main result 2.1. A review of prior investigations and our method. Up to few exceptions, models for incompressible multicomponent fluids have not been investigated in mathematical analysis. For a mathematical treatment in the case of the constraint ̺ = const, which corresponds to choosingV 1 = . . . =V N in (1), the reader might consult the papers [CJ15] and [MT15] (global weak solution analysis) and [HMPW17], and [BP17] (local-in-time well-posedness). 2 From the viewpoint of the mathematical structure, the case ̺ = const exhibits profoundly different features than the general relation (1). The principal difference is that (4), (5) and (6) imply the decoupling of the pressure and of the diffusion fluxes. The Navier-Stokes equations reduce to their single component solenoidal variant and can be solved independently. Of course, this does not mean that ̺ = const cannot be a good approximation under special circumstances. In [BS16] for instance, a class of multicomponent mixtures has been introduced for which the use of the incompressible Navier-Stokes equation is realistic: Incompressibility is assumed for the solvent only, and diffusion is considered against the solvent velocity. See also the discussion in the paragraph 4.8 of [PS14] on incompressible mixtures. In the case thatV is not parallel to 1 N , (4) implies that the pressure affects the diffusion fluxes via the chemical potentials. A corollary of this fact is that if we multiply the equations (2) with the constantsV i and sum up, we obtain (9) for the local change of volume. Moreover, (a) the viscous stress tensor does not simplify to the symmetric velocity gradient; (b) the total mass density is calculated from the continuity equation ∂ t ̺ + div(̺ v) = 0; (c) the pressure remains partly connected to the other variables by an algebraic formula. Our main method to approach the PDE problem is a switch of variables in the transport problem as already applied in [BDb]. Instead of the original variables (ρ 1 , . . . , ρ N ) and (p, v 1 , v 2 , v 3 ), we regard N − 1 linear combinations of the chemical potentials (µ 1 , . . . , µ N ), the mass density ̺ and the velocity field as main variables. After the transformation we obtain for the new free variables (̺, q 1 , . . . , q N −2 , ζ, v 1 , v 2 , v 3 ) -instead of (2), (3) -the equations (here without external forcing and chemical reactions) The nonlinear field R and the function P , the vector field A, the positive matrix M , and the positive coefficient function d will be constructed below, combining certain linear projection operators with the inverse map for the algebraic equations µ =V p + ∇ ρ k(ρ). We are then faced with a nonlinear PDE system of mixed parabolic-elliptic-hyperbolic type. All variables are unconstrained, but for the restriction ̺ min < ̺ < ̺ max on the total mass density. Here, the constants 0 < ̺ min < ̺ max < +∞ are the thresholds of the total mass for states ρ 1 , . . . , ρ N that satisfy the constraint (1): Comparing with the paper [BDb] on compressible class-one models based on a similar reformulation, we see that the incompressible limit corresponds structurally to the case that one of the relative chemical potentials is subject to an elliptic -instead of a parabolic -equation, and the total mass density is confined to a bounded interval.
For an overview of possible methods to study the transformed PDE system, we refer to our study [BDb]. We shall follow here the same principal road map, but profound transformations are necessary to deal with the constraint on ̺, since it implies that the nonlinear functions occurring in the transformed system are singular for dist(̺, {̺ min , ̺ max }) → 0. The solution operator to the continuity equation, however, does not 'see' these thresholds, which is the source of additional problems when we attempt to linearise. Moreover, we must construct a solution operator for the parabolic-elliptic subsystem of general form for (q, ζ), while the reduced transport problem in [BDb] was purely parabolic. Nontrivial extensions of the method are therefore necessary to deal with the incompressible case. We shall study the problem in the class proposed in the paper [Sol80] for Navier-Stokes: W 2,1 p with p larger than the space dimension for the velocity and W 1,1 p,∞ for the density. For the variable q 1 , . . . , q N −2 , we also choose the parabolic setting of W 2,1 p . For the elliptic component ζ, we choose the state space W 2,0 p . In these classes, we are able to prove the local existence for strong solutions. In general, we obtain only a short-time well-posedness result, and boundedness in the state space is not sufficient to guarantee that the solution can be extended to a larger time interval. This is due to the constraint ̺ min < ̺ < ̺ max : A strong solution with bounded state space norm might break down if the density reaches the thresholds. However, it is to note that for choices of the tensor M reflecting the physically expected behaviour that, in the dilute limit, a diffusion flux is linearly proportional with the mass density of the vanishing species, we are able to show that a sufficiently smooth solution (p > 5) bounded in the state space cannot reach the critical values in finite time. Thus, a kind of maximum principle is available for the system. We shall also prove the global existence under the condition that the initial data are sufficiently near to an equilibrium (stationary) solution. However, since this result relies on stability estimates in the state space, we need to assume higher regularity of the initial data in order to obtain some stability from the continuity equation. Therefore, these solutions exist on arbitrary large time intervals, but do not enjoy the extension property.
We shall not make use of the Lagrangian coordinates but employ the approach of controlled growth in time of the solution by means of a priori estimates. Let us finally mention also the paper [FLM16], devoted to binary mixtures. Starting from different modelling principles in the spirit of [JHH96], the authors derive for N = 2 a similar PDE system. The variable q does not occur, and the coefficient d is assumed constant. The authors prove for this system the global existence of weak solutions if the singularity of P (̺) at the thresholds is sufficiently strong. The weak solution analysis for the general system is considered in the paper [Dru19].
2.2. Main results. We denote Q = Q T = Ω×]0, T [ with a bounded domain Ω ⊂ R 3 and T > 0 a finite time. We use the standard Sobolev spaces W m,p (Ω) for m ∈ N and 1 ≤ p ≤ +∞, and the Sobolev-Slobodecki spaces W s p (Ω) for s > 0 non-integer. If Ω is a domain of class C 2 , the spaces W s p (∂Ω) are well defined for 0 ≤ s ≤ 2. With a further index 1 ≤ r ≤ +∞, we use the parabolic Lebesgue spaces L p,r (Q) (space index first: L p (Q) = L p,p (Q)). For ℓ = 1, 2, . . . and 1 ≤ p ≤ +∞ we introduce the parabolic Sobolev spaces and, with a further index 1 ≤ r < ∞, the spaces In these notations, the space integrability index always comes first. For r = +∞, W ℓ,ℓ p,∞ (Q) denotes the closure of C ℓ (Q) with respect to the norm above and, thus, We also encounter, for ℓ = 1, 2 and 1 ≤ p < +∞, We denote by C(Q) = C 0,0 (Q) the space of continuous functions over Q and, for α, β ∈ [0, 1], the Hölder spaces are defined by C α, β (Q) : Some brief remarks on notation: (1) All Hölder continuity properties are global. For the sake of notation we identify C α, β (Q) with C α, β (Q).
(2) Whenever confusion is impossible, we shall also employ for a function f of the variables x ∈ Ω and t ≥ 0 the notations f x = ∇f for the spatial gradient, and f t for the time derivative.
(3) For maps like R, M which depend on ̺ and q, the derivatives are denoted by R ̺ , M q . Due to (5), the matrix M(ρ) possesses only N −1 positive eigenvalues that moreover might degenerate for vanishing species. The orthogonal projection on the N −1 dimensional linear The vectorV occurring in (1) defines another singular direction in the model preventing parabolicity. We denote by We also introduce the notations The surface SV is the domain of existence for the incompressible state. It is readily seen that ρ ∈ SV implies for the variable ̺ := N i=1 ρ i the inequalities ̺ min = 1 max j=1,...,NVj < ̺ < ̺ max = 1 min j=1,...,NVj for all ρ ∈ SV .
Our first main result is devoted to the short-time existence of a strong solution. (In order to avoid notational confusion with the pressure field, the integrability index is called s in the next statements.) Theorem 2.1. We fix s > 3 and T > 0 and assume that (a) Ω ⊂ R 3 is a bounded domain of class C 2 ; , positively homogeneous, convex in its domain R N + , and lim inf m→+∞ |∇ ρ k(y m )| = +∞ for all sequences {y m } ⊂ S 1 approaching the relative boundary of S 1 ; (d) r : . For simplicity, we assume ν(x) · P {1 N } ⊥ b(x, t) = 0 for x ∈ ∂Ω and λ 1 −almost all t ∈]0, T [. (f ) The initial data ρ 0 1 , . . . ρ 0 N : Ω → SV are positive measurable functions satisfying the following conditions: • The initial total mass density ̺ 0 : Then, there exists T * ∈ (0, T ] such that the problem (2), (3) with closure relations (4), (6), incompressibility constraint (1) and boundary conditions (10), (11), (12), (13) possesses a unique solution (ρ, p, v) of class The solution can be uniquely extended to a larger time interval whenever the two following conditions are fulfilled: (ii) There is α > 0 such that the quantity stays finite as t ր T * . Here z = z(s) is defined via z = 3/(s − 2) for 3 < s < 5, z > 1 arbitrary for s = 5 and z = 1 if s > 5.
It is to note that the possibility to extend the solution is not -like in the compressible case -reducible to the smoothness criterion (ii). If (i) is failing, even a smooth solution can break down if its total mass density reaches the critical values {̺ min , ̺ max }. This singularity plays an important role also in the context of the weak solution analysis (see [Dru19]). However, we provide an important complement for physically motivated choices of the mobility matrix M and of the function k. Here the boundedness in the natural state space norm is sufficient to guarantee the extension property.
Theorem 2.2. In the situation of Theorem 2.1 we assume, in addition, that s > 5 and that k is the function defined in (7). We define a matrix B i,j (ρ) := M i,j (ρ)/ρ j for i, j = 1, . . . , N, and we assume that there is a continuous function C = C(|ρ|), bounded on compact subsets of R N + \ {0}, such that |B i,j (ρ)| + ρ k |∂ ρ k B i,j (ρ)| ≤ C(|̺|) for all i, j, k ∈ {1, . . . , N} and all ρ ∈ R N + . Then the strong solution of Theorem 2.1 can be extended beyond T * whenever Our second main result concerns global existence under suitable restrictions on the data. An equilibrium solution for (2), (3) is defined as a vector (ρ eq 1 , . . . , ρ eq N , p eq , v eq 1 , v eq 2 , v eq 3 ) of functions defined in Ω such that ρ eq ∈ W 1,s (Ω; SV ), p eq ∈ W 1,s (Ω), v eq ∈ W 2,s (Ω; R 3 ) , the vector µ eq := p eqV + ∇ ρ k(θ, ρ eq ) satisfies P {1 N } ⊥ µ eq ∈ W 2,s (Ω; R N ) and the relations div(ρ eq and div(̺ eq v eq ⊗ v eq − S(∇v eq )) are valid in Ω. The boundary conditions are v eq = 0 and ν(x) · M i,j (ρ eq ) (∇µ eq j − b j (x)) = 0 on ∂Ω . We show that the problem (2), (3) possesses a unique strong solution on an arbitrary large, but finite time interval if the distance of the initial data to an equilibrium solution is sufficiently small, and if both initial conditions and equilibrium solution are smooth enough.

Road map.
In sections 3 and 4 we show how to reformulate the original system such that it becomes easier to tackle via functional analytic methods. The functional setting is discussed in section 5. In section 6, we introduce two ways to linearise the PDE system and reformulate the initial-boundary-value problem as a fixed point problem in the state space. Both fixed point equations exploit the parabolic substructure for the variables (q, v) and treat the linear equations for (ζ, ̺) as side conditions. In the first method, used to prove the short-time well posedness, all lower-order nonlinearities are frozen. For the proof of Theorem 2.3 on small perturbations, a somewhat more elaborated linearisation principle is used in order to exhibit some stability estimates.
The estimates for the linearised principal part of the system are presented in section 7.
Here we can rely partly on our work in [BDb] for the compressible system, but have to discuss the additional problems caused by the presence of an elliptic equation and of a density constraint in the continuity equation. Section 8 shows the self mapping estimate for the first fixed point equation, which yields the well posedness result in section 9. The extension criteria proved for the solution in the same section 9 deserve attention in their own right. The proof of the global well-posedness result for small data, or rather small perturbations, is given in section 10. Finally, some reminder, tools, and purely technical statements are compiled in the Appendix.

The singular free energy function and its conjugate
In comparison to the analysis of compressible models in [BDb], a main specificity of the incompressible model concerns the bulk free energy density and the definition (6) of the chemical potentials. With k : R N + → R given, we introduce a bulk free energy density defined for ρ ∈ R N + of the form The function h ∞ is singular, but the subdifferential ∂h ∞ is non-empty for every ρ satisfying the incompressiblity constraint N i=1 ρ iVi = 1. If the function k is continuously differentiable, it can be shown that µ ∈ ∂h ∞ (ρ) if and only if there exists p ∈ R such that µ i = pV i + ∂ ρ i k(ρ) for i = 1, . . . , N. It can easily be verified that the number p can be characterised as follows: where (h ∞ ) * is the convex conjugate of h ∞ . For systematic discussions and a proof of these elementary statements, we refer to [BDa]. Our approach essentially relies on the properties of the dual free energy function f := (h ∞ ) * on R N . We shall recall three statements of the paper [BDa]. Proofs are provided in the Appendix, Section A for the reader's convenience. In the special case that the gradient of k is explicitly invertible on S 1 (see (14)), the statements can also be proved by direct algebraic computations yielding in many cases explicit formulae; see the Section 4 in [Dru19] for a complete characterisation of the example (7).
Lemma 3.1. We assume that k : R N + → R is a positively homogeneous convex function of class C 3 (R N + ). We moreover assume that the restriction of k to the surface S 1 is essentially smooth, meaning that |∇ ρ k(y m )| → +∞ for sequences {y m } m∈N ⊂ S 1 such that min i=1,...,N y m i → 0 as m → +∞. For µ ∈ R N , we define f (µ) := sup ρ∈SV {µ · ρ − k(ρ)}. Then the function f belongs to C 3 (R N ), and ∇ µ f maps onto SV .

Change of variables for the incompressible model
We propose a reformulation of the equations (2), (3) subject to the constitutive equations (4), (6) and to the volume constraint (1) in order to eliminate the positivity constraints on ρ, the singularity due to M 1 N = 0 (cf. (5)), and the singularity direction due to the incompressibility (1) -equivalently, the fact that the function f , interpreted as the dual of the free energy, is affine in the direction ofV (D 2 fV = 0, Lemma 3.2). Like in the investigations in [DDGG17], [BDb], [Dru19], the idea is to invert the algebraic relations (6) for µ, p, ρ and to combine this procedure with appropriate linear projections.
4.1. General ideas. We choose a basis of R N : {ξ 1 , . . . , ξ N −2 , ξ N −1 , ξ N } with ξ N = 1 N and ξ N −1 =V . We then choose η 1 , . . . , η N to be the dual basis, i. e. ξ i · η j = δ i j for i, j = 1, . . . , N. We define variables q 1 , . . . , q N −2 and ζ via We exploit the result of Lemma 3.3 saying that (6) implies ρ i = ∂ µ i f (µ 1 , . . . , µ N ) for i = 1, . . . , N. The vector µ is then decomposed according to expressed by the variables q and ζ, and its projection on span{1 N }. Next, the last coordinate µ · η N is eliminated using the equation The gradient ∇ µ f is invariant in the directionV (cf. Lemma 3.2) and, therefore, the variable ζ decouples from the latter equation, that now reads This representation is an algebraic equation F (µ · η N , q 1 , . . . , q N −2 , ̺) = 0. In view of Lemma 3.2, note that ∂ µ·η N F (µ · η N , q 1 , . . . , q N −2 , ̺) = −D 2 f (µ)1 N · 1 N < 0, due the fact that 1 N is not parallel toV . Thus, the last component µ · η N is defined implicitly as a differentiable function of ̺ and q. We call this function M and obtain the equivalent formulation where only the total mass density ̺ and the relative chemical potentials q 1 , . . . , q N −2 and ζ occur as free variables. Note, moreover, that ζ and ρ decouple. Similarly, we obtain a representation of the pressure as All this is summarised in the following Lemma, the proof of which is direct in view of the Lemmas 3.1 and 3.3.
Lemma 4.1. We adopt the assumptions of Theorem 2.1 for the function k.
In order to deal with the right-hand side (external forcing), we define in the same spirit:

This allows to express
For the reaction terms, we definer ℓ (̺, q) : 4.2. Reformulation of the partial differential equations and of the main theorem.
The relation (5) and the equivalence of Lemma 4.1 show that If we introduce the rectangular projection matrix Π j,ℓ = ξ ℓ j for ℓ = 1, . . . , N − 2 and j = 1, . . . , N, In the latter system, we have ρ = R(̺, q) and (̺, q 1 , . . . , q N −2 , ζ, v 1 , v 2 , v 3 ) are the independent variables. Next, we define for k = 1, . . . , N − 2 the maps Multiplying the mass transfer equations with ξ k i , we obtain that It can be checked easily that the matrix Π T M(ρ) Π ∈ R (N −2)×(N −2) is symmetric and strictly positive definite on all states ρ ∈ SV . The Jacobian is also strictly positive definite. Indeed, vectors of the form Π a in R N with nonzero a ∈ R N −2 can by construction never belong to span{1 N ,V }. We next multiply the mass balance equations withV i . Making use of the constraint (1) yields where we use the additional assumption that r maps into {V } ⊥ . Using that ρ = R(̺, q), we define Overall, we get for the variables (̺, q 1 , . . . , q N −2 , ζ, v) -instead of (2), (3) -the equations The problem (P ′ ) consisting of (25), (26), (27) and (28) for the variables (̺, q, ζ, v) might seem to exhibit more nonlinearities than the original problem for ρ, p and v. However, it has the advantage that -up to the restriction on the total mass density ̺ min < ̺ < ̺ max -it is completely free of constraints. Furthermore, the differential operator is linear in the variable ζ, which occurs only under spatial differentiation. Our first aim is now to show that, at least locally in time, the system (25), (26), (27) and (28) for the variables (̺, q 1 , . . . , q N −2 , ζ, v) is well posed. We consider initial conditions Due to the preliminary considerations in section 4.1, prescribing these variables is completely equivalent to prescribing initial values for the mass densities ρ i and the velocity. It suffices to define q 0 k = η k · ∂ ρ k(ρ 0 ) for k = 1, . . . , N − 2. For simplicity, we consider the linear homogeneous boundary conditions The conditions (31) and (13) are equivalent, because we assume throughout that the given forcing b satisfies ν(x) · P {1 N } ⊥ b(x, t) = 0 for x ∈ ∂Ω (see assumption (e) in the statement of Theorem 2.1). Under the assumptions of Theorem 2.1 for the function k, the coefficient functions R, M , A, d and P are of class C 2 in the domain of definitions I ×R N −2 as shown in the Lemma 4.1. We reformulate the Theorem 2.1 for the new variables. Since the thermodynamic pressure does not occur explicitly as a variable, we now switch to denoting p > 3 the integrability exponent (denoted s in the statement 2.1).
Theorem 4.2. Assume that the coefficient functions R, M , A, d and P are of class C 2 , andr is of class C 1 in the domain of definition I × R N −2 . Let Ω be a bounded domain with boundary ∂Ω of class C 2 . Suppose that, for some p > 3, the initial data are of class Then there is 0 < T * ≤ T , depending only on these data, such that the problem (25), (26), (27) and (28) with boundary conditions (29), (30) and (31) is uniquely solvable in the class . The solution can be uniquely extended within this class to a larger time interval whenever at least one of the following holds: (1) p > 5 and the state space norm stays finite as t ր T * ; (2) The two following conditions are valid as with α > 0 and z = z(p) defined by Theorem 2.1,

Functional analytic approach
Recall thatb,b andb are given coefficients. To get rid of the highest-order coupling in the time derivative of ̺, we shall employ the same approach as in [BDb], which is sketched below. Consider a solution u = (q, ζ, ̺, v) to A (u) = 0. Computing time derivatives in the equation A 1 (u) = 0, we obtain that Here the nonlinear functions R, R ̺ , R q , A and M ,r etc. are evaluated at (̺, q). Under the side-condition the first component being the differential operator defined by (32). Clearly, A (u) = 0 if and only if A (u) = 0. The functional setting was introduced in Section 2.2. Similar spaces were used in [BDb] to study the compressible system and, in order to save room, we shall refer to this paper for the trace and embedding theorems needed in the present analysis. For p > 3 and α := 1/2 + 3/(2p), we recall the interpolation inequality (see [Nir66], Theorem 1) valid for any function f in W 2,p (Ω), with certain constants C 1 , C 2 depending only on Ω. We consider the operator (q, ζ, ̺, v) → A (q, ζ, ̺, v) acting on The natural trace space at time zero is denoted Tr Ω×{0} X T . The functional setting does not allow to introduce traces for the variable ζ. Therefore, u(0) ∈ Tr Ω×{0} X T means that (q(0), We denote by 0 X T the space of functions fulfilling zero initial conditions. This only makes sense, of course, for the variables having traces at Ω × {0}. Thus Since the coefficients of A are defined only if ̺ has range in I, the domain of the operator is contained in the subset We shall moreover make use of a reduced state space containing only the parabolic components (q, v), namely The operator A is the composition of differentiation, multiplication and Nemicki operators. Therefore, the properties of the coefficients R,M etc. allow to show that A is continuous and bounded from X T,I into Since the coefficients R, M , A, d and P are twice continuously differentiable in their domain of definition I × R N −2 , the operator A is even continuously differentiable at every point of X T,I . We spare the proof of these rather obvious statements.
6. Linearisation and reformulation as a fixed-point equation We shall present two different manners to linearise the equation A (u) = 0 for u ∈ X T with initial condition u(0) = u 0 in Tr Ω×{0} X T . They correspond to the two main Theorems 2.1, 2.3 respectively. In both cases, we start considering the problem to find u = (q, ζ, ̺, v) ∈ X T,I such that A (u) = 0 and u(0) = u 0 , which after permuting rows, possesses the following structure The functions g, h and f stand for the following expressions: These expressions are independent on the component ζ. We can regard g, h and f as functions of x, t and of the vectors u and D x u and write g(x, t, u, D x u).
6.1. The first fixed-point equation. For (q * , v * ) given in W 2,1 p (Q T ; R N −2 )×W 2,1 p (Q T ; R 3 ) and for unknowns u = (q, ζ, ̺, v), we consider the following system of equations together with the initial conditions (29), (29), (29) and the homogeneous boundary conditions (30), (31). Note that the continuity equation can be solved independently for ̺. Once ̺ is given, we solve the linear parabolic-elliptic system (47), (48) for q and ζ. Here we must be careful, since the coefficients of this system are only defined as long as ̺(x, t) takes values in I. Thus, the solution (q, ζ) might exist only on a shorter time interval. We can solve the problem (49), which is linear in v, under the same restriction.
We will show that the solution map (q * , v * ) → (q, v), denoted T , is well defined from Y T into itself for T fixed and suitably small. The solutions are unique in the class Y T . Clearly, a fixed point of T is a solution to A (q, ζ, ̺, v) = 0.
6.2. The second fixed-point equation. Here we construct the fixed-point map comparing the solutions to a given reference vector (q 0 ,v 0 ) ∈ Y T that extends the initial data. We assume thatq 0 andv 0 satisfy the initial compatibility conditions. In order to find an extension for ̺ 0 ∈ W 1,p (Ω), we solve the problem For this problem, Theorem 2 of [Sol80] establishes unique solvability in W 1,1 p,∞ (Q T ) and, in particular, the strict positivity̺ 0 ≥ c 0 (Ω, v 0 We find the extensionζ 0 by solving, for all values of t such that the coefficientsb(t) and b(t) are defined, the elliptic problem with homogeneous Neumann boundary conditions and zero mean-value side-condition. Consider a solution u = (q, ζ, ̺, v) ∈ X T to A (u) = 0. We introduce the differences r := q −q 0 , χ = ζ −ζ 0 , w := v −v 0 and σ := ̺ −̺ 0 , and their vectorū := (r, χ, σ, w). Clearly,ū belongs to the space 0 X T of homogeneous initial conditions. Recall that this does not imply a trace condition for χ, cp. (35). The equations A (u) = 0 mean, equivalently, that A (û 0 +ū) = 0. The vectorū = (r, χ, σ, w) satisfies Herein, all non-linear coefficients R, R q , etc. are evaluated at (̺, q), while g, h and f correspond to (43), (44) and (45). We next want to construct a fixed-point map to solve (52), (53), (54), (55) by linearising g 1 , h 1 and f 1 defined in (52), (53) and (55). First, we expand as follows: Here, (·) θ applied to a function of x, t, u and D 1 x u stands for the evaluation at ( In short, in order to avoid the integral and the parameter θ, we write Obviously, the latter expressions make sense only if u, u * both belong to X T,I , in which case the entire convex hull {θ u + (1 − θ) u * : θ ∈ [0, 1]} is in X T,I . Following the same scheme as for (57), we write in short Similar expressions are obtained for h 1 and f 1 . In the case of h 1 , note however that x ) = 0 due to the construction (51) ofζ 0 . Now we construct the fixed-point map to solve (52), (53), (54) and (55). For a given vector (r * , w * ) ∈ 0 Y T , we define q * :=q 0 + r * and v * :=v 0 + w * . Then we define ̺ * to be the unique solution to We thus write ̺ * := C (v * ) where C is the solution operator to the continuity equation with initial data ̺ 0 . We employ the abbreviation Forū := (r, χ, σ, w), we next consider the linear problem with the boundary conditions ν · ∇r = 0 = ν · ∇χ on S T and w = 0 on S T , and with zero initial conditions for r, σ and w. The superscript * on a coefficient means evaluation at (C (v * ), q * ). We will show that the solution map T 1 : (r * , w * ) → (r, w) is well defined from 0 Y T into itself for T > 0 arbitrary, provided that the distance of the initial data to an equilibrium solution is sufficiently small. As to the latter restriction, note that the expressions (g 1 ) ′ (u * ,û 0 ) make sense only if the density components in both u * andû 0 map into the interior of the critical interval, which cannot be expected globally for the solutions to (50) and (59). Ifū = (r, w) is a fixed point of T 1 , then we can show that u :=û 0 +ū is a solution to A (u) = 0. This is verified exactly as in [BDb], Remark 6.1. 6.3. The self-mapping property. Assume that the map T : (q * , v * ) → (q, v) via the solution to (46), (47), (48), (49) is well defined in Y T , with image in YT for someT = T (q * , v * ) > 0. Then, we want to show that T maps some closed bounded set of Y T 0 into itself for a fixed T 0 > 0. Here, a major change occurs in comparison to the compressible case, since we do not expect that the linearised map T produces a solution defined globally up to T . This is due to the constraint ̺ ∈]̺ min , ̺ max [ which can by nature be enforced only locally for solutions to the continuity equation (46). We shall rely on continuous estimates expressing the controlled growth of the solution in time. We will show that there is a parameter a 0 depending on the distance of the initial density to the singular values {̺ min , ̺ max } such that, whenever t > 0 satisfies t 1− 1 p (q * , v * ) Yt < a 0 , the pair (q, v) = T (q * , v * ) is well defined in Y t and satisfies the estimate Here R 0 stands for the magnitude of the initial data q 0 , ̺ 0 and v 0 , and of the external forces b in their respective norms. The function Ψ is continuous, increasing in all arguments, and finite for t 1− 1 p (q * , v * ) Yt < a 0 . Hence we obtain a self mapping property with the help of the following Lemma.

The inequality (65) is valid by assumption and it yields
In the case of the map T 1 : (r * , w * ) → (r, w) defined via solution to (59), We will show that this can be ensured if the starting perturbation w * satisfies an inequality of type in which a 0 > 0 is a fixed number depending on the distance of the initial data to the critical values {̺ min , ̺ max }, and φ 0 is a continuous function on R 2 + , which increases in both arguments. We then prove a continuity estimate of the type Here R 0 stands for the magnitude of initial data (q 0 , ̺ 0 and v 0 ) and external forces b. The parameter R 1 expresses the distance of the initial data to a stationary/equilibrium solution (def. in (16), (17)). Defining η 0 to be the smallest positive solution to the equation φ 0 (T, η 0 ) η 0 = a 0 , we will show that T 1 maps the ball of radius η 0 in 0 Y T for initial data satisfying R 1 ≤ η 0 /Ψ(T, R 0 , η 0 ). In order to apply the contraction principle and prove the theorems, we shall therefore prove the continuity estimate (65), (66). This is the main object of the next sections.

Estimates of the linearised problems
In this section, we present the estimates on which our main results in Theorem 2.1, 4.2 are footing. The preliminary work done in the paper [BDb] shall, in many points, allow to abridge the calculations. The main novelty is the inversion of the parabolic-elliptic subsystem, which shall be dealt with in all details. To achieve more simplicity in the notation, we introduce both for a function or vector field f ∈ W 2,1 p (Q T ; R k ) (k ∈ N) and t ≤ T the notation Recall that W 2−2/p p (Ω) is the trace space for f ∈ W 2,1 p (Q T ), f → f (·, t). Moreover we will need Hölder half-norms. For α, β ∈ [0, 1] and f scalar valued, we denote The corresponding Hölder norms f C α (Ω) , f C α (0,T ) and f ∈ C α,β (Q T ) are defined by adding the corresponding L ∞ −norm to the half-norm.
7.1. Estimates of a linearised problem for the variables q and ζ. We first formulate some global assumptions and notations. Recall that I =]̺ min , ̺ max [. In this section, the maps R q , M : I × R N −2 → R (N −2)×(N −2) are assumed to be of class C 1 (I × R N −2 ) into the set of symmetric, positive definite matrices. Furtheron, A : I × R N −2 → R N −2 , and d : I × R N −2 → R + are of class C 1 too. We fix p > 3, and we consider given q * ∈ W 2,1 p (Q T ; R N −2 ) and ̺ * ∈ W 1,1 p,∞ (Q T ) such that ̺ * (x, t) ∈]̺ min , ̺ max [ for all (x, t) ∈ Q T .
We then denote R * q := R q (̺ * , q * ), M * := M (̺ * , q * ), A * := A(̺ * , q * ) and d * := d(̺ * , q * ). For t ≤ T , we introduce the positive functions We let g ∈ L p (Q T ; R N −2 ), q 0 ∈ W 2−2/p p (Ω; R N −2 ) such that ν · ∇q 0 (x) = 0 on ∂Ω in the sense of traces, and h ∈ W 1,0 p (Q T ; R 3 ). For a pair (q, ζ) : Q T → R N −2 × R we consider the linear parabolic-elliptic auxiliary problem and we want to obtain an estimate in the norm of W 2,1 p (Q T ; R N −2 ) × W 2,0 p (Q T ) for the solution. To this aim we first show that (70), (71) can be equivalently reformulated as a system coupled only in the lower order.
Lemma 7.1. We adopt the general assumptions and notations formulated at the beginning of this section. A pair (q, ζ) ∈ W 2,1 is a solution to the problem (70), (71) if the identity (71) and the initial and boundary condition are satisfied, and if instead of (70) we have Computing the derivatives in the elliptic equation (71), we obtain that Thus, under the side-condition (73), the parabolic equations (70) are equivalent to yields the claim. Using this lemma, we next prove an estimate for the solution to the linearised parabolicelliptic problem.  Ψ 1 (t, a 1 , . . . , a 5 ) and Φ = Φ(t, a 1 , . . . , a 5 ) defined for all t ≥ 0 and all numbers a 1 , . . . , a 5 ≥ 0, such that for all t ≤ T and for 0 < β ≤ 1 arbitrary: , ∇̺ * L p,∞ (Qt) ) .
The function Ψ 1 possesses moreover the following two properties: It is increasing in all arguments, and the value of Ψ 1 (0, a 1 , . . . , a 5 ) = Ψ 0 1 (a 1 , a 2 ) is a function independent on the three last arguments. The function Φ is increasing in all arguments.
Proof. The existence and uniqueness can be easily obtained by means of the uniform estimates. We thus suppose first that (q, ζ) ∈ W 2,1 p (Q T ; R N −2 ) × W 2,0 p (Q T ) is a given solution, and we prove the claimed estimate. In order to simplify the discussion, we adopt the following convention: When computing the derivative of a coefficient, like ∇d ) a generic continuous function depending only on M * (t) and q * L ∞ (Qt) , and increasing in these arguments. We then bound the L ∞ (Q T ) norms of all non-linear functions depending on ̺ * , q * by this generic c * 1 .
Step 1: First estimate for the variable ζ.
For almost all s ≤ t, the function ζ satisfies the weak Neumann problem By well-known weak elliptic theory, there is a unique solution ζ(s) ∈ W 1,p (Ω) with Ω ζ(x, s) dx = 0. Moreover, for all 0 < β < 1, perturbation techniques shortly recalled in the Appendix, Lemma B.5 yield the estimate We define φ * t := sup s≤t (1 + [d * (s)] C β (Ω) ) 1 β . We bound sup s≤t A * (s) L ∞ (Ω) with a generic c * 1 , and it follows that Step 2: First bound for the variable q.
Step 3: Main estimate for the variable ζ. Since ζ ∈ W 2,0 p (Q T ), we can employ the pointwise identity (73). Since ζ has mean-value zero for all times, the full W 2,p norm can be estimated by the Neumann-Laplacian, and we obtain that Computing the derivatives of the coefficients, and using the same conventions as above, we derive from (79) the inequality ζ W 2,0 p (Qt) ≤ c * 1 ( △q L p (Qt) + div h L p (Qt) ) + c * 1 ( ∇̺ * · ∇ζ L p (Qt) + ∇q * · ∇ζ L p (Qt) + ∇̺ * · ∇q L p (Qt) + ∇q * · ∇q L p (Qt) ) . We estimate △q L p (Qt) ≤ V (t; q), then we employ the inequality (78) to see that Step 4: Combined estimates. We add (78) to (80) to obtain that ) . In order to control the factors on the right-hand, we first apply (33) to find that , for all ǫ > 0 and a, b > 0. By these means, it follows that Here we use the abbreviation | · | r for · L r (Ω) . Just in the same way, we show that We let F * (t) := sup s≤t ( ∇q * (s) p L p (Ω) + ∇̺ * (s) p L p (Ω) ) and X * (t; ζ) := ∇̺ * · ∇ζ p L p (Qt) + ∇q * · ∇ζ p L p (Qt) . With the help of (82) and of (83), it follows that We choose ǫ = 2 −2−1/p (c * 1 (1 + 2CΨ 1,t )) −p , where c * 1 , C andΨ 1,t are the numbers occurring in the relation (81). Then . Due to our conventions, we can bound every power of c * 1 and the maximum of 1 and c * 1 again by another such function. Introducing a factor By means of (75), we bound ζ L p (Qt) ≤ c * 1 φ * t ( ∇q L p (Qt) + h L p (Qt) ). Raising (85) to the power 1/p, we show that We insert the latter result into (81), obtaining + c * 1 (1 + 2CΨ 1,t ) ( ∇̺ * · ∇q L p (Qt) + ∇q * · ∇q L p (Qt) ) + C p c * 1 φ * t Φ * 1,t ∇q L p (Qt) . In order to estimate X * (t, q), we apply the same steps as for X * (t, ζ) (cf. (85)). Hence which, after raising to the power 1/p, yields Since q W 2,0 p (Qt) ≤ V (t; q), the latter and (86) imply that . In order to finally get rid of the factors with q on the right-hand side, we introduce We raise (87) to the p th power. We use f (t) ≤ V p (t; q) and q p W 1,0 In this way, we obtain the inequality f (t) ≤ 2 p A(t) + 2 p B(t)´t 0 f (τ ) dτ . Using that A and B are monotone increasing by construction, the Gronwall Lemma yields f (t) ≤ 2 p A(t) exp(2 p t B(t)). In particular, we conclude that Combining the latter with (87), it follows that In order to verify that the factors occurring in the latter inequality possess the structure as claimed in the statement, we note that occurrences of B(t) in (88) are multiplied by a power of t, so that they do not occur at t = 0. Moreover, the factorΨ 1,t possesses the structure required for Ψ 1,t in the statement. In order to estimate the dependence of q * L ∞ (Qt) on the coefficients c * 1 , we apply the same strategy as in the section 7 of [BDb]: we are done.

7.2.
Estimates for linearised problems for the variables v and ̺. First we state the estimate for the linearised momentum equation. The proof follows the lines of the corresponding result in [BDb]. (Since we can assume ̺ * ∈ [̺ min , ̺ max ], the proof is actually simpler.) The function Ψ 2 is continuous and increasing in both arguments, and it can be chosen such that Ψ 2 (0, a) = (min{1, For the linearised continuity equation, we must acknowledge the main difference with respect to the analysis of the compressible models. Proposition 7.4. Assume that v * ∈ W 2,1 p (Q T ; R 3 ) and that ̺ 0 ∈ W 1,p (Ω) satisfies ̺ min < ̺ 0 (x) < ̺ max in Ω. We define Then the problem ∂ t ̺ + div(̺ v * ) = 0 in Q T with ̺(x, 0) = ̺ 0 (x) in Ω possesses a unique strictly positive solution of class W 1,1 p,∞ (Q T ). Define also M(t) := M(̺, t) (cf. (69)). Then, we can find a constant c depending only on Ω and a function Ψ 3 = Ψ 3 (t, a 1 , a 2 ) continuous and finite in the set {t, a 1 , a 2 ≥ 0 : c a 1 t 1− 1 p a 2 e c t 1− 1 p a 2 < 1} , For i = 3, 4, 5, Ψ i is continuous and increasing in all variables, and Ψ i (0, a 1 , a 2 ) = Ψ 0 i (a 1 ) is independent on the last variable. The identity Ψ 4 (0, a 1 , a 2 ) = a 1 , and the inequality Ψ 5 (0, a 1 , a 2 ) ≤ C a 1 , are also valid.
Proof. The existence statement as well as the construction of the functions Ψ 4 and Ψ 5 is proved in [BDb], Corollary 7.8. The critical point is the construction of the function Ψ 3 . We start from the well-known representation of the solution to the continuity equation (see a. o. [Sol80]) where y(τ ; x, t) is the characteristic curve with speed v * through (x, t). Therefore, Use of |1 − e b | ≤ e |b| |b| allows to bound Owing to the continuity of W 1,p (Ω) ⊂ L ∞ (Ω) and Hölder's inequality .
Thanks to a similar argument applied to ̺ min − ̺, we find that and define the function Ψ 3 to be the right-hand of the latter relation.

The continuity estimate for T
We now want to combine the Propositions 7.2 and 7.3 with the linearisation of the continuity equation in Proposition 7.4 to study the fixed point map T described at the beginning of Section 6 and defined by the equations (46), (47), (48), (49) for given v * ∈ W 2,1 p (Q T ; R 3 ) and q * ∈ W 2,1 p (Q T ; R N −2 ). We define V * (t) := V (t; q * ) + V (t; v * ). At first we state estimates for the lower-order nonlinearities (43), (45).
These estimates were proved in [BDb] for the case that the non-linear coefficients R, M are defined for ̺ * taking values in ]0, +∞[. The proof is exactly the same for ̺ * taking values in I, provided that we adapt the definition of m * (t), M * (t) via (69). Moreover, the arguments are very similar to the ones used to bound the right-hand vector field h. This statement, that we next prove in detail, might serve as an illustration.
We sum up the continuity estimates in the following statement.
Proof. We define a 0 > 0 to be the solution to the equation c M 0 x e c x = 1 associated with the numbers in (91). We apply the Lemma 6.1 with Ψ(t, R 0 , η) := Ψ 6 (t, R 0 , η) from Prop. 8.3, and the claim follows.
9. Proof of the theorem on short-time well-posedness 9.1. Existence and uniqueness. We choose T 0 , η 0 > 0 according to Proposition 8.4. Starting from (q 1 , v 1 ) = 0, we consider a fixed point iteration (q n+1 , v n+1 ) := T (q n , v n ) for n ∈ N. Recalling (67), we define V n+1 (t) := V (t; q n+1 ) + V (t; v n+1 ). Since obviously V 1 (t) ≡ 0, Proposition 8.4 guarantees that From Lemma 9.1 hereafter, we infer that the fixed-point iteration therefore yields strongly convergent subsequences in L 2 (Q T 0 ) for the components of q n , ζ n , ̺ n and v n and for the gradients q n x , ζ n x and v n x . The passage to the limit in the approximation scheme is then a straightforward exercise, since we can rely on a uniform bound in X T 0 . The proofs are almost identical with the fixed-point iteration in [BDb]. We leave the minor changes to the interested reader, and state without proof the following iteration lemma.
Then there are k 0 , p 0 > 0 and 0 < t 1 ≤ T 0 such that for all t ∈ [0, T 0 − t 1 ], the quantity satisfies E n+1 (t) ≤ 1 2 E n (t) for all n ∈ N. 9.2. Verification of continuation criteria. In order to complete the proof of the Theorems 2.1, 2.2 it remains to investigate the claimed characterisations of the maximal existence interval.
Proof. First criterion (1). We must show that the quantity V (t; q) + V (t; v) is bounded by a continuous function of t, M(̺, t), N (t). We will only sketch this point, which relies on going carefully through the proofs of the estimates in the Propositions 7.2, 7.3 in the spirit of [BDb].
To begin with, we notice that the components of v x have all spatial mean-value zero over Ω due to the boundary condition (30). Hence, for For the solution to the continuity equation, Theorem 2 of [Sol80] (see also Proposition 7.7 in [BDb]) implies that sup τ ≤t [̺(τ )] C α (Ω) is bounded by a function of´t 0 [v x (τ )] C α (Ω) dτ , thus also by a function of N (t). Moreover, as in the same references, we show for all t ≥ 0 that Here and throughout the proof, we denote by φ some generic continuous function increasing in its arguments, and R 0 stands for the initial data and the external forces. We next exploit the momentum balance equation for v. We apply Proposition 7.3, hence ). The function f obeys (45) and therefore Coefficients depending on ̺ and q can in general be bounded following the example of Therefore, we show that As shown, ∇̺ p L p (Qt) ≤ φ(R 0 , v x L ∞,1 (Qt) )´t 0 (1 + V (τ ; v)) p dτ , and ζ satisfies the weak Neumann problem (48), hence We define z = 3 p−2 if 3 < p < 5, z > 1 arbitrary if p = 5 and z = 1 if p > 5. Recalling the continuity of the embedding W 1−2/p p ⊂ L 3p/(5−p) + , we show by means of Hölder's inequality Therefore, combining the latter bounds yields We invoke the Gronwall Lemma, hence V p (t; v) ≤ φ(t, M(̺, t), N (t)) ( ∇q p L p (Qt) +A 0 (t)). Since ∇q L p (Qt) is also controlled by a function of t and N (t), so does V p (t; v). It follows that ∇̺ p L p,∞ (Qt) ≤ φ(t, R 0 , N (t)). For β = 1 − 3/p, the Proposition 7.4 yields that ̺ C β,β/2 (Qt) ≤ φ(t, R 0 , N (t)). Recalling that q satisfies (76), we can now finish the proof as in [BDb], Lemma 9.2. Second criterion (2). The more interesting point is to get rid of the dependence on the distance M(̺, t) to the density thresholds in the estimates. First we note that the relation (28) implies for the gradient of the pressure Clearly, F L p (Qt) is bounded by a function of b and the norms of ζ and v occurring in the quantity K(t). We notice in particular that this function is independent on M(̺, t). In order to obtain a bound on the entire pressure gradient, we employ the continuity equation (27). We compute that Define m(̺, t) := min Qt {1 − ̺/̺ max , ̺/̺ min − 1}. Thanks to Lemma B.1, the properties of the pressure function guarantee that |P ̺ (̺, q)| ≤ c 4 m(̺, t) −1 in Q t . Since |P ̺ (̺, q)| |∇̺| = |∇P (̺, q) − P q (̺, q) ∇q|, the same Lemma B.1 also implies that c 3 m(̺, t) −1 |∇̺| ≤ |∇P (̺, q)| + c 5 |∇q|) .
By means of (98), (95), (96) we see that also ∂ t P (̺, q) L p,p/2 (Q) is bounded by a function of t and K(t), independently of M(̺, t). Overall we have P x L p (Qt) + P t L p,p/2 (Qt) ≤ Ψ. For p > 5, we can show that this implies a bound P L ∞ (Qt) ≤ C(t) Ψ, where C(t) is the embedding constant of an anisotropic Sobolev space into L ∞ (Q t ). It remains to recall that for the choice (7), the function P satisfies (cf. [Dru19], Proposition 5.3) This implies that M(̺, t) ≤ C 1 e C 2 ( P (̺, q) L ∞ (Q t ) + q L ∞ (Q t ) ) , and the claim follows.
10. Global well-posedness 10.1. The map T 1 is well defined. We consider the equations (59), (61), (62), (63), (64) characteristic of the definition of the map T 1 . We recall that these equations are obtained by comparing a solution to some suitable extension (q 0 ,v 0 ) ∈ Y T , to be constructed here below, of the initial data. The initial density ̺ 0 is extended by a function̺ 0 obtained via the solution of (50). We moreover introduce the functionζ 0 , solution to (51). In order to define T 1 we must make sense of the linear operators (g 1 ) ′ (u * ,û 0 ), (h 1 ) ′ (u * ,û 0 ) and (f 1 ) ′ (u * ,û 0 ). The density components in the vectorsû 0 = (q 0 , 1,̺ 0 ,v 0 ) and u * (def. in (60)) must therefore assume values in I up to time T > 0! This property is to be expected if the initial data are close enough to an equilibrium solution (ρ eq , p eq , v eq ) defined by the relations (16), (17). The distance of the initial data to this solution is expressed by the number in which ̺ eq := N i=1 ρ eq i and q eq ℓ = η ℓ · ∇ ρ k(ρ eq ) for ℓ = 1, . . . , N − 2. Throughout this section, we moreover employ the abbreviation Observe the occurrence of the higher-order W 2,0 p -norm of̺ 0 in the definition of R 0 . To commence with, we recall a result of [BDb] for estimating the gradient of solutions to a perturbed continuity equation. The proof in [BDb] is given for zero initial conditions, but the extension to the nonzero case is completely straightforward.
10.2. Continuity estimates. We need at first an estimate for the operators (g 1 ) ′ , (h 1 ) ′ and (f 1 ) ′ . We shall prove it for general body forces b = b(x, t), even if the statement of Theorem 2.3 requires only b = b(x).
Next we prove the main continuity estimate. We apply Proposition 7.2 to (61), (62). Making use of the fact that r(0, x) = 0 in Ω, we get the estimate ) . HereΨ 1,T = max{Ψ 1,T , Φ T } depends continuously on the data. We then apply Proposition 7.3 to (64) and obtain again with someΨ 2,T depending on T and sup s≤t [̺ * (s)] C α (Ω) ). We estimate ∇χ L p (Qt) by means of (113). We next raise both (113) and (114) to the p th power, add both inequalities, and get for the function V (t) := V (t; r) + V (t; w) + χ W 2,0 p (Qt) the inequality ) . Then we make use of w p W 1,0 p (Qt) ≤´t 0 V p (s) ds, and we apply the Lemma 10.2 to find that The Gronwall inequality implies that ) . We thus have proved the following continuity estimate.
Proposition 10.3. We define R 0 via (100). Suppose that (r * , w * ) ∈ 0 Y T satisfy the condition (108). Then (r, w) = T 1 (r * , w * ) is well defined in 0 Y T . Moreover, there is a continuous function Ψ 7 = Ψ 7 (T, R 0 , η), increasing of all arguments and finite for all ) . 10.3. Existence of a unique fixed-point of T 1 . We are now in the position to prove a self-mapping property for sufficiently 'small data'. We recall the definitions (99), (100) of the critical norms R 0 , R 1 . We denote u eq = (q eq , ζ eq , ̺ eq , v eq ) and letû 0 := (q 0 ,ζ 0 ,̺ 0 ,v 0 ) andû 1 := u eq −û 0 . In (104), (107), we just proved that û 1 X T ≤ C R 1 . Recalling that the operator A is continuously differentiable into the space Z T defined in (38), and that A(u eq ) = 0 by the definition of an equilibrium solution, we can verify that These considerations allow to state and prove the main properties of T 1 .
Proof. If w * satisfies (108) and if R 1 satisfies (105), T 1 (r * , w * ) is well defined in Y T . We apply Proposition 10.3, use (116) and obtain We consider the iterationū n+1 := T 1 (ū n ), starting at zero. The sequence (q n , ζ n , ̺ n , v n ) is then uniformly bounded in X T . We show the contraction property with respect to the same lower-order norm than in Lemma 9.1. There are k 0 , p 0 > 0 such that the quantities E n (t) := p 0ˆt +t 1 t {|∇(r n − r n−1 )| 2 + |∇(χ n − χ n−1 )| 2 + |∇(w n − w n−1 )| 2 } dxds Appendix A. Properties of the free energy density In this section we prove the statements of the section 3 devoted to the convex conjugate of the free energy density: Lemma 3.1, Lemma 3.2 and Lemma 3.3. We assume that k satisfies the assumptions of Lemma 3.1. Notice that requiring k essentially smooth on S 1 , while positive homogeneous, induces that k is also essentially smooth on SV . To see this, we consider any sequence {r m } ⊂ SV such that r m →r for m → ∞, and r belongs to the relative boundary of SV , which means that there is i ∈ {1, . . . , N} such thatr i = 0. Then we define y m :=r m / N i=1 r m i which belongs to S 1 for all m, and satisfies y m i → 0 for m → ∞. Since k is positively homogeneous, we have ∇ ρ k(r m ) = ∇ ρ k(y m ). Thus, by the assumptions of Lemma 3.1, we see that |∇ ρ k(r m )| → +∞, which is the essential smoothness on SV . Consider now µ ∈ R N arbitrary. Then we claim first that there exists a uniquer ∈ SV such that Since SV is bounded, we first notice that sup r∈SV {µ · r − k(r)} = max r∈SV {µ · r − k(r)}. Thus, there isr ∈ SV such that sup r∈SV {µ · r − k(r)} = µ ·r − k(r). We want to show that r is an interior point. Since SV is a convex set, we can find for every a ∈ SV a h > 0 such thatr + h (a −r) ∈ SV . Due to the choice ofr which yields k(r + h (a −r)) − k(r) ≥ h µ · (a −r) and lim hց0 k(r+h (a−r))−k(r) h > −∞. The latter however contradicts the fact that k is essentially smooth on SV (cf. [Roc70], Lemma 26.2). Thus,r ∈ SV is an interior point. The uniqueness ofr follows from the strict convexity of k on SV . Since k is differentiable, and since r → µ · r − k(r) attains its maximum inr, we must have (∇k(r) − µ) · ξ = 0 for every tangential vector ξ ∈ R N such that ξ ·V = 0. Thus, there is p ∈ R such that µ = ∇ ρ k(r) + pV . Multiplying withr, use of the homogeneity of degree one implies thatr · ∇ ρ k(r) = k(r), hence sup r∈SV {µ · r − k(r)} = µ ·r − k(r) = pr ·V = p , showing that p = f (µ). Due to the structure f (µ) = µ ·r − k(r) = max r∈SV {µ · r − k(r)}, we easily show that f is differentiable in µ with ∇ µ f (µ) =r. In order to show the differentiability of higher order, we can exploit the identities For a system of orthonormal vectors ξ 1 , . . . , ξ N −1 for {V } ⊥ , and ξ N :=V /|V |, we then have The latter can be viewed as an algebraic system of the form F (X) = (µ·ξ 1 , . . . , µ·ξ N −1 , 1 |V | ) for the unknowns X := (ξ 1 · ∇ µ f (µ), . . . , ξ N · ∇ µ f (µ)) ∈ R N . The Jacobian of this system obeys where D 2 k is evaluated at ∇ µ f (µ). We can easily verify that {D 2 k(∇ µ f (µ))ξ i ·ξ j } i,j=1,...,N −1 is strictly positive definite: A vector of the form N −1 j=1 ξ j a j , a = 0 can never be parallel to ∇ µ f (µ), since multiplying withV yields a contradiction. On the other hand, the properties of k guarantee that the kernel of D 2 k(∇ µ f (µ)) is the one-dimensional span of ∇ µ f (µ). Thus, the equations F (X(µ)) = (µ · ξ 1 , . . . , µ · ξ N −1 , 1 |V | ) define implicitly a map µ → X(µ) of class C 1 (R N ). This clearly implies that f ∈ C 2 (R N ), and we obtain the formula If k ∈ C 3 (R N + ), we then differentiate again to obtain that f is C 3 (R N ). This proves the claims of Lemma 3.1. The claims of Lemma 3.2 and 3.3 are also readily established (use (118) and (117)).

Appendix B. Auxiliary statements
For the proof of the following Lemma, we need the variable transformation in the section 4.1.
Moreover, by the same means Arguing the same for the other choice of u, it follows that |∂ q ℓ M i,j (R(̺, q))V j | ≤ C m(̺).
The other estimates claimed have been verified for this special case of the function k in the Section 4 of [Dru19].
Remark B.2. In the case that the matrix M results from inversion of the Maxwell-Stefan equations, we notice that the matrix B of Lemma B.1 is nothing else but the pseudo-inverse of the Maxwell-Stefan matrix. It is shown in the paper [BDc] that natural assumptions on the binary diffusivities are sufficient for proving that the entries of B consist of regular functions of the state variables. In particular, they satisfy the assumptions of Lemma B.1.
We also recall some estimates of Hölder norms. This is also proved in [BDb].
Finally we have a perturbation Lemma for elliptic problems. This property ought to be well known, and we only mention details for more convenience on reading.
Here m 0 is some geometric constant associated with the covering of Ω. It remains to estimate where we employ the preliminary consideration at the beginning of this proof to show that u L ∞ ≤ c F L p . Recalling the choice of r, we are done.