{'fraction_non_alphanumeric': 0.06699043860980422, 'fraction_numerical': 0.02249203217483685, 'mean_word_length': 3.624649073554183, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Incompressible limit of the Navier-Stokes model with a growth term

Starting from isentropic compressible Navier-Stokes equations with growth term in the continuity equation, we rigorously justify that performing an incompressible limit one arrives to the two-phase free boundary fluid system.


Introduction
The purpose of this work is to analyze the Navier-Stokes equations that generalize the fluid-based models of tumors. In the mathematical literature, tumor growth has been modelled using various microscopic and macroscopic models [5]. At the macroscopic level, we may distinguish between models which describe the tumor growth through the dynamics of its cell density, and free boundary models in which tumor is described by its geometric domain subjected to mechanical constrains [12]. From the mechanical viewpoint living tissues may be considered as fluids [7]. In the simplest approach the dynamics of cell density is governed by cell division and mechanical pressure. Depending on the modelling assumption, and the complexity of the model, mechanical pressure is incorporated in the fluid velocity through Darcy's law, Stokes' law, Brinkman's law or Navier-Stokes' law (see e.g. [29,28,17,16,6]). Notice that Darcy's law, Stokes's law or Brinkman's law may be derived at least formally from Navier-Stokes' law, and so, the latter may be considered as a generalization of the other models.
In this paper we perform mathematical analysis of the Navier-Stokes model with the growth term as for the models of tumor. We are particularly interested in the stiff pressure law limit, often referred to as incompressible limit. The limiting model is a free boundary compressible/incompressible system of fluid equations. Derivation of the free boundary models from cell mechanical models has been the subject of many recent contributions in the field of tumor growth modelling [24,25,26,14,15,21]. These models identify tumor with an area of the incompressible (constrained) fluid, while the surrounding healthy tissue can be viewed as a compressible (unconstrained) fluid. In almost all aforementioned works, for simplicity, Darcy's law is used as a closure relation for the system. This means that the velocity is proportional to the gradient of the mechanical pressure, which results in a porous-like type of system. In [26], Birkman's law was used to model the tumor as a visco-elastic medium, see also [2], and [27] for the model of growth of tissue in which cell division and apoptosis introduce stress sources that, in general, are anisotropic. The aim of this work is to extend these works to more general relation between the velocity and the pressure, namely the Navier-Stokes equation. The unknowns are ρ(t, x), the cell density, and u(t, x), the macroscopic velocity field, depending on the time t > 0 and the position x ∈ R d . Our starting system reads as follows where p is the pressure and µ > 0, µ + ξ > 0 are the viscosity coefficients. For future use, we provide also the nonconservative form of the equation (1b), which reads ρ (∂ t u + u · ∇u) − µ∆u − ξ∇ div u + ∇p = 0.
The right hand side in (1) represents the growth term depending on the pressure, we assume that the quantity P M is often refered to as the homeostatic pressure [27]. As in [24,25,26], we choose the barotropic pressure law with the exponent γ that might be very big.
The system (1) is complemented with initial data ρ(0, x) = ρ 0 (x), (ρu)(0, x) = m 0 (x), which are chosen such that for any large enough γ, uniformly with respect to γ. Moreover, we prescribe the values of u and ρ at infinity: with the relevant compatibility condition for the initial data. When γ is fixed, the system (1) is the compressible Navier-Stokes system with additional terms on the right hand side of the continuity equation (1a) and in the momentum equation (1b). The purpose of this paper is to rigorously justify the so-called stiff pressure law limit, i.e. γ → +∞ which leads to the two phase compressible/incompressible system This system is complemented with the same initial data ρ 0 , m 0 as system (1). The limit of this type was first considered for the compressible Navier-Stokes equations without any additional growth terms by Lions and Masmoudi [20]. Similar limit passage was also recently investigated for polymeric fluids [10]. We would also like to remark that the two-phase models of the type (7) can be obtained as the limit of the compressible Navier-Stokes system with the singular pressure, see [23,9]. Compared to the case without growth term, the main difficulty lies in obtaining strong convergence for the density. Indeed classical approach developed by Lions [19] and Feireisl [11] fails precisely due to the presence of the growth term. Therefore, we follow a recent strategy proposed by Bresch & Jabin [3] for the compressible Navier-Stokes equation (see also [4]) and adapt it to the case at hand.
Before stating our main result, let us explain formally how the system (7) may be obtained from (1). We assume existence of a sequence denoted by n, such that for n → ∞, γ n → ∞, and ρ n → ρ ∞ , u n → u ∞ , ρ γn n → p ∞ strongly. Writing (4) as p n = ρ n p γn−1 γn n and letting n → ∞ we check that ρ ∞ , p ∞ satisfy the relation (7d).
Let us now introduce the set Ω = {p ∞ > 0} ⊂ R d , we have two cases: which is the compressible pressureless Navier-Stokes system with the source term.
• On Ω, we deduce from (7d) that we have ρ ∞ = 1. Then (7a)-(7c) reduces to div u ∞ = G(p ∞ ), which might be seen as the incompressible Navier-Stokes system. Note that from the expression of G in (3), we may rewrite the last system as Therefore the limit system (7) reveals the features of both: compressible and incompressible fluid equations with the free interphase separating Ω from R d \ Ω.
We conclude the introduction by explaining the link between the system (7) and the Hele-Shaw system for tumor growth. Neglecting the acceleration term and assuming that the viscous resisting force is proportional to the velocity, then the momentum equation in system (9) reduces to This is the so-called Darcy's law. Inserting this equation into the first equation in (9), we recover the Hele-Shaw system for tumor growth, −∆p ∞ = ν 0 G(p ∞ ) on Ω.

The main result
Our main result concerns the convergence of weak solutions of system (1) to weak solutions of the system (7). Before formulating the main theorem let us first specify the notion of solutions.
Definition 2.1 (Weak solution of the primitive system) Suppose that the initial conditions are as in (5). Let T > 0. We say that the couple (ρ, u) is a weak solution of the problem (1), and we have: for all test functions ψ ∈ C 1 c ([0, T ] × R d ); (iii) there exists a constant C > 0 such that the energy inequality holds for a.a t ∈ (0, T ), where Definition 2.2 (Weak solution to the limiting system (7)) Let and we have: The compactness of the sequence of weak solutions to system (1) with γ n → ∞ is guaranteed by our main theorem.
n=1 be a sequence of weak solutions to system (1) with p(ρ n ) = ρ γn n , in the sense of Definition 2.1. Then, up to extraction of a subsequence, the limit of {(ρ n , u n , ρ γn n )} ∞ n=1 for n → ∞ solves (7) in the sense of Definition 2.2. More precisely, there exist ρ ∞ , u ∞ , p ∞ such that: The existence of solutions to the primitive system (1) is a combination of nowadays classical techniques and compactness argument used in the proof of Theorem 2.3, therefore it is postponed and only roughly discussed in the end of the paper in Section 6. Otherwise, the paper is organized as follows. In Section 3 we derive the a-priori estimates, i.e. the estimates that can be obtained for the weak solutions of system (1) and are uniform with respect to parameter γ. Then, in Section 4 we present the main compactness argument implying the pointwise convergence of the sequence ρ n . In Section 5, we show that the a-priori estimates and the compactness argument are sufficient to pass to the limit in (1) to obtain (7) which finishes the proof of Theorem 2.3.

A-priori estimates
The estimates presented in this section are derived using the assumption that (ρ n , u n ) is sufficiently smooth solution of (1). This is not necessarily true for the weak solutions from Definition (2.1). However, the calculations can be made rigorous on certain level of approximation discussed in Section 6.

The energy estimate
Let us denote the energy and the energy dissipation (13), (14) corresponding to ρ n , u n , and p(ρ n ) = ρ γn n by E n , J n , respectively. The following a-priori estimates are then uniform with respect to n.
Lemma 3.1 Under assumptions (5) and (3), let T > 0 be fixed, then we have the following a priori estimates, uniform in n ∈ N: with E n (t) and J n (t) defined in (13), (14).
Proof. The proof of (i) is standard. By Stampacchia method we show that the nonnegativity principle holds; since ρ 0 (x) ≥ 0, we deduce that ρ n (t, x) ≥ 0 for any time t > 0. Thus p n (t, x) = ρ γn n (t, x) ≥ 0. Then by a simple integration of (1a) where we use (3). We then conclude by integration in time and the Gronwall inequality. For part (ii), we compute where we used (1a) for the first two terms and (2) for the last term. Noticing that ∇p n = γ n ρ γn−1 n ∇ρ n , and u n · ∇|u n | 2 = 2(u n · ∇u n ) · u n , we may cancel the second integral with the last two terms of the last integral. Then, integrating by parts, we deduce Since p n ≥ 0, and G satisfies (3), we have G(p n ) ≤ G 0 P M , then Moreover, still using assumption (3), we have that G(p n ) ≤ 0 if p n ≥ P M ⇐⇒ ρ n ≥ P 1/γn Thus, using the bound on the L 1 norm of ρ n from part (i), we deduce that there exists a nonnegative constant C such that uniformly in n we have and we conclude using the Gronwall lemma.
(iii) As a consequence of the point (ii) above, we deduce Then, sup By interpolation, for any q ∈ (1, γ n ), we have with 1 q = θ n + 1−θn γn , then when n → +∞, we have θ n → 1 q . We deduce that for N large enough, the sequence {ρ n } n≥N is uniformly bounded in L ∞ (0, T ; L q (R d )).
Lemma 3.1 implies that we may extract a subsequence (still labelled by n), such that ρ n ⇀ ρ ∞ weakly as n → +∞. Then, by lower semi-continuity of norm, passing into the limit n → +∞ in (19), we have Since this latter estimate is true for any q ≥ 1, we may let q going to +∞ to find In addition, we can prove Lemma 3.2 Under the same assumptions as in Lemma 3.1, we have (ρ n − 1) + → 0 as n → +∞.
Proof. Let us introduce φ n = (ρ n − 1) + . From the energy estimate (17), we deduce It has been proved in [20, p.24] that, for any q > 1 and any x ≥ 0, there exists a q > 0 such that (1 + x) k ≥ 1 + a q k q x q , for k large enough. Thus, for sufficiently large n, we have Therefore φ n → 0 strongly in L q ((0, T ) × R d ) for any q ≥ 1 .

The estimate of the pressure
Note that the energy estimate from the previous section does not provide any estimate of the pressure independent of n, only the estimate of the density. In the following lemma, we state an L 2 estimate on the pressure.
Proof. From the renormalization property [8] for equation (1a), we have that for any C 1 function Let K > 0 be a nonnegative constant. We define, for γ n > 1, β K the function For all y ≥ 0 and γ n > 1, we have Using (20) with β = β K and inserting the assumption (3) on the growth function G, we deduce where we used (21) to get the last inequality. On the set {ρ n ≤ 1}, we clearly have ρ n β ′ K (ρ n ) ≤ γ n ρ n . Then, Integrating and using the Cauchy-Schwarz and the Young inequalities, we deduce that for any ǫ > 0, there exists a nonnegative constant C ǫ such that We may fix ǫ > 0 such that, from (21), Integrating in time, we obtain that there exists a nonnegative constant C such that Using estimates in Lemma 3.1, we deduce that there exists an uniform (with respect to n and K) constant C > 0 such that Therefore, for all K ≥ 0, we deduce that We may now let K go to +∞ and, by the monotone convergence theorem, we conclude the proof.
Then we may apply Lemma 6.9 from [22] and deduce that (20) holds with β(y) = y γ . We therefore obtain the evolution equation for the pressure p n = ρ γn n , ∂ t p n + u n · ∇p n + γ n p n div u n = γ n p n G(p n ).
Remark 3.5 The fact that we can derive the uniform estimates for the pressure is one of the main advantages of the growth term in the continuity equation (1a). Not having it would require more laborious estimates with the application of Bogovski type of operator, see for example [20], [23].

The estimate of the nonlinear terms in the momentum equation
Before letting n → ∞, we need to provide the uniform estimates of the rest of nonlinear terms from the momentum equation (1b). Applying the operator (−∆) −1 div to both sides of (2), we deduce where we use the notation for the total derivative D(ρ n u n ) = −(−∆) −1 div(ρ n ∂ t u n + ρ n u n · ∇u n )) .
Using the L 2 estimate of the pressure and the L 2 estimate of div u n following from the energy estimate we deduce the following fact.

Compactness
The purpose of this section is to establish the compactness of the density sequence {ρ n } ∞ n=1 . To do it, we follow the strategy proposed by Bresch & Jabin [3] (see also [4]) in the context of compressible Navier-Stokes equations with the non-monotone pressure law. We adapt their approach to whole space R d case, with a nonzero growth term in the right hand side of the continuity equation, and consequently, the conservative form of the momentum equation. Application of nowadays classical approach developed by Lions [19] and Feireisl [11] fails precisely due to the presence of this additional term.
The main result of this section is the following (1), (4) with assumptions (3), (5), such that the estimates from Lemma 3.1 and in Lemma 3.3 hold.
Then the sequence . The rest of this section is dedicated to the proof of this fact.

A compactness criterion
In order to prove local compactness for the density sequence {ρ n } ∞ n=1 we use a compactness criterion, for the proof of which we refer the reader to [1, Lemma 3.1], or [3,Proposition 4.1]. This criterion was applied to the study of Navier-Stokes equations with non-monotone pressure and anisotropic stress tensor in the aforementioned papers [3,4].
Let us first introduce the necessary notations.
We define a family {K h } h>0 of nonnegative function by ) and is compactly supported in B(0, 2). Moreover K h is equal to some function K(x) independent on h outside B(0, 3/2). We will also make use of the inequality which holds for some nonnegative constant C independent of h, thanks to our choice for K h . We also denote Then the compactness criterion states what follows.
, then the above lim sup converges to 0 as h goes to 0.

Definition of the weights
Let us define the weights w n as solutions of the transport equation complemented with the initial data w(t = 0) = 1. Here λ is some nonnegative constant which will be fixed later on. To simplify the notations, we drop the index n denoting the weight simply by w.
By M we denote the maximal operator, defined by Recall that we have the following inequality (see e.g. [30]) for any Φ in W 1,1 (R d ). Note that, thanks to Lemma 3.1 and Lemma 3.3, we have that B n defined in (25) is uniformly bounded in L 2 ([0, T ] × R d ). This allows us to deduce the following properties of the weight w.
Proposition 4.3 Let us assume that u n is given and that it is bounded in uniformly with respect to n. Then, there exists a unique solution to (25). Moreover, we have (ii) If we assume moreover that the pair (ρ n , u n ) is a solution to (1a) and ρ n is uniformly bounded , by standard theory of renormalized solutions to the transport equations [8], we may construct a nonnegative solution to (25). Moreover, since B n is nonnegative, we have clearly that w ≤ 1, since it is true initially.
We integrate it in space and use (3) to deduce Using the Gromwall lemma, we obtain Finally, since B n and ρ n are uniformly bounded in L 2 ((0, T ) × R d ), we conclude using the Cauchy-Schwarz inequality.
This computation can be made rigorous using renormalization technique [8]. We observe that thanks to our pressure law in (4), we have that sign (ρ n (x) − ρ n (y)) = sign (p n (x) − p n (y)). Then, we can rearrange the last term of the right hand side as (ρ n (x)G(p n (x)) − ρ n (y)G(p n (y)) (ρ n (x) − ρ n (y)) where we use the definition of G (3). Moreover, since p n is nonnegative, G(p n ) ≤ G 0 P M . We arrive at We then introduce where the weights w satisfy (25). Using (27) and the symmetry of K h , we deduce where − ρ n (y)) 2 (∂ t w(y) + u n (y) · ∇w(y)) dx dy, K h (x − y)(div u n (x) − div u n (y)) (ρ n (x) − ρ n (y)) ρ n (x)w(x) dx dy.

Estimate of A 1
The term A 1 is the same as in [3,4]. For the sake of completeness we recall how to estimate it below. First, we make use of the following inequality dz. Recall that D n u n ≤ M |∇u n |. For the proof we refer the reader to [13,Lemma 3.1]. Then, using inequality (24) and the symmetry of K h we get Next, we integrate in h on (h 0 , 1). Using that and changing the variables z = x − y, we may apply the Cauchy-Schwarz inequality and the uniform L 4 bound on ρ n to deduce We may bound D |x−y| u n by the Maximal operator M |∇u n |, thus The second term on the right hand side of (29) will be controlled by the term A 2 .

Estimate of A 2
From (25), we have Therefore, combining the latter equality with (29), we deduce From the definition of B n in (25), we can find λ large enough such that

Estimate of A 3
To estimate the A 3 term, we first recall the link between div u n and p n (23), and the notation D(ρu) = −(−∆) −1 div(ρ∂ t u + ρu · ∇u)) . Then, Note that since p n = ρ γn n is increasing with respect to ρ n , we have (p n (x)−p n (y)) (ρ n (x) − ρ n (y)) ≥ 0. Therefore, the first term in (31) has a good sign when moved to the left hand side.
Thus, departing from (28) and integrating in h, we use (30) and (31) to deduce The estimate of the second term in (32) follows from the following Lemma:

Remark 4.6
The only change in the statement of the above lemma with respect to Lemma 8.3 in [3] is that in our case the continuity equation has an extra production term. Note however, that the operator D(ρu) is the Riesz operator applied to the nonconservative form of the momentum transport, see (23). However, the momentum equation in the nonconservative form (2) does not include any extra contribution from G(p). This makes the proof of Lemma 4.5 the same as the proof of Lemma 8.3 from [3].
In order to apply Lemma 4.5, we need to truncate the integrant in the last integral of (32). We introduce a smooth truncation function φ : [0, ∞) → [0, 1] such that 0 ≤ φ ≤ 1, φ(x) = 1 for x ≤ 1 2 , and φ(x) = 0 for x > 1. We then split the last term in (32) into two parts Note that for some α > 0, we have since the left hand side vanishes when ρ n (t, x) ≤ L/2 and ρ n (t, y) ≤ L/2. Therefore, for the same α > 0 upon using the Cauchy-Schwarz inequality, the uniform bounds on D(ρ n u n ) in L 2 ([0, T ] × R d ) (see Corollary 3.6) and on ρ n in L ∞ (0, T ; L q (R d )) for q ∈ (1, γ n ), we obtain Then, we may apply Lemma 4.5 with the function By definition of the truncation φ, we have that Φ L ∞ ≤ CL 2 . For the control on the time derivative of Φ, we notice that Φ is a combination of functions ρ n and w which satisfy a transport equation with the same velocity field, but different right hand sides. Then, ∂ t Φ + div x (u n (x)Φ) + div y (u n (y)Φ) = f 1 div x u n (x) + f 2 div x u n (y) + f 3 B n (x) + f 4 B n (y) + f 5 ρ n (x)G(p n (x)) + f 6 ρ n (y)G(p n (y)), where B n is defined in (25) and G(p n ) is defined in (3). Every function f i contain as a factor φ(ρ n /L) or a derivative of φ. Then f i L ∞ ≤ CL 2 for i = 1, . . . , 4. We deduce that the constant C Φ in Lemma 4.5 is bounded by CL 2 . Thus, Optimizing in L, i.e. choosing L = h −θ/(α+2) , we deduce that there exists θ 0 > 0 such that Finally, integrating in time (32) and inserting (33) and (36), we obtain for all t ∈ [0, T ]

Removing the weights and compactness argument
Let η < 1. We define ω η = {x : w ≤ η} and denote by ω c η its complementary. We have with where we used the symmetry of K h and the fact that K h L 1 = 1. To treat the last integral we recall an interpolation inequality for ρ n ∈ L ∞ (0, T ; L q (R d )), where τ = q−2 2(q−1) . Therefore since for η < 1, | log w(x)| ≥ | log η| for all x ∈ ω η , and the last inequality follows by (26). Inserting these estimates on I 1 and I 2 into (38), we arrive at Finally, from (37), we deduce Finally, we obtain the compactness of the sequence {ρ n } n , as stated in Proposition 4.1, by applying the compactness criterion in Lemma 4.2. Indeed the estimate on the time derivative is a direct consequence of the conservation equation (1a) and of the energy estimate in Lemma 3.1.

Limiting system
This section is dedicated to the limit passage n → ∞ in the definition of the weak solutions to the approximate system (Definition 2.1). We will first gather together all the uniform estimates for the sequence of solutions {ρ n , u n , ρ γn n } ∞ n=1 and pass to the limit in the continuity and the momentum equation. Then we prove the complementary relation (7d). Finally, we also prove the complementary relation div u = G(p ∞ ) on the set {ρ ∞ = 1}.

Convergence in the continuity and the momentum equations
Following the estimates of Section 3 and the compactness result in Section 4, there exists (ρ ∞ , p ∞ , u ∞ ) such that, for n → +∞, up to a subsequence, we have ρ n is uniformly bounded in L ∞ (0, T ; L q (R d )), q ≥ 1, (Lemma 3.1) (42) From (41) and (42) by interpolation of the Lebesgue spaces, we deduce that In addition, the time derivative of ∂ t ρ n can be expressed by means of equation (1a), therefore the Arzelá-Ascoli theorem and the uniform estimate (42) imply that Moreover, uniformly with respect to n we have √ ρ n u n L ∞ (0,T ; and so, using also (42) we get Combining (46) with (44) we check that and therefore from the uniqueness of the weak limit ρu = ρ ∞ u ∞ , and also Using the estimates of p n , ρ n and u n , we deduce that ∂ t (ρ n u n ), given by (1b), is uniformly bounded in L 2 (0, T ; W −1,q 3 (R d )) + L 2 (0, T ; W −1,2 (R d )) + L ∞ (0, T ; L q (R d )) + L 2 (0, T ; L p (R d )), for 1 ≤ p < 2. This estimate might be used to identify the limit in (52). To this purpose, we recall the following compensated-compactness lemma, see [20,Lemma 3.3].
Then f n g n converges to f g weakly in D ′ ([0, T ] × R d ).
Taking g n = ρ n u n and f n = u n in this lemma, we justify that (52) is in fact The last task is to pass to the limit in the production term of the momentum equation ρ n G(p n )u n . To this purpose we first note that this sequence is weakly convergent in L p ((0, T ) × R d ) for some p > 1 to a limit denoted by ρ ∞ G(p ∞ )u ∞ . To identify this limit, we will use (43) and the strong convergence of the sequence {ρ n u n } ∞ n=1 . To deduce the latter, we first note that (41) and (44) imply that √ ρ n u n ⇀ √ ρ ∞ u ∞ weakly in L 2 (0, T ; L q loc (R d )) for q < 2d d−2 . Next, as in (54) we show that for any compact set K ⊂ R d we have The weak convergence of √ ρ n u n and the convergence of the L 2 -norm implies that √ ρ n u n converges to √ ρ ∞ u ∞ strongly in L 2 (0, T ; L 2 loc (R d )). From (46) and from the uniform bounds on ρ n u n in (49) it then follows that ρ n u n converges to ρ ∞ u ∞ strongly in L 2 (0, T ; L 2 loc (R d )). This concludes the proof of the passage to the limit in the continuity and in the momentum equations leading to the weak solution from Definition 2.2.

Passage to the limit in the congestion relation
Here we follow a similar argument from [20]. In order to recover relation (7d) we first see that for any δ > 0, there exists n 0 sufficiently large such that for n ≥ n 0 we have ρ γn+1 n ≥ ρ γn n − δ.
Thus, passing with n to the limit we obtain ρ γn+1 n ≥ p ∞ − δ.
The limit on the left hand side can be immediately identified with ρ ∞ p ∞ , due to the strong convergence of ρ n and weak convergence of p n . Therefore, letting δ → 0, we get Note however, that due to (45), ρ ∞ ≤ 1, therefore ρ ∞ p ∞ ≤ p ∞ , which implies that ρ ∞ p ∞ = p ∞ .

Consistency relation
In the following lemma, we show that conditions (8) and (9) are compatible. This is provided by the equivalency of the following conditions.
where ρ ≥ 0 a.e. in (0, T ) × R d satisfy the transport equation Then the following two assertions are equivalent (i) div u = G(p) a.e. on {ρ ≥ 1} and 0 ≤ ρ 0 ≤ 1, Proof. We follow the idea from [20, Lemma 2.1]. We first prove the implication (i) ⇒ (ii). From the renormalization property, we have that for any C 1 function β from R to R such that |β(t)| ≤ C(1+t), We choose for β the function Then we get (after regularization and passing to the limit for the rigorous justification): Denoting σ = β(ρ) − ρ and subtracting from the latter equation (55), we obtain where we used the assumption div u = G(p) on {ρ ≥ 1}. Moreover, thanks to our choice of function β, we have σ = β(ρ) − ρ = (1 − ρ)1 {ρ≥1} . Therefore, we arrive at It is classical to deduce that |σ| satisfies the same equation. Integrating it over R d , we obtain Note that σ(0) = 0, since by (i) 0 ≤ ρ 0 ≤ 1. Therefore, using the Gronwall lemma we conclude that 0 = |σ| = (1 − ρ)1 ρ≥1 which implies (ii). For the reverse implication, (ii) ⇒ (i) we proceed as follows. Since ρ is bounded, equation (56) holds for any C 1 function β. In particular, for β(ρ) = ρ k , for any integer k, we get . This means that kρ k (G(p) − div u) is a distribution bounded uniformly with respect to k. We deduce that we can pass into the limit k → ∞ we therefore obtain in the sense of distributions.

About existence of solutions
In this section we explain the main steps leading to the construction of the weak solutions from Definition 2.1. We will explain how this solution can be obtained by chain of approximations of system (1), inluding parabolic regularization of the continuity equation and the Faedo-Galerkin approximation of the momentum equation.

Existence of solutions to system with additional dissipation
The weak solution from Definition 2.1 will be obtained as a limit (ρ, u) as ε → 0 + of the weak solutions (ρ ε , u ε ) to the following system with artificial viscosity The existence of solutions to system (57) is guaranteed by the following theorem.
Then, there exists a weak solution (ρ ε , u ε ) to the system (57) with the boundary conditions (6), the pressure given by (4) and G given by (3). More precisely, the following norms on ρ ε and u ε are bounded uniformly in ε: and ρ ε , u ε satisfy the equations (57) in the sense of distributions.
Proof. The solution to system (57) can be constructed using the invading domains approach described in [22,Chapter 7]. This means to find the solution to (57) on a bounded domain Ω R = B(0, R) first and then to let R → ∞. To prove that (57) has a weak solution on Ω R , we need to supplement the system with Dirichlet boundary conditions for u ε and the zero Neumann boundary condition for ρ ε . The weak solutions to such problem can be constructed by the Faedo-Galerkin discretization of the momentum equation (57b) and the fixed point argument. The details of the last two steps are only slight modification of the procedure from [22] as all the additional terms related to G(p ε ) are of lower order and the basic a-priori estimates are still valid. Saying this, let us recall that at the level of Faedo-Galerkin approximation u ε is a suitable test function for the momentum equation and the continuity equation is satisfied pointwisely. Therefore, the energy estimate can be justified rigorously and it implies the following uniform in ε bounds. Lemma 6.2 Under assumptions (5) and (3), let T > 0 and ε > 0 be fixed, then there exists a nonnegative constant C (uniform in ε) such that the weak solution (ρ ε , u ε ) of Theorem 6.1 satisfies, for all t ∈ [0, T ], with E ε (t) and J ε (t) defined in (13), (14).
Proof. The proof of this fact follows exactly the proof of the energy estimate (17). The extra term in the momentum form ε∇ρ ε · ∇u ε allows to cancel the extra term coming from multiplication of the continuity equation by |uε| 2 2 .
We can also easily check that the estimate of the pressure from Lemma 3.3 is valid. Indeed, multiplying (57a) by γρ γ−1 ε , we deduce the equation for the pressure (60) Lemma 6.3 Let γ ≥ 2 and let the initial conditions satisfy (5). Then there exists a positive constant C such that uniformly with respect to ε we have Moreover, uniformly with respect to ε we have Proof. The proof of the first estimate (61) follows directly by an integration of (60) over Ω R and by letting R → ∞. The proof of the first part in estimate (62) follows directly by integration of (57a) over the space. To prove the second bound in (62), we multiply the continuity equation (57a) by ρ ε . Integrating by parts we obtain The last two terms can be bounded using (59) and (61), on account of the fact that γ ≥ 2.
With these estimates at hand, the proof of Theorem 6.1 is complete.
6.2 Passage to the limit ε → 0 Existence of weak solutions to our initial system (1) is then obtained by passing to the limit ε → 0.
Theorem 6.4 Let T > 0, and γ large enough be fixed. Let the initial conditions be given by (5).
Then, there exists a weak solution (ρ, u) to the system (1) in the sense of Definition 2.1, with the boundary conditions (6), the pressure given by (4) and G given by (3).
To pass to the limit in the rest of the terms of system (57), one needs to combine the arguments from Section 5 with the compactness of the sequence approximating the density {ρ ε } ε>0 . Note, that in Section 5 we were using the property (42) which is not available for γ fixed. However, taking γ sufficiently large one can still repeat all of the steps. The important changes concern solely the compactness argument for the sequence {ρ ε } ε>0 . Then in the rest of the proof, we only explain how to modify the method presented in Section 4 to handle the extra ε-related terms and get compactness for the sequence {ρ ε } ε>0 .

Modified definition of the weights
We first modify the weight by replacing the equation (25) into complemented with the initial data w ε (t = 0) = 1. Here λ is some nonnegative constant which will be fixed later on. We establish a similar property as Proposition 4.3 for this weight.
Lemma 6.5 Let us assume that u ε is given and uniformly bounded with respect to ε in L 2 . Then, there exists a unique solution to (63). Moreover, we have (ii) If we assume moreover that the pair (ρ ε , u ε ) solves (57a) and ρ ε is uniformly bounded in Proof. (i) Since (63) is a parabolic equation with B ε nonnegative and with initial data w ε (t = 0) = 1, we have that 0 ≤ w ε (t, x) ≤ 1.
Inequality (68) is the equivalent to (28) derived in Section 4.3 for no artificial viscosity case. We estimate the new term A 4 by noticing that by definition of K h we have ∆K h ≤ C h 2 K h . Then we may bound The terms A 1 and A 2 may be estimated as before. For the term A 3 , the estimate should be adapted since the relation (23) is not valid anymore. Indeed, there is an extra term where F ε = ε(−∆) −1 (div(div(u ε ⊗ ∇ρ ε ))). Hence, we arrive at the following equivalent of (32), The second term on the right hand side may be controlled as before thanks to (33). To control the third term on the right hand side of (70), we truncate using the function φ as in Section 4.3. Since D(ρ ε u ε ) + F ε is uniformly bounded in L 2 ([0, T ] × R d ), we may write as before (see (34)), x, y) dx dy dt where the function Φ h is defined, similarily as in (35), by Φ h (t, x, y) = (ρ ε (y) − ρ ε (x)) ρ ε (x)K h * w ε (x)φ ρ ε (t, x) L φ ρ ε (t, y) L .
By definition of the truncation φ, we have that Φ h L ∞ ≤ CL 2 . In particular, it allows us to use Lemma 4.5 to bound the second term on the right hand side. Here we actually use the extension of Lemma 4.5 to the case when ρ ε satisfies the continuity equation with additional dissipation term (57a). On account of Remark 4.6 and Lemma 8.3 in [3] the resulting estimate is the same. The last term in (71), thanks to the truncation, may be bounded by From Lemma 6.2 and Lemma 6.3, we deduce that the sequence { 1 √ ε F ε } ε>0 is uniformly (with respect to ε) bounded in L 1 loc ([0, T ] × R d ). Therefore the sequence {ε −1/4 F ε } ε>0 converges to 0 strongly, and therefore is compact in L 1 loc ([0, T ] × R d ). On account of Lemma 4.2 it implies that for Thus integrating in time (70), using (33) and Lemma 4.5, we arrive at Choosing L = h −θ/(α+2) , we deduce that there exists θ 0 = αθ α+2 such that This estimate is the equivalent to estimate (37).

Removing the weights
The last step consists in removing the weight. Introducing ω η = {x : K h * w ε ≤ η}, we use the same idea as in Section 4.4 to remove the weight w ε , using Lemma 6.5, and arrive at a similar estimate as (39)