Existence analysis of a stationary compressible fluid model for heat-conducting and chemically reacting mixtures

The existence of large-data weak solutions to a steady compressible Navier-Stokes-Fourier system for chemically reacting fluid mixtures is proved. General free energies are considered satisfying some structural assumptions, with a pressure containing a $\gamma$-power law. The model is thermodynamically consistent and contains the Maxwell-Stefan cross-diffusion equations in the Fick-Onsager form as a special case. Compared to previous works, a very general model class is analyzed, including cross-diffusion effects, temperature gradients, compressible fluids, and different molar masses. A priori estimates are derived from the entropy balance and the total energy balance. The compactness for the total mass density follows from an estimate for the pressure in $L^p$ with $p>1$, the effective viscous flux identity, and uniform bounds related to Feireisl's oscillations defect measure. These bounds rely heavily on the convexity of the free energy and the strong convergence of the relative chemical potentials.


Introduction
This paper is concerned with the mathematical analysis of a stationary model for fluid mixtures, coupling Maxwell-Stefan diffusion with a Navier-Stokes-Fourier model. A key feature of the present paper is that the fluid mixture is modeled in a thermodynamically consistent way. Compared to [6], we include temperature effects, and compared to [31,32], our constitutive equations are different and we allow for temperature gradients inside the diffusive fluxes, which yields a cross-diffusion coupling between the equations for the partial mass densities and the equation for the energy. The most important issue which allows for better results than in previous papers for steady compressible models of chemically reacting mixtures is the convexity of the Helmholtz free energy, similarly as in [6]- [8], where an evolutionary model was studied, however, under the assumption that the temperature is constant. On the other hand, we neglect chemical reactions on the boundary which is an important issue in the aforementioned papers.

Variable
Physical meaning ρ i partial mass density of the ith species ρ = N i=1 ρ i total mass density m i molar mass of the ith species n i number density, n i = ρ i /m i n total number density, n = N i=1 n i J i partial flux of the ith species v barycentric velocity r i reaction term for the ith species p pressure S viscous stress tensor b momentum force term J entropy flux Ξ entropy production Q internal heat flux µ i chemical potential of the ith species θ thermodynamic temperature s specific entropy e specific internal energy E = e + |v| 2 /2 specific total energy ψ specific Helmholtz free energy Table 1. Physical meaning of the variables.
We also prescribe the total mass where |Ω| denotes the measure of Ω, and ρ is assumed to be positive. Note that the total mass is in fact |Ω|ρ and the quantity ρ has the meaning of the average density. In some situations, especially when the integrability of the density is low, we replace the total energy balance by the entropy inequality (1.8) − div J + Ξ ≤ 0.
Note that we replaced one equality by one inequality. Hence, we may obtain too many solutions to our problem which are non-physical. However, due to mathematical reasons (and physically, it is not surprizing either), as explained below, we cannot expect equality for the balance of entropy. To avoid this problem, similarly as for the single-component steady Navier-Stokes-Fourier system (see, e.g., [18]), we also add the total energy balance (1.3) integrated over Ω, We will discuss this issue later; at this point, we just note that (1. 1.2. Notation. The unit matrix in R m×m is denoted by I, where m ∈ N. Given two matrices A, B ∈ R m×n , we define A : B ≡ m i=1 n j=1 A ij B ij . We denote by a bold letter a vector a ∈ R 3 with the components (a 1 , a 2 , a 3 ), by a blackboard bold letter a matrix A ∈ R 3×3 with the coefficients (A ij ) or (A ij ), and by a the vector (a 1 , . . . , a N ) in R N . Recall that N denotes the number of species in equation (1.1). We also set R + = [0, ∞).

Constitutive equations.
We specify the entropy flux J , entropy production Ξ, viscous stress tensor S, heat flux Q, diffusion fluxes J i , and reaction terms r i .
The entropy flux J is the sum of the contributions of free transport, diffusion fluxes, and heat flux, The entropy production Ξ keeps into account the contributions from the diffusion fluxes, the heat flux, the viscous stress and the reaction terms, According to the second law of thermodynamics [21], the entropy production Ξ must be nonnegative. This is achieved by requiring that The following definitions of S, Q, J i , and r i are chosen in such a way that these requirements are satisfied.
The viscous stress tensor is assumed to be a linear function of the symmetric velocity gradient, where λ 1 (θ) and λ 2 (θ) are the temperature-dependent shear and bulk viscosity coefficients, respectively.
The heat flux consists of the Fourier law and the molecular diffusion term, where κ(θ) is the thermal conductivity and the coefficients M i depend on ρ and θ. The molecular diffusion term plays an important role in the energy identity.
To fulfill the mass conservation equation div(ρv) = 0, the sum of the diffusion fluxes and the sum of the reaction terms should vanish, i.e. N i=1 J i = 0, N i=1 r i = 0. Hence, the diffusion matrix has a nontrivial kernel, and we assume that (1.14) To be consistent with the first inequality in (1.10) and relations (1.12) and (1.13), we assume that κ(θ) > 0 and that the mobility matrix (M ij ) is nonnegative. More precisely, (1.15) ∃C M > 0 : where Π = I − 1 ⊗ 1/N is the orthogonal projector on span{ 1} ⊥ . For the structure of the mobility matrix and its connection to the Maxwell-Stefan theory, we refer to [1]. According to the third inequality in (1.10), the entropy production due to reactions, − N i=1 r i µ i /θ, must be nonnegative. Therefore, we suppose for the reaction terms that (1.16)      r i = r i (Π( µ/θ), θ), i = 1, . . . , N, ∃C r > 0, ζ > 0, β > 0 : − N i=1 r i (Π( q), θ)q i ≥ C r |Π q| 2 , |r i (Π( q), θ)| ≤ C r (|Π q| 5(6−ζ)/6 + θ 5(3β−ζ)/6 ) for all q ∈ R N , θ > 0, for some ζ > 0 possibly very small, β > 0, and all i = {1, 2, . . . , N}. This condition is satisfied for the reaction terms used in [6], (1.17) r γ j k µ k θ , i = 1, . . . , N, j = 1, . . . , l, where Φ: R l → R is a convex potential with suitable growth and γ j ∈ R N is a vectorial coefficient associated to the jth reaction; see Remark 1.3. The remaining variables-chemical potentials µ i , pressure p, total internal energy e, and entropy s-are determined from the Helmholtz free energy density ρψ which is assumed to be a function of ρ and θ: The definition of the pressure is known as the Gibbs-Duhem relation. In this paper, we allow for general free energies satisfying Hypothesis (H6) below. A specific example is given in Remark 1.2.
(H7) Pressure: p is defined by (1.18) and satisfies Remark 1.1 (Discussion of the hypotheses). (H1) If we approximate Ω by smooth domains, it is possible to extend our existence result to less regular domains as, e.g., Lipschitz ones. As this would technically complicate the paper, we skip this and only point out the paper [27], where such a technique is used in the case of compressible Navier-Stokes equations. (H2) The force term b is assumed to be dependent only on x just for simplicity. It may also depend on ( ρ, θ); then b(x, ρ, θ) needs to be measurable in the first variable, continuous in the last N + 1 variables, and bounded. We may assume that α 1 = 0 (but not α 2 = 0); then we need that Ω is not axially symmetric to apply the Korn inequality [28, (4.17.19)]. In case α 1 > 0, we have the Korn inequality also for axially symmetric domains, but the estimate of the velocity gradient depends additionally on the tangential trace of the velocity. However, as shown below in Lemma 2.2, the estimate of the trace of v depends also on the density and therefore, the problem becomes slightly more complex and leads to additional restrictions on the exponents β and γ (cf. [18]). To avoid such problems, we prefer to assume that the domain is not axially symmetric. (H3) The linear growth of the viscosities leads to optimal a priori estimates. The lower bound avoids degeneracies in the coefficients. (H4) The growth assumption on M ij , M i is necessary to have strong relative compactness M ij ( ρ δ , θ δ ), M i ( ρ δ , θ δ ) in L 3 (Ω) (for a suitable subsequence of ( ρ δ , θ δ )). (H5) We show in Remark 1.3 that the reaction terms (1.17) satisfy condition (1.16).
This condition excludes vanishing reaction terms. In fact, it is needed to derive an L 2 (Ω) bound for Π q and, together with (1.15), an H 1 (Ω) bound for Π q (see Lemma 2.2). Condition (1.16) can be replaced by a Robin-type boundary condition for the diffusion fluxes involving the chemical potentials, as done in [5,Hypothesis (H8)]. (H6) The fact that ρ → h θ ( ρ) is strictly convex implies that ∇ ρ h θ ( ρ) is a strictly monotone (and therefore invertible) operator from R N + to R N . We do not impose any assumptions directly on the growth of ∇ ρ h θ ( ρ), since we require some growth conditions on the pressure stated in Hypothesis (H7). Condition (1.22), fulfilled by our example of a free energy from Remark 1.2 presented below, is rather technical and is used only in the construction of a solution, not in the a priori estimates or compactness part. Thus, it can be viewed as a technical assumption, and we expect that it can be avoided by using a better approximation scheme. (H7) We present below an example of the free energy for which p( ρ, θ) = nθ + (γ − 1)n γ .
This pressure satisfies Hypothesis (H7). The condition γ > 3/2 is needed to derive a bound of ρ in L γ+ν (Ω); see Lemma 2.3. This allows us to control the pressure in a better space than L 1 (Ω); see Lemma 2.4. Remark 1.2 (Example of a free energy). An example of the free energy density fulfilling Hypothesis (H6) is given by where γ > 1 is some exponent and c W > 0 is the heat capacity; see Appendix B for the proof. The first term is related to the mixing of the components, the second one is needed for the mathematical analysis (to obtain an estimate for the total mass density; a certain physical justification of this term can be found in [12,Chapter 1]), and the third one describes the thermal energy. We have assumed that the specific volumes of all components are the same. In the paper [6,Formula (23)], the first term is defined in a different way: This expression does not contribute to the pressure, while our mixing entropy term gives the pressure contribution nθ of an ideal gas. With the free energy density (1.25), we obtain n i log n i + c W ρ(log θ + 1).
The definition of r i differs slightly from that one in [6] because of the role played by the temperature. In fact, the entropy inequality provides us with an estimate for γ j · µ/θ which in turn yields an L 2 (Ω) bound for Π( µ/θ) (see Lemma 2.2). We require that N i=1 γ j i = 0 for all j = 1, . . . , l in order to achieve N i=1 r i = 0, needed for mass conservation. As a consequence, γ j ∈ span{ 1} ⊥ . We assume that the linear hull of all γ j is in fact equal to span{ 1} ⊥ , which implies that l ≥ N − 1. This condition is necessary since Theorem 8.3 in [7], which gives an L 2 (Ω) estimate for Π( µ/θ), cannot be generalized in a straightforward way to the non-isothermal case.
For all p ∈ [1, ∞], we introduce the following spaces: For a certain choice of parameters, we always obtain the latter while the former will be satisfied only if some terms from the total energy balance are integrable, hence for a smaller set of the parameters. Let us explain the definition of our solutions more precisely. Before doing so, we detail the weak formulation of our equations. We consider • the weak formulation of the mass balance, for all φ 1 , . . . , φ N ∈ W 1,∞ (Ω); • the weak formulation of the momentum balance, • the weak formulation of the total energy balance, for all ϕ ∈ W 1,∞ (Ω); and • the weak formulation of the entropy inequality, (1.29) Φ ds, in Ω. We also use the term the global total energy equality, which is nothing but equality (1.28) with ψ ≡ 1, i.e.
Finally, recall that if we choose φ i ≡ 1 in (1.26) and add the weak formulation for all constituents, we obtain in the sum the weak formulation of the continuity equation, for all Φ ∈ W 1,∞ (Ω). However, we shall work with a slightly stronger definition of a solution to this equation, with the renormalized solution to the continuity equation, which is an important tool in the theory of weak solutions to the compressible Navier-Stokes equations to show compactness of the sequence of densities. Hence, we consider only renormalized solutions in what follows, i.e., solutions satisfying in addition to the weak formulation (1.31) for u ∈ W 1,∞ (Ω) and b ∈ C 0,1 (R) with b ′ having compact support, Definition 1.4 (Weak and variational entropy solutions). We call the functions  Remark 1.6. Replacing the estimates of the total density computed from the Bogovskii operator estimates by a different technique used, e.g., in [30] (see also [18,26,31,32]), we could treat the problem with arbitrary γ > 1 and certain bounds on β depending on γ. However, it would significantly complicate and extend this quite technical and long paper. Therefore we prefer not to do it here. Remark 1.7. Note that in (1.26) we may use test functions φ i = 0 and φ j being a non-zero constant for all j = i. This leads to a sort of compatibility condition, Ω r j dx = 0, j = 1, 2, . . . , N.
We do not require this condition directly, but due to our assumptions on r j and due to the construction, we are able to fulfil these conditions; see Remark 5.3.
1.6. Key ideas of the proof. We first prove the existence of solutions to an approximate Navier-Stokes-Fourier system, then derive a priori estimates from the entropy and total energy balances, and finally pass to the limit of vanishing approximation parameters. The approximate system is obtained from the mass densities and momentum balance equations We use several levels of approximation. They are described in detail in Section 4 dealing with existence of solutions. In the construction of a solution for the approximate problem, we use for the velocity a Galerkin approximation with dimension n ∈ N; we add lower-order regularizations with parameter ε > 0 and higher-order regularizations with parameters (χ, δ, ξ) to the other equations, leading to H 2 (Ω) regularity; and we regularize the free energy with parameter η > 0 to obtain higher integrability of the total mass density. Moreover, we construct not the temperature, but rather its logarithm which allows us to deduce the a.e. positivity of the temperature.
We cannot use the artificial viscosity regularization in the total mass continuity equation as in [12, Section 3.3.1], yielding a linear elliptic equation with drift, because of the cross-diffusion terms, which significantly complicates our analysis. In fact, we need two approximation levels, distinguishing between higher-order and lower-order regularizations.
The existence of a solution to the approximate scheme is proved in Section 4. The positivity of the temperature is obtained from the bound log θ ∈ W 1,4 (Ω) ֒→ L ∞ (Ω). The uniform estimates from the entropy and total energy balance allow us to pass to the limits δ → 0 and n → ∞ (in this order). For the limits χ → 0, ε → 0, and ξ → 0 (in this order), we need an estimate for the total mass density in L γ+ν (Ω) for some ν > 0. It yields a uniform bound for the pressure in a space better than L 1 (Ω). This is shown by using a test function involving the Bogovskii operator in the weak formulation of the momentum balance equation (see [28,Section 7.3.3]).
The most difficult part of the proof is the strong convergence of the approximate mass densities (ρ ε ). It is based on an effective viscous flux identity or weak compactness identity for the pressure [12, Section 3.7.4] (Lemma 3.2) and some properties related to Feireisl's oscillations defect measure. More precisely, we shall prove that converges strongly in L 1 (Ω) to some function, where the "overline" signifies the weak limit and T k is some truncation operator (Lemma 3.3), and that ) converges strongly in L 1 (Ω) to zero (Lemma 3.5). The proofs of the strong limits of (1.33) and (1.34) rely heavily on the convexity of the free energy density and the strong convergence of the relative chemical potentials. These limits are necessary to deal with the cross-diffusion terms, and the proofs seem to be new. Another ingredient is the proof that the weak limit ρ of (ρ ε ) is a renormalized solution to the mass continuity equation. Using a special test function and renormalization function, we are able to control the oscillations defect measure and to prove the strong convergence of (ρ ε ) to ρ.
We finally underline that the method used in this paper is essentially new in the aspect how we treat the strong convergence of the total density. Note that this allows us to prove the result under minimal conditions and in fact also to improve the known results for the steady compressible Navier-Stokes-Fourier equations. 1.7. State of the art and originality of the paper. The literature on compressible Navier-Stokes and Navier-Stokes-Fourier systems is very extensive. First results on the existence of solutions to the stationary compressible Navier-Stokes equations (for a single species) without assumptions on the size of the data goes back to P.-L. Lions [22]. He assumed the pressure relation p(ρ) = ρ γ with γ > 5/3 for the stationary flow. The most difficult part of the proof, the strong convergence of the sequence of approximate densities, is based on a weak continuity property of the effective viscous flux p(ρ) + (2λ 1 + λ 2 ) div v and the theory of renormalized solutions applied to the continuity equation. A first improvement on the pressure exponent γ is due to Březina and Novotný [3], who assumed that γ > (1 + √ 13)/3. Their proof is based on some ideas due to Plotnikov and Sokolowski [33] and on Feireisl's idea of the oscillations defect measure estimate (see [9] in the evolutionary case), described in the steady case in [28] for a special class of non-volume forces. Further improvements of the lower bound for the exponent γ are due to Frehse, Steinhauer, and Weigant [15]. The existence of weak solutions to the steady compressible Navier-Stokes equations for any γ > 1 was shown in [20] (for space periodic boundary conditions), in [19] (for slip boundary conditions), and in [34] (for Dirichlet boundary conditions). In the case of evolutionary Navier-Stokes equations, the existence of a solution was proved in [22] for γ ≥ 9/5 and the lower bound was decreased to γ > 3/2 in [13].
The existence theory for the Navier-Stokes-Fourier system employs the techniques of both Lions and Feireisl. The first result for the stationary compressible Navier-Stokes-Fourier system was proved by P.-L. Lions [22] under the assumption that the density is bounded in some L q space for sufficiently large values of q. Without this assumption, when the density is a priori controlled only in L 1 (Ω), the existence of weak solutions was shown in [24] for γ > 3 and in [25] for γ > 7/3. These results were improved in a series of papers, see [29,30] for Dirichlet boundary conditions and [18] for the Navier boundary conditions, showing the existence of a variational entropy solution (satisfying the entropy inequality and the global total energy balance) for any γ > 1 and the existence of a weak solution (satisfying the total energy balance) for any γ > 5/4 (Navier boundary conditions) or γ > 4/3 (Dirichlet boundary conditions). We refer to [26] for further information.
For results on the time-dependent compressible Navier-Stokes and Navier-Stokes-Fourier equations, we refer to the monographs [11,12,28]. One difficulty is the proof of the strong convergence of the sequence of approximate temperatures which makes necessary the application of the div-curl lemma and the effective viscous flux relation by using a cancellation property different from the isothermal model. The transient Navier-Stokes equations with density-dependent viscosities satisfying a certain structure condition allow for new a priori estimates thanks to the Bresch-Desjardin entropy [2]. However, these estimates are not available for the steady problem. The evolutionary problem for a heat conducting fluid with basically the same restriction on the adiabatic coefficient γ was considered in [10,12] for different formulations of the energy balance (internal energy inequality and entropy inequality, respectively).
All these results concern the single-species case. The theory of fluid mixtures requires some careful modeling to maintain thermodynamic consistency; we refer to [1,4,16] for the thermodynamic theory of fluid mixtures. One of the first results was proved in [36], namely the existence of weak solutions to the stationary Navier-Stokes equations assuming Fick's law and γ > 7/3. This result was improved in [17] for Maxwell-Stefan-type fluxes and γ > 5/3. Another improvement was achieved in [31,32] for variational entropy solutions with γ > 1 and for weak solutions with γ > 4/3. These results are based on the assumption of same molar masses for each constituent. Concerning the evolutionary problem, the first global in time result for arbitrarily large data is due to [14] for Fick's law. An existence result for a general thermodynamically consistent transient Navier-Stokes model with γ > 3/2 was recently presented in [6]. Electrically charged dynamical incompressible mixtures were analyzed in [5]. We also mention the work [23] in which a multicomponent viscous compressible fluid model with separate velocities of the components was studied.
Our existence result generalizes previous works on fluid mixtures. Indeed, the mobility matrix in [17,31,32] has no contributions to the temperature gradients, the pressure and internal energy are not defined through the free energy, and the molar masses are assumed to be equal. These restrictions were removed in [5], but the authors consider incompressible fluids. In the works [6]- [8], a very general compressible fluid model is analyzed but no temperature effects have been taken into account.
In this paper, we combine all the features studied in the above-mentioned works, namely we allow for temperature gradients, a thermodynamically consistent modeling starting from the Helmholtz free energy, compressible fluids, and different molar masses. Moreover, we obtain a new proof for the strong convergence of the sequence of approximate densities by exploiting the convexity of the free energy.
The paper is organized as follows. A priori estimates for smooth solutions are derived in Section 2. We prove in Section 3 the compactness of the sequence of total mass densities satisfying the Navier-Stokes-Fourier mixture model. This step highlights the key features and novelties of the proof without obstructing it by the numerous approximating terms. The construction of smooth approximate solutions and the deregularization limits is presented in Sections 4 and 5, respectively. In the Appendix, we recall some auxiliary results needed in this paper and we show that the free energy density in Remark 1.2 satisfies Hypothesis (H6).

A priori estimates for smooth solutions
This section is devoted to the derivation of suitable a priori estimates for smooth solutions. Although we consider later weak solutions only, the computations help us to identify the key estimates of the existence proof. In fact, we need several regularizations for the full model which may abstruse the main arguments.
We assume the existence of a smooth solution to the following Navier-Stokes-Fourier system with chemically reacting species: Note that on this level, we can freely switch from the internal energy balance (2.3) to the total energy balance or the entropy equality; see Lemma 2.1 below. In order to use the procedure also during the construction of the solution, we immediately obtain the integrated entropy and total energy balances. To obtain the weak formulation of the entropy and total energy balances, we can proceed as in the proof of the lemma below, we just multiply the corresponding strong formulation additionally by a smooth test function ψ, which gives several additional terms containing the derivatives of ψ. In case of the entropy (in)equality, later on, we also require that the test function is nonnegative.
Lemma 2.1 (Entropy and total energy balances). Let ( ρ, v, θ) be a smooth solution to (2.1)-(2.6) such that θ > 0. Then the entropy balance and the total energy balance Proof. We multiply (2.1) by µ i /θ and (2.3) by −1/θ, sum both equations, sum over i = 1, . . . , N, integrate over Ω, and integrate by parts (up to one term). The terms involving the coefficients M i cancel, and, taking into account the boundary conditions (2.4)-(2.6), it follows that We claim that some of the terms cancel, namely To prove this, we use the thermodynamic relations (1.18) to deduce that (we denote H( ρ, θ) = ρψ( ρ, θ) for a moment) Hence, we deduce (2.10) after an integration by parts, and (2.9) simplifies to (2.7). Next, we multiply (2.2) by v and add to the resulting equation the energy balance (2.3), integrate over Ω, and integrate by parts. The integrals involving S : ∇v and p cancel and we end up with The first term vanishes. Indeed, the sum of (2.1) from i = 1, . . . , N yields div(ρv) = 0. Multiplying this equation by |v| 2 /2 and integrating over Ω, an integration by parts gives By (2.5), it yields the second identity (2.8), finishing the proof of the lemma.
The entropy and total energy balances yield some a priori estimates. We define Lemma 2.2 (Estimates from the entropy balance). The following a priori estimates hold: where here and in the following, C > 0 denotes a generic constant dependent only on the given data.
Proof. We claim that every term on the left-hand side of (2.7) is nonnegative. In view of (1.15), we need to consider only the last two terms. We deduce from Hypothesis (H3) and the Korn inequality (Lemma A.3 in Appendix A), taking into account that Ω is not axially symmetric thanks to Hypothesis (H1), that for all v ∈ H 1 ν (Ω; R 3 ), The L 2 (Ω) bound for ∇Π( q) is a consequence of (1.15), and (1.16) gives an L 2 (Ω) bound for Π( q) i , i = 1, . . . , N. At this point, we need the nonvanishing reaction terms. Thus, Π( q) is bounded in H 1 (Ω).

Lemma 2.3 (Estimate for the total mass density). Let
Then there exists a constant C > 0 depending only on the given data such that Proof. The proof is based on estimates from the momentum balance (2.2) using the Bogovskii operator B. We refer to Theorem A.1 in Appendix A for some properties of this operator. We multiply Recall that, by Hypothesis (H7), Ω pρ ν dx ≥ c p Ω (ρ γ+ν + ρ 1+ν θ) dx. We estimate the right-hand side of (2.14) term by term. We start with the two delicate terms which lead to the restrictions on the exponent ν. We have for α > 3/2, Next, in view of Hypothesis (H3), As we control the L 1 (Ω) norm of the density (see (1.7)), we can control Ω ρ γ dx Ω ρ ν dx by C ρ λ L γ+ν (Ω) for some λ < γ + ν, by interpolating between the L 1 and L γ+ν norms. Hence, we only need to deal with the part of the first integral containing the temperature.
Let us first consider the case ν ≤ 1. Then Ω ρ ν dx is bounded by a constant and The first term on the right-hand side is bounded, and the last term can be absorbed by the left-hand side of (2.14) for sufficiently large K. Next, note that for ν > 1 we have 2γ − 3 > 1, i.e. γ > 2. Then, by Hölder's inequality, It follows from γ > 2 and β > 2/3 that γ/(γ − 1) < 3β. Hence, using once more γ > 2, Collecting all estimates, we deduce from 2 + ν < γ + ν that where λ < γ + ν. This leads to the desired estimate of ρ.

Weak sequential compactness for smooth solutions
In this section, we focus on the weak sequential stability of a weak solution and formulate it as an independent result. Then, in Section 4, we will just adapt the method introduced here and use it to prove the existence of a weak solution. The main result of this part is the following theorem.
We shall prove the theorem in several steps and each step is described in one of the following subsections.

3.1.
Convergence results based on a priori estimates. Based on Lemmata 2.2-2.4, we collect all weak convergence results that allow us to pass to the limit in the weak formulation. It remains to show that the densities ρ δ converge pointwise to identify the pressure and the chemical potentials as functions of the densities and the temperature.
Limit in the mass balance. First, using the facts that (up to subsequences) ρ δ → ρ weakly in L γ+ν (Ω; R N ), θ δ → θ strongly in L r (Ω), 1 ≤ r < 3β, together with Hypothesis (H4), we see that where i, j = 1, 2, . . . , N as δ → 0 and a bar over a quantity denotes its weak limit. Since the partial densities converge only weakly, we cannot generally identify the weak limits with M ij ( ρ, θ) and M i ( ρ, θ), respectively. Furthermore, recalling that we deduce from Hypothesis (H5) that We do not know at this moment whether q = µ/θ, where µ is given by (1.18). Therefore, letting δ → 0 in the weak formulation (1.26) of the mass balance, we infer from Hypotheses (H4) and (H5) that, for all φ 1 ,. . . , Limit in the momentum balance. By Lemma 2.4, At this point, it is not clear whether p( ρ δ , θ δ ) =: p = p(ρ, θ) and this will be proved later. The weak convergence of (a subsequence of) v δ in H 1 (Ω) and the strong convergence of (θ δ ) in L r (Ω), r ≥ 2, implies that where S δ and S are the stress tensors (1.11) associated to θ δ , v δ and θ, v, respectively. We use Hypothesis (H3) to control the viscosities. Therefore, we can perform the limit δ → 0 in the weak formulation (1.27) of the momentum balance; hence for all u ∈ W 1,∞ ν (Ω), Limit in the entropy inequality. In view of (1.18), Hypothesis (H6) (in particular (1.23) and (1.24)), and the bounds on the temperature and densities, we obtain Using the weak lower semicontinuity in several terms, the previous weak limits, Hypotheses (H3)-(H6), and Lemma A.4, we conclude from (1.29) in the limit δ → 0 that Limit in the total energy balance. The problem with the total energy balance is more complex. We can easily pass to the limit if the test function is constant, yielding the global total energy equality (1.30). To obtain a suitable limit in the weak formulation (1.28) of the total energy balance, we have to assume that γ > 5/3 and β > 1. This ensures that γ + ν > 2 and for some r > 1. Moreover, in view of Hypothesis (H6), for some r > 1. Therefore, letting δ → 0 in (1.28), it follows for all ϕ ∈ W 1,∞ (Ω) that It remains to verify that in Ω as well as to identify Π( q) with Π( µ/θ), where µ is given by (1.18). The rest of this section is devoted to the proof of (3.5). For the sake of notational simplicity, we introduce the following notation for the double limit (δ, We also do not mention explicitly that we are working with subsequences and therefore, we do not relabel any sequence. Since we use only countably many relabelings, such procedure can be made rigorous by the standard diagonal procedure. Note that in all cases considered below, the order of the limit passages is not important; for the sake of clarity, we will assume that we let first δ → 0 and afterwards ε → 0. To end this first part, we also introduce the truncation function T k : R + → R + of class C 1 (R + ), which will be needed later. For arbitrary k ∈ N and z ≥ 0, we set T k (z) := kT 1 (z/k).

3.2.
Effective viscous flux. We first focus on an effective viscous flux identity. We follow the procedure developed in [9,22] very closely, and the proof is presented here for the sake of completeness.
Then it holds for every k ∈ N and T k , defined in (3.6), that Proof. Thanks to our assumptions and since we consider the proper, not relabeled subsequence, all terms in (3.7) are well defined. We introduce an auxiliary function φ δ as the solution to As ∂Ω is of class C 2 , we have φ δ ∈ W 2,q (Ω) for all q < ∞. (The proof would also work for open bounded domains Ω by arguing locally, since the regularity holds true away from the boundary.) We set ψ ψ ψ δ := ∇φ δ and ψ ψ ψ := ∇φ. The convergence properties of φ δ yield the strong convergence ψ ψ ψ δ → ψ ψ ψ in L q (Ω; R 3 ) and the weak convergence ∇ψ ψ ψ δ ⇀ ∇ψ ψ ψ in L q (Ω; R 3×3 ) for all q < ∞.
and therefore, we can apply the div-curl lemma (see Lemma A.2) to the matrix-valued functions T δ − ρ δ v δ ⊗ v δ and ∇ψ ψ ψ δ . Since the divergence of the former sequence is bounded in L r (Ω; R 3 ) for some r > 1 and the curl of ∇ψ ψ ψ δ vanishes, we infer that Since (Ω) for q < 6 and ρ δ ⇀ ρ weakly in L γ+ν (Ω) for γ > 3/2, the product converges weakly to the product of the limits, i.e. Hence, which is the starting point of further investigations. First, we focus on the term involving the tensorial product of velocities. Note that the sum of (1.26) gives div(ρ δ v δ ) = 0 in the weak sense and that for all q < 6. Therefore, using the div-curl lemma again, (and thus a.e. in Ω), where the second equality follows from the a.e. convergence of (v δ ). Hence, we deduce from (3.10) that Recall that T δ = S δ − p δ I = S δ − pI. In view of the definitions of ψ ψ ψ δ and ψ ψ ψ, this shows that Finally, by the definitions of φ δ and φ (see (3.8) and (3.9), respectively), we obtain The left-hand side corresponds to that one of (3.7). It remains to identify the terms on the right-hand side. The right-hand side is uniquely defined, so we just need to identify it almost everywhere in Ω. Using convergences (3.2) and the Egorov theorem, we can find for any ε > 0 a measurable set Ω ε ⊂ Ω such that |Ω\Ω ε | ≤ ε and θ δ → θ strongly in L ∞ (Ω ε ). Consequently, using definition (1.11) of the viscous stress tensor, the previous convergence result, and convergences (3.2) again, we can identify the weak limits in Ω ε and conclude (with the help of (3.8) and (3.9)) that Then, substituting these expressions into (3.11), But since |Ω \ Ω ε | ≤ ε for any ε > 0, relation (3.12) holds a.e. in Ω.
Thus, it remains to identify the first term, i.e., we wish to relate the difference involving ∇ 2 φ δ and ∇ 2 φ with a difference involving ∆φ δ and ∆φ. For this, let χ ∈ C ∞ 0 (Ω) be a test function. In the following, we use only formal manipulations which can, however, be justified by approximation by smooth functions. It follows that The second term vanishes because of the strong convergence of v δ . Thus, integrating by parts twice, and the second term again vanishes because of the strong convergence of v δ . As the test function χ is arbitrary, we infer that Thus, inserting this expression into (3.12) and using the properties of φ δ and φ, i.e., (3.8) and (3.9), respectively, we deduce (3.7) a.e. in Ω.

3.3.
Estimates based on the convexity of the free energy. In this part, we show that the left-hand side of the effective viscous flux identity (3.7) gives us the important description of the possible oscillations of the total density ρ δ , provided we assume the convexity of the free energy with respect to the partial densities. This is summarized in the following lemma. Lemma 3.3. Let the free energy satisfy the hypothesis (H6) and the sequence ( ρ δ , θ δ , µ δ , p δ ) with θ δ > 0 a.e. in Ω fulfil ρ δ ⇀ ρ weakly in L 1 (Ω; R N ), (3.13) p δ ⇀ p weakly in L 1 (Ω), (3.14) θ δ → θ strongly in L 1 (Ω), (3.15) ln θ δ → ln θ strongly in L 1 (Ω), in Ω, (3.19) where K 2 > 0 is given in (1.20). Proof.
Furthermore, introducing the linear mapping P : R N → R N by P( u) i := u i − u N for i = 1, . . . , N, we know that |P( u)| ≤ C|Π u| (see Lemma A.5 and Remark A.6), and due to the linearity of P, we deduce from (3.20) that Therefore, we just need to show that (3.18)-(3.19) holds true a.e. in Ω η , and since |Ω\Ω η | ≤ η, the equalities will necessarily hold also a.e. in Ω.
Step 2: Using the conjugate h * θ . We characterize W k as in Ω.
Thus, we need to compute the differences p( ρ ε , θ ε ) − p( ρ δ , θ δ ) and T k (ρ ε ) − T k (ρ δ ) for ε > 0, δ > 0. As the pressure is equal to the Legendre transform of h θ , we can write The first term Y 1 ε,δ is formulated as where z t := t µ ε + (1 − t) µ δ . Using the decomposition In a very similar way, using the fact that ρ ε = ∇h * θε ( µ ε ), we find that and rewrite Z 1 ε,δ as where we defined for arbitrary 0 ≤ s ≤ 1 and k ∈ N, The sum over i, j = 1, . . . , N in (3.25) can be rewritten, using (3.22), as Step 3: Proof of (3.18). We restrict our analysis to Ω η , since for every fixed k ∈ N, which follows from assumption (3.14). Furthermore, we define the sets Similarly as before, we again have for every fixed k ∈ N, Therefore, we focus on the behavior of the sequence on the sets Ω ε,δ η,R . We use Hypothesis (H6) to show that some terms Y j ε,δ and Z j ε,δ vanish in the limit. Since we assume that ρ ε (x) + ρ δ (x) ≤ R for x ∈ Ω ε,δ η,R , it follows on this set that Consequently, since we know that θ δ , θ ε , and θ are bounded from below and above in Ω η , Hypothesis (H6) implies that where ω refers to the modulus of continuity of h * θ , introduced in (1.19). Finally, thanks to the continuity, for any t ∈ (0, 1) and z t = t µ ε + (1 − t) µ δ , and Hypothesis (H6) again implies that in Ω ε,δ η,R .
With these auxiliary results, we can now focus on the limiting process. Let the function χ ∈ L ∞ (Ω) be arbitrary and nonnegative. It follows from definition (3.21), the uniform convergence of (θ ε ) in Ω η , estimate (3.29), and Hypothesis (H6) that Similarly, again thanks to Hypothesis (H6) and the proper definition of the set Ω ε,δ η,R , lim sup η,R ) = 0. Hence, using the weak convergence (3.14) and definitions (3.21) and (3.24), where we identify the terms I j ε,δ for j = 1, 2, 3 with the help of (3.23) and (3.26) as We start with the easiest term, which is I 3 ε,δ . We deduce from (3.30) and (3.31) and the uniform convergence of Π µ ε in Ω η (and consequently also of P( µ ε )) that and then it obviously follows that I 1 ε,δ ≥ 0. It remains to analyze the integral I 2 ε,δ . For this, we use the Cauchy-Schwarz and Young inequalities for some κ > 0, and the positive semi-definiteness of (∂ 2 jℓ h * θ ): where in the last line we used (3.30)-(3.31). As κ > 0 can be taken arbitrarily small, we use the uniform convergence of (P( µ ε )) and (3.33) to deduce that for any χ ∈ L ∞ (Ω), This inequality, together with (3.27) and (3.28), shows that W k ≥ 0. Moreover, since for all k we know that Λ k (s) ≤ Λ k+1 (s) (it follows from T ′ k ≤ T ′ k+1 ), we see from the definition of I 1 ε,δ that W k ≤ W k+1 .
Step 4: Proof of (3.19). Here, we can repeat the arguments from Step 3 almost step by step. Indeed, we have Again, we just need to identify the inequality on the set Ω ε,δ η,R , since the remaining parts vanish due to (3.13). Proceeding in the same way as in (3.24), we write Hence, repeating the procedure from the previous step, we deduce that for any bounded nonnegative function χ, The next lemma combines the results coming from the convexity of h θ and the effective viscous flux identity, which finally lead to a uniform bound on (W k ), defined in (3.17).

Using (3.7) and (3.19) again, we observe that
It follows from the Cauchy-Schwarz inequality that We infer that finishing the proof.
The following lemma is the key step in the proof. It shows that we can use the truncation function T k (ρ) instead of ρ in all estimates, and this change creates only a small error, which can be neglected in a suitable topology.

Lemma 3.5. Let all assumptions of Lemmata 3.3 and 3.4 be satisfied and let
The quantities are nonnegative and satisfy for all k ∈ N, Proof. For the identification of Q k , we can repeat the proof of Lemma 3.3 step by step, where instead of (3.24)-(3.25), we use a similar computation for ( ). Heuristically, this means that we replace Λ k , which corresponds to the derivative of T k (s), by Λ k , which corresponds to the derivative of T k (s) − sT ′ k (s). This leads to the identification of the limit (compare with (3.34)): where the term corresponds to (3.32), just replacing the function Λ k by Since T 1 is concave, we have Λ k ≥ 0. An elementary computation shows that there exists C > 0 such that for all s ≥ 0, 2s)), which implies that Λ k (s) ≤ C(Λ 2k (s) − Λ k/2 (s)). Therefore, Finally, by Lemma 3.3, (W k ) is monotone and, by Lemma 3.4, (W k /(λ 1 (θ) + λ 2 (θ))) is bounded in L 1 (Ω) and monotone as well. We deduce from the monotone convergence theorem that (W k /(λ 1 (θ) + λ 2 (θ))) is strongly converging in L 1 (Ω) which implies that Inequality (3.41) then directly leads to Now, we focus on O k , defined in (3.39). Since T k is concave, O k is nonnegative. First, we show that O k → 0 a.e. in Ω. As ρ ε ⇀ ρ weakly in L 1 (Ω; R N ), the sequence (ρ ε ) is uniformly equiintegrable. Therefore, by the weak lower semicontinuity of the L 1 (Ω) norm, Because of the weak convergence of ( ρ ε ), it holds that |{ρ > k}| + |{ρ ε > k}| ≤ C/k and consequently, the uniform integrability of (ρ ε ) shows that which implies for a subsequence that O k → 0 a.e. in Ω. Second, we show the proper bound on O k , which will enable us to use the dominated convergence theorem. We deduce from the weak lower semicontinuity of the L 2 (Ω) norm that Substituting the algebraic inequality Consequently, (3.19) implies that a.e. in Ω.
Arguing as before, the right-hand side is convergent in L 1 (Ω), and we can apply the dominated convergence theorem to conclude the second part of (3.40).

3.4.
Renormalized continuity equation. We prove that the weak limit ρ of (ρ δ ) is a renormalized solution to the mass continuity equation. Note that the proof is different to the standard ones. Indeed, for our purpose, we just need the pressure having at least linear growth with respect to the density, which is not the case in other works.
Proof. We use (3.42) with b = T k : pass to the limit δ → 0, and use the fact that v δ → v strongly in L 2 (Ω) to infer that which is the weak formulation of div(T k (ρ δ )v) + (ρ δ T ′ k (ρ δ ) − T k (ρ δ )) div v δ = 0 in the sense of distributions. Since T k (ρ δ ) ∈ L ∞ (Ω) and v ∈ W 1,2 (Ω), this equation can be renormalized; see, e.g., [28,Theorem 10.29]. It follows for any Lipschitz continuous function β with compactly supported derivative and any function ψ ∈ W 1,∞ (Ω) that (3.43) To identify and estimate the term on the right-hand side, we use definition (3.38) of Q k , the effective viscous flux identity (3.7), and a similar identity with ρ δ T ′ k (ρ δ ) instead of T k (ρ δ ): We insert this relation into (3.43) to obtain Our goal is to perform the limit k → ∞ in (3.44) to recover (3.42) with (ρ δ , v δ ) replaced by (ρ, v). We first show that T k (ρ δ ) → ρ as k → ∞ a.e. in Ω. To this end, we recall the weak lower semicontinuity of the L 1 (Ω) norm leading to Since (ρ δ ) is weakly convergent in L 1 (Ω) (assumption (3.13)), we have |{ρ δ ≥ k}| ≤ C/k, and since the sequence is uniformly equiintegrable, we deduce that which shows the pointwise convergence of T k (ρ δ ) (at least for a subsequence). Similarly, The concavity of T k implies that |T k (s) − sT ′ k (s)| ≤ T k (s) and consequently, in Ω. Thus, since b is bounded and its derivative is compactly supported, there exists C > 0 only depending on b (but not on k) such that

This estimate and (3.45)-(3.46) lead to
for any q < ∞. These convergence results, estimate (3.40) for Q k , and the boundedness of b ′ allow us to pass to the limit k → ∞ in (3.44) which gives (3.42) for (ρ, v).

3.5.
Strong convergence of the densities. We finally show that ρ δ → ρ a.e. in Ω. First, we prove that the total mass density ρ δ converges pointwise and then, exploiting the compactness of Π µ δ and the convexity of the free energy, we deduce the compactness of the sequence ( ρ δ ). We start with the convergence of the total mass densities.  Proof. The idea of the proof is to use the effective flux identity (3.7) and combine it with the properties of W k , defined in (3.17). We prove the key property We claim that if (ρ, v) or (ρ δ , v δ ), respectively, are renormalized solutions, then For this purpose, we choose ψ = 1 in (3.42) such that for any bounded b with compactly supported derivative, we have Let (f m ) be a sequence of compactly supported functions f m satisfying f m ր T k pointwise as m → ∞. Choosing f = f m in the previous identity and passing to the limit m → ∞ leads to (3.49). Of course, the same argument applies to (ρ δ , v δ ). Thus, and (3.48) holds true. Thus, using (3.7) and (3.17) in (3.48), it follows that By Hypothesis (H3), this yields lim k→∞ Ω W k dx = 0. It follows from W k ≥ 0 that lim k→∞ W k = 0, and the monotonicity W k ≤ W k+1 implies that Hence, (3.19) shows that θ(ρ δ T k (ρ δ ) − ρT k (ρ δ )) = 0 in Ω such that, since θ > 0 a.e. in Ω, Finally, using the weak lower semicontinuity of the L 1 (Ω) norm, the Cauchy-Schwarz inequality, and (3.50), The left-hand side is independent of k and therefore, we may perform the limit k → ∞ also on the right-hand side. Then, thanks to uniform equiintegrability of (ρ δ ), we deduce (3.47), which finishes the proof.
We end this subsection, and also the proof of Theorem 3.1, by the following lemma, which shows the compactness of the vector ρ of mass densities. In addition, we show that either the total mass density vanishes, i.e., all partial mass densities are zero, or all partial mass densities are strictly positive a.e. in Ω.  Moreover, for a.e. x ∈ Ω and all i = 1, . . . , N, Proof. Taking into account the strong convergence results, derived in this section, Egorov's theorem implies that for any η > 0, there exists a set Ω η such that |Ω \ Ω η | ≤ η and We collect some facts about the free energy density h θ . Since it is strictly convex with respect to ρ and ∇ ρ h θ is a strictly monotone invertible mapping from R N + to R, for every κ > 0, there exists C(κ) > 0 such that for all θ ∈ (κ, κ −1 ), (3.53) Finally, we focus on the strong convergence of the densities. We consider two cases: the total density vanishes or the total density remains positive. The first case is easy, so we assume that the total density is positive. We define the set Ω η,κ := {x ∈ Ω η ; ρ(x) ≥ 2Nκ}.
To finish the proof of Theorem 3.1, recall that Since we cannot exclude that ρ = 0 on a set of positive measure, we redefine µ and set µ := 0 on {ρ = 0}. On the other hand, we know that in Ω, and therefore also in Ω.

The approximate scheme
4.1. Auxiliary results. The following property is needed in the construction of the approximate solution.
The following lemma plays an important role in the construction of the approximate solutions.

4.2.
Formulation of the approximate equations. In this subsection, we specify the approximate system, leading to a sequence of smooth, approximate solutions.

4.2.1.
Internal energy balance. We replace the total energy balance (1.28) by the internal energy balance (2.3), whose weak formulation reads as

4.2.2.
Parameters of the approximate problem. For the sake of clarity, we state here the parameters employed in the approximate scheme. A more detailed explanation is given in the part that follows: ε ∈ (0, 1): lower-order regularization in all equations; ξ ∈ (0, 1): higher-order regularization in the internal energy equation; δ ∈ (0, ξ/2): higher-order regularization in the partial mass densities equations; χ ∈ (0, 1): quasilinear regularization in the partial mass densities equations; n ∈ N: Galerkin approximation in the momentum equation; η ∈ (0, 1): regularization in the free energy and heat conductivity.
We point out that h (η) θ satisfies the property from Lemma 4.1 and Hypothesis (H6); in particular, (1.22) holds with the same constants, since the mixed second-order derivatives of the free energy are not changed by the terms added in (4.4).
We also need to regularize the heat conductivity:

4.3.
Existence of a solution to the approximate equations. We formulate (4.6)-(4.8) as a fixed-point problem for a suitable operator.
To this end, we apply the Leray-Schauder theorem.
Indeed, the compactness of F follows from the fact that the image of F is bounded in V 0 (consequence of Lax-Milgram's Lemma), the compact Sobolev embedding H 2 (Ω) ֒→ W 1,4 (Ω) (which holds, e.g., for d = 3), and the fact that X n is finite dimensional.
It remains to prove that the set 1]. If this is shown, the existence of a fixed point for F in V follows from the Leray-Schauder theorem. 1]. This means that the following relations hold: for all φ 1 , . . . , φ N ∈ H 2 (Ω), u ∈ X n , and φ 0 ∈ H 2 (Ω). We need to show that ( q, v, w) can be uniformly bounded in V with respect to σ. To find such estimate, we derive global entropy and energy inequalities.

4.3.2.
Global entropy inequality for approximate solutions. Let us choose φ i = q i − q 0 (i = 1, . . . , N) and φ 0 = − exp(−w) = −1/θ in (4.13) and (4.15), respectively, and sum the equations. By arguing like in the derivation of (2.7), we find that where we defined We deduce from the construction of q 0 and R i (Lemma 4.2) and the convexity of h Furthermore, straightforward computations and the assumption δ ≤ ξ/2 show that for some constant c > 0. Putting (4.16)-(4.17) together gives the global entropy inequality: However, Ω N i=1 R i dx = |Ω|ρ by construction (see Lemma 4.2), thus (4.18) shows that

4.3.4.
Global total energy inequality for approximate solutions. We choose now φ i = − 1 2 |v| 2 (i = 1, . . . , N), u = v, φ 0 = 1 in (4.13)-(4.15) and sum the equations. By arguing like in the derivation of (2.8), we obtain It follows from (4.18) that while the Cauchy-Schwarz inequality, definition (4.3), and (4.18) yield Therefore, (4.20) leads to the global total energy inequality (for approximate solutions): The right-hand side of (4.21) is estimated as  We prove an estimate for ∇ ρ ∈ L 2 (Ω; R N ×3 ). On the one hand, Young's inequality gives On the other hand, we have We conclude from the these inequalities, Hypothesis (H6), more precisely (1.22), and (4.4) that As q and log θ are essentially bounded functions, also ρ is essentially bounded (this bound is uniform with respect to δ and n). Therefore, we deduce from the previous estimate and bound (4.28) that In view of Lemma 4.2 and the ε-uniform L ∞ (Ω) bound for log θ, q 0 = q 0 (log θ) is bounded in ε, so we infer from (4.25) the desired gradient bound for ρ: where the constant C > 0 is independent of ε.

5.
Weak sequential compactness for approximate solutions 5.1. Limits δ → 0, n → ∞. The bounds (2.11)-(2.13) and (4.23)-(4.32) allow us to extract subsequences (not relabeled) such that, as δ → 0, q (δ) ⇀ q weakly in W 1,4 (Ω; R N ), q (δ) → q strongly in L ∞ (Ω; R N ), log θ (δ) ⇀ log θ weakly in W 1,4 (Ω), log θ (δ) → log θ strongly in L ∞ (Ω), The L ∞ bound of log θ (δ) implies that θ (δ) is bounded away from zero and infinity. The W 1,4 (Ω) bound yields the strong convergence of (a subsequence of) θ (δ) and thus the convergence log θ (δ) → log θ a.e. and uniformly. As a consequence of the strong L ∞ (Ω) convergence of q (δ) and log θ (δ) , we have The previous convergences are sufficient to pass to the limit δ → 0 in (4.6)-(4.8) and to obtain an analogous system of equations without the terms proportional to δ. The terms that are nonlinear in ∇ q, ∇ log θ are monotone, so we can use Minty's monotonicity technique to identify the limit of these terms. The limit n → ∞ is carried out in the same way as the limit δ → 0. We have chosen not to take the two limits simultaneously because of the term δK(n) Ω |v| 2 v · u dx in (4.14), which vanishes as δ → 0 (n ∈ N fixed) but not necessarily when both δ → 0, n → ∞. The two limits could be carried out simultaneously provided that one assumes that δK(n) → 0.

5.2.
Limit χ → 0. Since the estimate of ∇ q (χ) in L 4 (Ω; R N ×3 ) blows up when χ → 0, we lose the control of ρ in L ∞ (Ω; R N ). Therefore, we need to establish χ-independent estimates of ρ, which is the goal of the following lemma.

5.2.1.
Estimate for ρ via the Bogovskii operator. The following result is the analogue of Lemma 2.3 for the approximate equations.
As before, it is not difficult to justify that θ (χ) → θ a.e. in Ω. We deduce from the continuity of the free energy h θ that ρ (χ) is also a.e. convergent in Ω. Bound (5.4) then implies that ρ (χ) → ρ strongly in L r (Ω; R N ) for all r < 12 11 L.

5.2.3.
Entropy and total energy balance equations. We will derive the balance equations for the entropy and total energy. It was not possible to derive the latter equation as long as the p-Laplacian regularization in terms of q was included in the mass densities equation. Let ψ ∈ W 1,∞ (Ω) with ψ ≥ 0 be arbitrary. By choosing φ i = (q i − q 0 )ψ in (5.6) and φ 0 = − 1 θ ψ in (5.8) and proceeding as in the derivation of (4.18), we obtain Next, let ϕ ∈ W 1,∞ (Ω). By choosing φ i = − 1 2 |v| 2 ϕ in (5.6), u = vϕ in (5.7), and φ 0 = ϕ in (5.8) and proceeding as in the derivation of (4.21), it follows that

5.3.2.
Estimates based on the convexity of the free energy. Lemmata 3.3-3.5 also hold for the approximate solutions, as they only require the properties of the free energy stated in the introduction, which still hold true for the regularized free energy introduced in (4.4), and the estimates in Lemmata 2.2-2.4, which are satisfied by the approximate solutions.
We replace ψ by b ′ (T k (ρ (ε) ))ψ in (5.12) and exploit the definition of the truncation operator T k to find that It is proved in Lemma 3.6 that the left-hand side of (5.13) converges to the left-hand side of the renormalized continuity equation (3.42) as k → ∞. Therefore, it is sufficient to prove that the right-hand side of (5.13) tends to zero as k → ∞. Indeed, since b ′ is bounded and (4.25) and (4.34) hold, it follows that meaning that the right-hand side of (5.13) tends to zero as k → ∞. We conclude that the weak limit (ρ, v) of (ρ (ε) , v (ε) ) satisfies the renormalized continuity equation (3.42).
This finishes the part related to the limit ξ → 0.

Appendix A. Auxiliary results
For the convenience of the reader, we collect some results needed in this paper. For the first lemma, which is proved in [28,Theorem 10.11], we introduce the space L p 0 (Ω) = u ∈ L p (Ω) : Ω u dx = 0 , 1 < p < ∞.
Theorem A.1 (Bogovskii operator). Let Ω ⊂ R d (d ≥ 2) be a bounded domain with Lipschitz boundary and let 1 < p < ∞. Then there exists a bounded linear operator B : L p 0 (Ω) → W 1,p 0 (Ω) such that div B(u) = u in the sense of distributions for all u ∈ L p 0 (Ω), and there exists C > 0 such that for all u ∈ L p 0 (Ω), B(u) W 1,p (Ω) ≤ C u L p (Ω) .
We need the following version of the Korn inequality, proved e.g. in [18, Lemma 2.1].
Lemma A.5. Let P ∈ R N ×N be a matrix such that where p i is the i-th column of the matrix P. Introduce the matrix A ∈ R N ×N with A ij = δ ij − δ N j for i, j = 1, . . . , N. Then there exists a constant C > 0 such that for all v ∈ R N , where | · | denotes the Euclidean norm in R N .
The matrix M ε := P T P − εA T A is clearly symmetric and M ε w = 0 for all w ∈ V . However, the matrix P is nonsingular on V ⊥ , thus inf v∈V ⊥ , | v|=1 |P v| 2 > 0. Now, it is enough to choose ε < inf v∈V ⊥ , | v|=1 |P v| 2 sup | v|=1 |A v| 2 . Then M ε is positive definite, and the claim follows.
Remark A.6. Note that if the vector v depends on a spatial variable x, then |P∇ v| ≥ C|A∇ v| for all v ∈ R N .

Appendix B. Example for the free energy
We show that the free energy density where γ > 1, satisfies Hypothesis (H6). We compute the partial derivatives where ∂ i = ∂/∂ρ i , ∂ 2 ij = ∂ 2 /(∂ρ i ∂ρ j ), and ∂ θ = ∂/∂θ. An induction argument with respect to N shows that the Hessian is positive definite for ρ ∈ R N + and θ ∈ R + . Moreover, h θ ( ρ) is smooth with respect to θ in this region and satisfies growth conditions (1.22), (1.23), (1.24) as well as (1.21). We write the Hessian D 2 h θ ( ρ) as the sum of the two matrices A and B with coefficients By construction of h * θ , we have D 2 h * θ ( µ) = D 2 h θ ( ρ) −1 | ρ=Dh * θ ( µ) . Therefore, (B. We set M := B −1/2 AB −1/2 . We are able to compute M explicitly by using the power series of the function (1 + z) −1 . Since where n = N i=1 n i = N i=1 ρ i /m i , the coefficients (M ℓ ) ij of the matrix M ℓ with ℓ ∈ N read as Therefore, we obtain