Charles Elbar
email: charles.elbar@sorbonne-universite.fr

Alexandre Poulain
email: alexandre.poulain@univ-lille.fr

Analysis and numerical simulation of a generalized compressible Cahn-Hilliard-Navier-Stokes model with friction effects

Keywords: 2010 Mathematics Subject Classification. 35B40, 35B45, 35G20, 35Q35, 35Q92, 65M08 Cahn-Hilliard equation, Navier-Stokes equation, Asymptotic analysis, Mathematical modeling, Numerical simulations, Scalar Auxiliary Variable method

come

Analysis and numerical simulation of a generalized compressible Cahn-Hilliard-Navier-Stokes model with friction effects

Introduction

We derive, analyze and simulate numerically the generalized compressible Navier-Stokes-Cahn-Hilliard variant (G-NSCH in short)

∂ρ ∂t + div (ρv) = 0, (1.1)
∂(ρc) ∂t + div (ρcv) = div (b(c)∇µ) + F c , (1.2)
ρµ = -γ∆c + ρ ∂ψ 0 ∂c , (1.3)
∂(ρv) ∂t + div (ρv ⊗ v) = -∇p + γdiv ∇c ⊗ ∇c - 1 2 |∇c| 2 1 + div ν(c) ∇v + ∇v T - 2 3 ∇ (ν(c) (div (v))) -κ(ρ, c)v, (1.4)
stated in (0, T ) × Ω, where T > 0 is finite time, and

Ω ⊂ R d (d = 1, 2, 3
) is an open bounded domain with a smooth boundary ∂Ω.

Interested in the modeling of invasive growth of tumors in healthy tissues, we motivate the different terms of the model with this biological application in mind. System (1.1)-(1.4) models the motion of a diphasic fluid composed of two immiscible components,i.e. the cells of the two different types, in a porous matrix and comprises viscosity effects, surface tension, and friction on the rigid fibers constituting the medium. In System (1.1)- (1.4), ρ is the total density of the mixture (i.e. the sum of the two partial densities), c is the relative mass fraction of one component (e.g. the cancer cells), v is the mass averaged total velocity, µ is called the chemical potential, p is the pressure. The coefficient γ is related to the surface tension and is equal to the square of the width of the diffuse interface existing between the two populations. The friction coefficient κ(•) is a monotone increasing function of the density and takes into account the possible difference of friction strength between the two populations. We use this friction term to model possible adhesive effects on the extracellular matrix (ECM in short). The coefficient ν(•) represents the viscosity of the mixture and again possible differences of viscosities could be considered for the two populations. The function ψ 0 represents the separation of the two components of the mixture and phenomenologically models the behavior of cells (i.e. cells tend to form aggregates of the same cell type). The function F c (•) accounts for the possible proliferation and death of cells. The non-negative function b(•) models the mobility of cells and is assumed to be doubly degenerated to, again, correspond to the behavior of cells. This latter assumption models the probability for a cell of any of the two populations to find an available neighboring spot to which it can move. More details about the general assumptions and precise forms of the different functions will be given in the next sections.

The motivation of our model stands from the modeling of tumor progression and invasion in healthy tissues. Indeed, our model can be viewed as a representation of a proliferating population of cells (i.e. the tumor cells) in a domain filled with a non-proliferating population (i.e. the healthy cells). However, we emphasize that this article concerns the analysis and numerical simulation of the general G-NSCH model (1.1)- (1.4). This latter comprises effects that are negligible in biological situations, e.g. inertia effects. Since the general model is of interest to material sciences, physics and fluid mechanics, we focus propose here an analysis of the model and a structure preserving numerical scheme for the G-NSCH model while keeping in mind our initial application, i.e. invasive tumor growth modeling. We also emphasize that the G-NSCH is the basis of a reduced model that takes into account only biologically relevant physical effects that play a role in invasive tumor growth. Therefore, this work has to be seen as the first part. The second will concern numerical simulations, and sensitivity analysis of the reduced model, presented here in Appendix B as Problem 2, that will rely heavily on the present work.

Literature review The motion of a binary mixture of two immiscible and compressible fluids can be described by the Navier-Stokes equation coupled to the Cahn-Hilliard model. The well-known incompressible variant of the compressible NSCH model has been denominated model H (see e.g. [START_REF] Gurtin | Two-phase binary fluids and immiscible fluids described by an order parameter[END_REF][START_REF] Hohenberg | Theory of dynamic critical phenomena[END_REF]). Model H has been proposed to represent viscous fluid flow in an incompressible binary mixture under phase separation. This model assumes matching densities, i.e. ρ 1 = ρ 2 and, hence, constant total density ρ. To consider non-matching densities, Lowengrub and Truskinovsky [START_REF] Lowengrub | Quasi-incompressible Cahn-Hilliard fluids and topological transitions[END_REF] proposed the compressible Navier-Stokes Cahn-Hilliard model (NSCH model in short). Expanding the divergence term in the mass balance equation, the authors found a relation denoting the quasi-compressible nature of the fluid. Concomitantly, Anderson, FcFadden, and Wheeler [START_REF] Anderson | Diffuse-interface methods in fluid mechanics[END_REF] proposed a similar system and we use this latter in the present work.

We also remark that a very recent work [START_REF] Eikelder | A unified framework for Navier-Stokes Cahn-Hilliard models with non-matching densities[END_REF] proposed a unified framework for the incompressible NSCH system and shows that the different NSCH models found in the literature only differ from their general model by specific constitutive hypotheses.

Under some simplifying assumptions compared to the system proposed in [START_REF] Lowengrub | Quasi-incompressible Cahn-Hilliard fluids and topological transitions[END_REF] but being closer to the system in [START_REF] Anderson | Diffuse-interface methods in fluid mechanics[END_REF], the analysis of the compressible NSCH model with no-flux boundary conditions has been realized by Abels and Feireisl [START_REF] Abels | On a diffuse interface model for a two-phase flow of compressible viscous fluids[END_REF]. Their analysis requires simplifying the model proposed in [START_REF] Lowengrub | Quasi-incompressible Cahn-Hilliard fluids and topological transitions[END_REF] to avoid zones with zero density which would make this analysis a lot more difficult since the control from certain estimates would be lost. In another article, for the same system, Abels proved the existence of strong solutions for short times [START_REF]Strong well-posedness of a diffuse interface model for a viscous, quasi-incompressible two-phase flow[END_REF]. Considering the same assumptions and dynamic boundary conditions, Cherfils et al. [START_REF] Cherfils | The compressible Navier-Stokes-Cahn-Hilliard equations with dynamic boundary conditions[END_REF] proved the well-posedness of the compressible NSCH model with these special boundary conditions. These latter allow to model the interaction of the fluid components and the walls of the domain.

Results on the analysis of the incompressible variant of the NSCH model, i.e. the model H, are numerous and we here mention only a few of them since a complete review would be out of the scope of the present article. With a non-degenerate mobility coefficient (b(c) in our notation) and a physically relevant choice of potential, the well-posedness and regularity analysis of model H has been performed by Abels [START_REF] Abels | On a diffuse interface model for two-phase flows of viscous, incompressible fluids with matched densities[END_REF] using tools both from the analysis of Navier-Stokes model and the Cahn-Hilliard model. It is worth mentioning that the non-degeneracy of the mobility coefficient leads to non-physical effects, i.e. Ostwald ripening effects (see [START_REF] Abels | Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities[END_REF]). For this reason, Abels, Depner and Garcke studied model H with a degenerate mobility [START_REF] Abels | On an incompressible Navier-Stokes/Cahn-Hilliard system with degenerate mobility[END_REF]. Their analysis relies on a regularization of the mobility and singular potential into, respectively, a non-degenerate and non-singular potential. Then, suitable a-priori estimates uniform in the regularization parameter allow to pass to the limit in the regularization and show the existence of weak solutions to the degenerate model H. We now focus on the Cahn-Hilliard equation alone and its use for the modelling of tumors. The Cahn-Hilliard equation has been initially used to represent the phase separation in binary mixtures and has been applied to the spinodal decomposition of binary alloys under a sudden cooling [START_REF] Cahn | On spinodal decomposition[END_REF][START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF]. The model represents the two phases of the fluids as continua separated by a diffuse interface. This equation has been used later in many different applications and we do not intend here to give an overview of all these. However, we refer the reader interested in the topic to the presentation of the Cahn-Hilliard equation and its applications to the review book [START_REF] Miranville | The Cahn-Hilliard equation[END_REF]. We are interested here in the application of the Cahn-Hilliard framework to tumor modelling (see e.g. [START_REF] Lowengrub | Analysis of a mixture model of tumor growth[END_REF][START_REF]Analysis of a mixture model of tumor growth[END_REF]). Latter, different variants of the Cahn-Hilliard model appeared: e.g. , without giving a complete overview again, its coupling to Darcy's law [START_REF] Garcke | A cahn-hilliard-darcy model for tumour growth with chemotaxis and active transport[END_REF], Brinkman's law [START_REF] Ebenbeck | Cahn-hilliard-brinkman systems for tumour growth[END_REF], chemotaxis [START_REF] Rocca | On a cahn-hilliard-keller-segel model with generalized logistic source describing tumor growth[END_REF]. Recently, a new variant has been used to better represent the growth and organization of tumors. The main change is to consider a single-well logarithmic degenerate potential instead of a double-well potential [START_REF] Agosti | A Cahn-Hilliard-type equation with application to tumor growth dynamics[END_REF][START_REF] Chatelain | Emergence of microstructural patterns in skin cancer: a phase separation analysis in a binary mixture[END_REF][START_REF] Perthame | Relaxation of the Cahn-Hilliard equation with singular single-well potential and degenerate mobility[END_REF]. This type of potential has been proposed in [START_REF] Ambrosi | On the closure of mass balance models for tumor growth[END_REF] to represent the action of the cells depending only on their own local density, i.e. attraction at low cell density and repulsion for large cell density representing the tendency of cells to avoid overcrowding.

The numerical simulation of Model H for binary fluids with non-matching densities has been the subject of numerous works (see e.g. [START_REF] Hosseini | Isogeometric analysis of the Navier-Stokes-Cahn-Hilliard equations with application to incompressible two-phase flows[END_REF] and references therein). However, in part due to its complexity, the numerical simulation of the compressible NSCH system has been less explored. A C 0 finite element numerical scheme for a variant of the quasi-compressible NSCH model proposed in [START_REF] Lowengrub | Quasi-incompressible Cahn-Hilliard fluids and topological transitions[END_REF] has been proposed in [START_REF] Guo | A numerical method for the quasi-incompressible Cahn-Hilliard-Navier-Stokes equations for variable density flows with a discrete energy law[END_REF]. Around the same time, Giesselmann and Pryer [START_REF] Aki | A quasi-incompressible diffuse interface model with phase transition[END_REF][START_REF] Giesselmann | Energy consistent discontinuous Galerkin methods for a quasi-incompressible diffuse two phase flow model[END_REF] designed a discontinuous Galerkin finite element scheme to simulate the quasi-incompressible NSCH system which preserves the total mass and the energy dissipation. A numerical method has also been proposed in [START_REF] He | Numerical study of compressible Navier-Stokes-Cahn-Hilliard system[END_REF] in the case of constant mobility b(c) and smooth polynomial potential ψ(c). Furthermore, the system simulated in [START_REF] He | Numerical study of compressible Navier-Stokes-Cahn-Hilliard system[END_REF] is a simplification of the compressible NSCH system since the pressure does not appear in the definition of the chemical potential µ in their system.

The previous works we presented for the simulation of the compressible or quasi-compressible NSCH systems deal with constant mobility combined with a smooth polynomial potential. We aim to simulate the compressible NSCH model with choices of mobility and potential relevant for biology (but also relevant for material sciences and fluid mechanics), i.e. degenerate mobil-ity combined with a logarithmic potential. We now review briefly some relevant discretization method for the Cahn-Hilliard equation alone with degenerate mobility and singular potentials. Considering a degenerate mobility and a double-well logarithmic potential, we mention the work of Barrett, Blowey and Garcke [START_REF] Barrett | Finite element approximation of the Cahn-Hilliard equation with degenerate mobility[END_REF]. In this article the authors proposed a finite element scheme with a variational inequality to preserve the bounds of the solution. Based on these ideas, Agosti et al. [START_REF] Agosti | A Cahn-Hilliard-type equation with application to tumor growth dynamics[END_REF] proposed a similar finite element scheme for the case of single-logarithmic potential. The difficulty in this latter case lies in the fact that the degeneracy and the singularity sets do not coincide and negative considering an order parameter that must remain within the bounds [0, 1) negative solutions can appear if a standard discretization method is used. The method proposed in [START_REF] Agosti | A Cahn-Hilliard-type equation with application to tumor growth dynamics[END_REF] solves this issue but does not preserve the exact mass. In a more recent work, Agosti [START_REF] Agosti | Discontinuous Galerkin finite element discretization of a degenerate Cahn-Hilliard equation with a single-well potential[END_REF] proposed a discontinuous Galerkin finite element scheme that preserves the bounds [0, 1) and preserves the exact mass. However, the main drawback of the previously mentioned methods is that they are computationally expensive and solve a strongly coupled nonlinear system and use iterative algorithms. Since the Cahn-Hilliard equation is a gradient flow (see e.g. [START_REF] Lisini | Cahn-Hilliard and thin film equations with nonlinear mobility as gradient flows in weighted-Wasserstein metrics[END_REF]), a structure-preserving linear scheme can be constructed using the Scalar Auxiliary Variable (SAV in short) method [START_REF] Shen | The scalar auxiliary variable (SAV) approach for gradient flows[END_REF]. This has been successfully used in [START_REF] Huang | Bound/positivity preserving and unconditionally stable schemes for a class of fourth order nonlinear equations[END_REF]. In this latter work, the scheme is structure-preserving from the use of a scalar variable that represents the discrete energy, and an additional equation is solved to ensure dissipation at the discrete level. The bounds of the order parameter are ensured using a transformation that maps R to the physical relevant interval ((0, 1) in the case of a double-well potential). To the best of our knowledge, the SAV method has not been applied to the compressible NSCH system.

Objectives of our work

The first objective of our work is to study the well-posedness of the G-NSCH model under some simplifying assumptions (i.e. smooth potential and positive mobility).

The second objective is the design of an efficient and structure-preserving numerical scheme for the G-NSCH model with singular double-well potential and degenerate mobility. The third focus of the present work concerns the rigorous derivation of the G-NSCH model that is presented in the Appendix.

Outline of the paper Section 2 presents the notation, functional spaces and assumptions we use in our work for the analytical part but also for the numerical part. Section 3 deals with the proof of the existence of weak solutions for the G-NSCH system (1.1)-(1.4) under simplifying assumptions. A structure preserving numerical scheme based on the SAV method is then proposed in Section 4 and some numerical results are presented in Section 5. Our model's equations come from a thermodynamically consistent derivation of the compressible Navier-Stokes-Cahn-Hilliard model including friction effects and source terms. The derivation is described in Appendix A. From a general model, we propose in Appendix B two reductions:

The G-NSCH studied and simulated in the present work and one biologically relevant reduction that will be the focus of a forthcoming work.

General assumptions, notations and functional setting

The equations are set in a domain

Ω T = Ω × (0, T )
with Ω an open and bounded subset of R 3 . We assume that the boundary ∂Ω is sufficiently smooth. We indicate the usual Lebesgue and Sobolev spaces by respectively L p (Ω), W m,p (Ω) with H m (Ω) := W m,2 (Ω), where 1 ≤ p ≤ +∞ and m ∈ N. For q ∈ [1, +∞], we indicate the Bochner spaces by L q (0, T ; X) (where X is a Banach space). Finally, C denotes a generic constant that appears in inequalities and whose value can change from one line to another. This constant can depend on various parameters unless specified otherwise.

Assumptions on functionals

We divide the assumptions on the different terms appearing in system (1.1)-(1.4) into two parts: analytical and numerical assumptions. Indeed we are not able to prove the existence of weak solutions in the general setting used in the numerical simulations. For instance, the case of the usual logarithmic double-well potential in the Cahn-Hilliard equation is not treated but can be implemented in our numerical scheme. However, we can analyze our system with a polynomial approximation of the double well. We also consider non-degenerate mobilities to obtain estimates on the chemical potential µ directly. The case of degenerate mobility, see for instance [START_REF] Elliott | On the Cahn-Hilliard equation with degenerate mobility[END_REF], seems unavailable as we do not have anymore the classical "entropy" estimates of the Cahn-Hilliard equation that provide bound on second-order derivatives of the mass fraction c.

Framework for numerical simulations We assume that the viscosity ν(c) and permeability κ(ρ, c) coefficients are smooth non-negative functions of the mass fraction c. The mobility is a non-negative function of the order parameter (mass fraction) c. Hence, we assume that

b ∈ C 1 ([0, 1]; R + ), and b(c) ≥ 0 for 0 ≤ c ≤ 1. (2.1)
In agreement with the literature (see e.g [START_REF] Cherfils | The compressible Navier-Stokes-Cahn-Hilliard equations with dynamic boundary conditions[END_REF]), the homogeneous free energy ψ 0 (ρ, c) is assumed to be of the form

ψ 0 (ρ, c) = ψ e (ρ) + ψ mix (ρ, c), (2.2)
with ψ mix (ρ, c) = H(c) log ρ + Q(c) and Q(c) is a double-well (or single-well) potential. Then, using the constitutive relation for the pressure, we have

p(ρ, c) = ρ 2 ∂ψ 0 ∂ρ = p e (ρ) + ρH(c), (2.3)
where p e = ρ 2 ψ ′ e (ρ) and is assumed to satisfy

p 1 ρ a-1 -p 2 ≤ p ′ e (ρ) ≤ p 3 (1 + ρ a-1
), for a > 3/2, p 1 , p 2 , p 3 > 0.

(2.4)

We assume that the source term F c (that can depend on the mass fraction and the density) is bounded,

|F c (ρ, c)| + F c (ρ, c) ρ ≤ C, ∀(ρ, c) ∈ R 2 .
(2.5)

Remark 2.1 (Double-well logarithmic potential). In the present work, we aim to use a doublewell logarithmic potential in the definition of the mixing potential. A relevant example of potential is

ψ mix = 1 2 (α 1 (1 -c) log(ρ(1 -c)) + α 2 c log(ρc)) - θ 2 (c - 1 2 ) 2 + k. (2.6)
This potential gives

H(c) = α 1 (1 -c) + α 2 c, Q(c) = 1 2 (α 1 (1 -c) log(1 -c) + α 2 c log(c)) - θ 2 (c - 1 2 ) 2 + k,
where θ > 1 and k is an arbitrary constant.

Additional assumptions for the existence of weak solutions Concerning the existence of weak solutions, we need to strengthen our assumptions. The viscosity coefficient ν(c) is assumed to be bounded from below by a positive constant and the friction coefficient κ(ρ, c) is assumed to be nonnegative. Moreover both ν(c) and κ(ρ, c) are two functions bounded in L 2 (0, T ; L 2 (Ω)) whenever c is bounded in L ∞ (0, T ; H 1 (Ω)) and ρ is smooth (for instance C(0, T ; C 2 (Ω)). We consider a > 2 the exponent of the pressure law. In the numerical simulations, we take degenerate mobilities of the form b(c) = c(1 -c) α . However, in the analysis, we consider a non-degenerate mobility by truncating the previous mobility. For instance, using a small parameter 0 < ε b << 1, we approximate the mobility b(•) by

b ε b (c) =      b(1 -ε b ), if c ≥ 1 -ε b , b(ε b ), if c ≤ ε b , b(c), otherwise,
and consider the case of a fixed ε b . We obtain that b ∈ C 1 (R; R + ), and b(c) ≥ C > 0 ∀c ∈ R.

(2.7)

Concerning the functionals appearing in the definition of the free energy ψ 0 we assume that H and H ′ are bounded and that Q is a polynomial approximation of the double well potential. More precisely we take

H 1 ≤ H ′ (c), H(c) ≤ H 2 , c ∈ R, H 1 , H 2 > 0, Q(c) = 1 4 c 2 (1 -c) 2 .
(2.8)

The case of the double-well logarithmic potential has not been tackled yet even though this is the main motivation for the decomposition of ψ mix as in the works [START_REF] Abels | On a diffuse interface model for a two-phase flow of compressible viscous fluids[END_REF] and [START_REF] Cherfils | The compressible Navier-Stokes-Cahn-Hilliard equations with dynamic boundary conditions[END_REF]. Also, to make the computations simpler, we assume that

• a > 6 where a is the pressure exponent,

• ψ e (ρ) = ρ a-1 a-1 and therefore p e (ρ) = ρ a . These two assumptions are not necessary but simplify the analysis. We refer for instance to [START_REF] Abels | On a diffuse interface model for a two-phase flow of compressible viscous fluids[END_REF][START_REF] Feireisl | Dynamics of viscous compressible fluids[END_REF] for the more general setting. For instance, the condition a > 6 is used to not introduce another parameter in the approximating scheme which would make the article even longer. Note that the assumptions on ψ 0 imply in particular the following lemma which is essential to obtain estimates on the energy dissipation: Lemma 2.2. There exists a constant C such that

ρ ∂ψ 0 ∂c ≤ Cρψ 0 + C.
Its proof uses the assumption on H and the fact that for c large,

Q ′ (c) ≈ c 3 ≤ c 4 + 1 ≈ Q(c) + 1.

Existence of weak solutions

We now turn to the proof of the existence of weak solutions for the G-NSCH model (1.1)-(1.4) subjected to boundary conditions

v = ∂c ∂n = b(c) ∂µ ∂n = 0 on ∂Ω, (3.1)
and initial conditions

ρ(0, x) = ρ 0 ≥ 0 ∈ L a (Ω), c(0, x) = c 0 ∈ H 1 (Ω) ρ 0 v(0, x) = m 0 , with |m 0 | 2 ρ 0 ∈ L 1 (Ω). (3.2)
Also, we suppose ρ 0 ̸ = 0. The proof of the result is quite long and technical. Therefore, when necessary and for the sake of clarity, we omit some proofs and give instead appropriate references.

Outline of the analysis For readability reasons, we here present the plan we use for the analysis of the G-NSCH model. We first start with the analysis of a "truncated" version of G-NSCH model in the sense that the double-well is truncated for large values of c with a parameter ε Q . Then, for this fixed truncation, we prove the existence of weak solutions using the ideas of [START_REF] Abels | On a diffuse interface model for a two-phase flow of compressible viscous fluids[END_REF][START_REF] Cherfils | The compressible Navier-Stokes-Cahn-Hilliard equations with dynamic boundary conditions[END_REF][START_REF] Feireisl | Dynamics of viscous compressible fluids[END_REF][START_REF] Lions | Compressible models[END_REF]. Then, we pass to the limit ε Q → 0. Namely, recalling that

Q(c) = 1 4 c 2 (1 -c) 2 we first consider Q ε Q (c) a smooth truncated approximation of Q which satisfies |Q ε Q |, |Q ′ ε Q |, |Q ′′ ε Q | ≤ C 1 ε Q . (3.3)
In the first subsections, we work with the regularized problem and we drop the ε Q notation. We will use the ε Q notation when we pass to the limit, and for the moment we benefit from the properties of the regularization.

Energy estimates

The G-NSCH system comes with an energy structure which is useful to obtain first a priori estimates.

Proposition 3.1. Smooth solutions of the system (1.1)-(1.4) satisfy the following energy relation

d dt E + D = Ω µF c dx, (3.4)
where E is the energy, and D is the dissipation defined as

E = Ω ρ |v| 2 2 + ρψ 0 + γ 2 |∇c| 2 dx, (3.5)
D = Ω ν(c) 2 ∇v + ∇v T - 2 3 div(v)I 2 + b(c)|∇µ| 2 + κ(ρ, c)|v| 2 dx. (3.6)
This yields a priori estimates on the solution i.e. there exists a positive constant C such that

E(t) + t 0 D(s) ds ≤ C + CE(0).
Note that the energy is bounded from below since ρ log ρH(c) is bounded from below with (2.8). Also, the purpose of the assumptions ν(c) and b(c) bounded from below by a positive constant become clear, they are crucial to obtain estimates on the H 1 (Ω) norm of µ and v.

Proof. We recall the formula

∇c∆c = div(∇c ⊗ ∇c) - 1 2 ∇|∇c| 2 . (3.7)
We denote by T the tensor ν(c)(∇v + ∇v T -2 3 div(v)I). Then we multiply Equation (1.1) by

|v| 2
2 and sum it with the scalar product of Equation (1.4) with v. We obtain

∂ ∂t ρ |v| 2 2 + div 1 2 ρ|v| 2 v + p(ρ, c)v -T • v + T : ∇v + κ(ρ, c)v 2 = p(ρ, c)div(v) + γdiv( 1 2 |∇c| 2 I -(∇c ⊗ ∇c)) • v,
which is equivalent to

∂ ∂t ρ |v| 2 2 + div 1 2 ρ|v| 2 v + p(ρ, c)v -T • v + T : ∇v + κ(ρ, c)v 2 = p(ρ, c)divv -γ∆c∇c • v. (3.8)
Then, we multiply Equation (1.2) by µ and obtain using also (1.1)

ρµ(∂ t c + v • ∇c) = div(b(c)∇µ)µ + µF c .
And, using (1.3) we obtain

ρ ∂ψ 0 ∂c (∂ t c + v • ∇c) = div(b(c)∇µ)µ + γ∆c(∂ t c + v • ∇c) + µF c .
The previous equation can be rewritten using the chain rule as

∂ t (ρψ 0 ) + div(ρψ 0 v) -ψ 0 (∂ t ρ + div(ρv)) -ρ ∂ψ 0 ∂ρ (∂ t ρ + v • ∇ρ) = div(b(c)∇µ)µ + γ∆c(∂ t c + v • ∇c) + µF c .
We have

ρ ∂ψ 0 ∂ρ (∂ t ρ + v • ∇ρ) = ρ ∂ψ 0 ∂ρ (-ρdiv(v)) = -pdiv(v) (see Equation (2.
3) for the definition of the pressure). Moreover, we know that

∆c∂ t c = div(∂ t c∇c) -∂ t |∇c| 2 2
and, hence,

∂ t (ρψ 0 ) + div(ρψ 0 v) + pdiv(v) = div(b(c)∇µ)µ + γ div(∂ t c∇c) -∂ t |∇c| 2 2 + ∆cv • ∇c + µF c . (3.9)
Summing (3.8) and (3.9) we obtain

∂ ∂t ρ |v| 2 2 + ρψ 0 + γ 2 |∇c| 2 +div ρψ 0 v + 1 2 ρ|v| 2 v + p(ρ, c)v -T : v -γ∂ t c∇c -div(b(c)∇µ)µ + T : ∇v + κ(ρ, c)|v| 2 = µF c .

Now we use the fact that

T : ∇v = ν(c) 2 ∇v + ∇v T - 2 3 div(v)I 2 .
(3.10)

Integrating in space and using the boundary conditions (3.1) ends the proof of the first part of the proposition. To prove the second part, we integrate the equation in time and control the right-hand side. Indeed, due to the assumption on the source term (2.5), we have

t 0 Ω µF c dx dt ≤ C t 0 Ω |µ|.
We want to use Lemma 3.6 to control the L 1 norm of µ. Integrating the equations on ρ to obtain Ω ρ dx = Ω ρ 0 dx > M 0 we satisfy the first assumption of the lemma. For the second, we notice that we can consider a variant of this lemma such that instead of asking ρ to be in L 1+ε we have the inequality

u - 1 Ω Ω ρu L 2 ≤ C∥∇u∥ L 2 + ∥ρ∥ L 1+ε .
Using Young's inequality, the fact that in the energy ρψ 0 contains a term of the form ρ a+1 we obtain for C small enough

t 0 Ω |µ| dx ≤ C + C t 0 Ω |µ| 2 dx ≤ C + CE(t) + inf c b(c) 2 Ω |∇µ| 2 dx + C Ω ρµ dx .
Since the energy dissipation controls the third term of the right-hand side, it remains to control the last term of the right-hand side. We recall that ρµ = ρ ∂ψ 0 ∂c -γ∆c. Using the Neumann boundary conditions on c, it remains to control Ω ρ ∂ψ 0 ∂c . Using Lemma 2.2, we obtain

Ω ρ ∂ψ 0 ∂c dx ≤ C + CE(t).
We conclude using Gronwall's lemma.

Existence of weak solutions for fixed ε Q

The weak solutions of system (1.1)-(1.4) are defined as follows Definition 3.2. We say that (ρ, v, c, µ) is a weak of system (1.1)-(1.4) provided:

• ρ ≥ 0 and we have the regularity

ρ ∈ L ∞ (0, T ; L a (Ω)), v ∈ L 2 (0, T ; H 1 0 (R 3 )), √ ρv ∈ L ∞ (0, T ; L 2 (Ω; R 3 )), T : ∇v ∈ L 1 (0, T ; L 1 (Ω)), c ∈ L ∞ (0, T ; H 1 (Ω)), µ ∈ L 2 (0, T ; H 1 (Ω)).
• Equations (1.1)-(1.4) are satisfied in the distributional sense.

• The initial conditions (3.2) are satisfied a.e. in Ω.

• The boundary conditions (3.1) are satisfied.

In order to prove the existence of weak solutions, we use an approximating scheme with a small parameter ε > 0 borrowing the idea from [START_REF] Feireisl | Dynamics of viscous compressible fluids[END_REF][START_REF]Compressible models[END_REF]. More precisely, let X n = span{η i } i=1,...,n be the set of the first n vectors of a basis of

H 1 0 (Ω; R 3 ) such that X n ⊂ C 2 (Ω; R 3
). We consider the following problem for (ρ, v n , c) with v n ∈ X n (with coordinates depending on time):

∂ t ρ + div(ρv n ) = ε∆ρ, (3.11)
and for every η ∈ X n ,

Ω ρv n (t) • η dx - Ω m 0 • η dx - t 0 Ω ρv n ⊗ v n : ∇η dx ds - t 0 Ω p(ρ, c)div(η) dx ds + ε t 0 Ω (∇v n ∇ρ) • η dx ds + t 0 Ω T : ∇η dx ds + γ t 0 Ω ( 1 2 |∇c| 2 I -(∇c ⊗ ∇c)) : ∇η dx ds + Ω t 0 κ(ρ, c)v n • η dx ds = 0. (3.12)
And for the equation on the mass fraction

∂ t c + v n • ∇c = 1 ρ div(b(c)∇µ) + F c ρ , µ = ∂ψ 0 ∂c -γ ∆c ρ . (3.13)
We consider Neumann boundary conditions

∇ρ • n = b(c)∇µ • n = ∇c • n = 0 on ∂Ω, (3.14)
and the Dirichlet boundary condition for v n is included in the definition of X n . Finally, we consider the initial conditions

ρ(0, •) = ρ 0,ε > 0, c(0, •) = c 0,ε , ρv n (0, •) = m 0 , (3.15)
where ρ 0,ε , c 0,ε satisfy the Neumann boundary conditions and they are smooth approximations of ρ 0 , c 0 (when ε → 0). We now comment on the scheme used above and detail the strategy of the proof. We add the artificial diffusion in (3.11) with the parameter ε > 0. Here, v n is fixed and we can conclude the existence of classical solutions to (3.11) which are positive since the initial condition is positive (and using maximum principle). Using this positivity, we conclude the existence of a strong solution to Equation (3.13) which is in fact a fourth-order parabolic equation. Having obtained c, we focus on Equation (3.12) and we prove existence for a small time with Schauder's fixed point theorem. Note the presence of the additional term ε (∇v n ∇ρ) • η which is useful to cancel energy terms introduced by ε∆ρ in (3.11). Having obtained existence on a short time interval we compute the energy of the system and obtain global existence. Then, we pass to the limit n → ∞. It remains to send ε and ε Q to 0 and obtain solutions of system (1.1)-(1.4). We first turn our attention to Equation (3.11). From [START_REF] Feireisl | Dynamics of viscous compressible fluids[END_REF], we obtain the following proposition, and lemma Proposition 3.3. Let Ω ⊂ R 3 be a bounded domain of class C 2+β for some β > 0. For a fixed v n ∈ X n , there exists a unique solution to Equation (3.11) with Neumann boundary conditions (3.14) and initial data conditions (3.15). Furthermore, the mapping

v n → ρ[v n ],
that assigns to any v n ∈ X n the unique solution of (3.11), takes bounded sets in the space C(0, T ; C 2 0 (Ω, R d )) into bounded sets in the space

V := {∂ t ρ ∈ C(0, T ; C β (Ω)), ρ ∈ C(0, T ; C 2+β (Ω))}.
Lemma 3.4. The solutions of (3.11) satisfy

( inf x∈Ω ρ(0, x)) exp - t 0 ∥div v n (s)∥ L ∞ (Ω) ds ≤ ρ(t, x) ≤ (sup x∈Ω ρ(0, x)) exp t 0 ∥div v n (s)∥ L ∞ (Ω) ds , for all t ∈ [0, T ] and x ∈ Ω.
Using the latter lemma, if the velocity field is in W 1,∞ , the density is bounded from below by a positive constant (provided the initial condition is positive). We now focus on Equation (3.13).

W := {c ∈ L ∞ (0, T ; H 1 (Ω)) ∩ L 2 (0, T ; H 3 (Ω))}. (3.16)
The existence of a strong solution is based on the remark that the highest order term of this equation is -γ b(c) ρ ∆ 2 c. Using b(c), ρ ≥ C > 0 we obtain a fourth-order parabolic equation with smooth coefficients and with zero Neumann boundary conditions. Therefore, we can admit the existence of a strong solution and we focus on the estimates (3.16). In the proof, we need the following two lemmas Lemma 3.6 (Lemma 3.2 in [START_REF] Feireisl | Dynamics of viscous compressible fluids[END_REF]). Let Ω ∈ R 3 be a bounded Lipschitz domain and let M 0 > 0, K > 0. Assume that ρ is a nonnegative function such that

0 < M 0 ≤ Ω ρ dx, Ω ρ a dx ≤ K, with a > 1.
Then, there exists a positive constant C = C(M 0 , K, a) such that the inequality

u - 1 |Ω| Ω ρu L 2 (Ω;R 3 ) ≤ C∥∇u∥ L 2 (Ω;R 3×3 ) ,
holds for any u ∈ W 1,2 (Ω; R 3 ).

Lemma 3.7 (Theorem 10.17 in [START_REF] Feireisl | Singular limits in thermodynamics of viscous fluids[END_REF]). Let Ω ⊂ R 3 be a bounded Lipschitz domain, and let

1 < p < +∞, M 0 > 0, K > 0, a > 1.
Then there exists a postive constant C = C(p, M 0 , K, a) such that the inequality

∥u∥ W 1,p (Ω;R 3 ) ≤ C ∥∇u + ∇ T u - 2 3 divuI∥ L p (Ω;R 3×3 ) + Ω ρ|u| dx ,
holds for any u ∈ W 1,p (Ω; R 3 ) and any non-negative function ρ such that

0 < M 0 ≤ Ω ρ dx, Ω ρ a dx ≤ K.
Proof of Proposition 3.5. We admit the existence of solutions and focus on a priori estimates. We multiply Equation (3.13) by -∆c. Using the boundary conditions and integrating in space yields

∂ t Ω |∇c| 2 2 dx + γ Ω b(c) ∇ ∆c ρ 2 dx = Ω 1 2 div(v n )|∇c| 2 -∇v n : ∇c ⊗ ∇c dx + Ω b(c)∇ ∂ψ 0 ∂c • ∇ ∆c ρ dx - Ω F c ρ ∆c.
Here, we have also used the formula (3.7). We use the

L ∞ bounds on v n , div(v n ), b(c), ρ the fact that Fc ρ is also bounded in L ∞ , properties on ∂ c ψ 0 (3.
3), and obtain

∂ t Ω |∇c| 2 2 dx + γ Ω b(c) ∇ ∆c ρ 2 dx ≤ C Ω |∇c| 2 dx + C Ω ∇ ∆c ρ dx + C Ω |∆c|.
We want to control the last term on the right-hand side. We use Lemma 3.6 with u = ∆c ρ (1, 0, 0) T and obtain, together with Neumann boundary conditions on c,

∆c ρ L 2 (Ω) ≤ C ∇ ∆c ρ L 2 (Ω;R 3 ) . (3.17)
Then, writing ∆c = ρ ∆c ρ and using the L ∞ bound on ρ,

Ω |∆c| ≤ C ∇ ∆c ρ L 2 (Ω;R 3 )
.

Finally, using Young's inequality and Gronwall's lemma, we obtain

sup t∈(0,T ) Ω |∇c| 2 dx + γ T 0 Ω b(c) ∇ ∆c ρ 2 dx ≤ C. (3.18)
With Lemma 3.6 (and integrating the equation on ρc using also the boundary conditions) we obtain the bound

c ∈ L ∞ (0, T ; H 1 (Ω)) ∩ L 2 (0, T ; H 3 (Ω)). (3.19)
Having defined ρ and c, we now solve Equation (3.12) with a fixed point argument. We define the operator

M[ρ] : X n → X * n , ⟨M[ρ]v, w⟩ := Ω ρv • wdx, v, w ∈ X n .
This operator ( [START_REF] Feireisl | Dynamics of viscous compressible fluids[END_REF]) M[ρ] is invertible, and

∥M -1 [ρ]∥ L(X * n ;Xn) ≤ 1 inf Ω ρ , ∥M -1 [ρ 1 ]-M -1 [ρ 2 ]∥ L(X * n ;Xn) ≤ C(n, ρ)∥ρ 1 -ρ 2 ∥ L 1 (Ω) , (3.20)
for any ρ 1 , ρ 2 ≥ ρ. Finally, Equation (3.12) can be reformulated as

v n (t) = M -1 [ρ(t)] m * 0 + t 0 N [v n (s), ρ(s), c(s)] ds , (3.21)
with

⟨m * 0 , η⟩ = Ω m 0 • η dx,
and

⟨N [v n , ρ, c], η⟩ = Ω ρv n ⊗ v n -T - γ 2 |∇c| 2 I + γ∇c ⊗ ∇c : ∇η + p(ρ, c)div(η) -(ε∇v n ∇ρ + κ(ρ, c)v n ) • η dx.
To prove that Equation (3.21) has a solution, we apply Schauder's fixed-point theorem in a short time interval [0, T (n)]. Then, we need uniform estimates to iterate the procedure.

Lemma 3.8 (Schauder Fixed Point Theorem). Let X be a Hausdorff topological vector space and S be a closed, bounded, convex, and non-empty subset of X. Then, any compact operator A : S → S has at least one fixed point.

With notation of the lemma 3.8, we call A the operator from Equation (3.21) and S = B(u 0,n ) the unit ball with center u 0,n in C([0, T ]; X n ), u 0,n is defined by

Ω ρ 0 u 0,n • η dx = Ω m 0 • η dx, ∀η ∈ X n .
More precisely, we consider

A : S → C([0, T ]; X n ), u → M -1 [ρ(t)] m * 0 + t 0 N [u(s), ρ(s), c(s)] ds .
Lemma 3.9. There exists a time T = T (n) small enough such that the operator A maps S into itself. Moreover, the mapping is continuous.

Proof. By definition of A and m * 0 , we need to prove that (3.20), it is sufficient to prove that there exists a final time T small enough such that

∥M -1 [ρ(t)] t 0 N (s)ds∥ C(0,T ;Xn) ≤ 1. With properties
t 0 N (s) ds C(0,T ;X * n ) ≤ inf Ω T ρ.
Note that the infimum of ρ needs to be taken over the set Ω T = (0, T ) × Ω as ρ depends on time.

But, since we only consider small times, using Lemma 3.4 we see that this infimum is bounded by below. More precisely, for every T 0 , there exists C(T 0 ) > 0 such that for every

T ≤ T 0 , inf Ω T ρ ≥ C(T 0 ). We recall that X n ⊂ C 2 (Ω; R 3 ) is finite-dimensional and we estimate t 0 Ω (ρu ⊗ u -T - γ 2 |∇c| 2 I + γ∇c ⊗ ∇c) : ∇η + p(ρ, c)div(η) -(ε∇u∇ρ + κ(ρ, c)u) • η dx ds ≤ √ T (∥η∥ Xn + ∥∇η∥ Xn )(∥ρ∥ L ∞ ∥u∥ 2 L ∞ + C∥ν(c)∥ L 2 ∥∇u∥ L ∞ + C∥∇c∥ 2 L 4 + ∥ρ∥ a L ∞ + ∥ρ∥ L ∞ ∥H(c)∥ L ∞ + ε∥u∥ Xn ∥∇ρ∥ L ∞ + ∥u∥ L ∞ ∥κ(ρ, c)∥ L 2 ).
Using assumptions of the subsection 2.1 and Propositions 3.3-3.5, we prove that all the quantities on the right-hand side are bounded, except ∥∇c∥ L 4 which needs an argument. Note that from (3.16), we deduce ∇c is bounded in L 2 (0, T ; H 2 (Ω)) ∩ L ∞ (0, T ; L 2 (Ω)) (by a constant which depends on ρ, and also on ∥u∥ L ∞ , ∥∇u∥ L ∞ ). By Sobolev embedding with d = 3, and interpolation, we obtain an L 4 (0, T ; L 4 (Ω)) bound. With the previous estimates, and for T small enough, we obtain the result.

Lemma 3.10. The image of S under A is in fact a compact subset of S. Therefore, A admits a fixed point.

Proof. We want to apply the Arzelà-Ascoli theorem to deduce the relative compactness of A(S). From the previous computation, and using the fact that X n is finite-dimensional, we can prove that A(S) is pointwise relatively compact. It remains to prove its equicontinuity. We want to estimate for t

′ ≤ t the X n norm of M -1 [ρ(t)] m * 0 + t 0 N [u(s), ρ(s), c(s)] ds - M -1 [ρ(t ′ )] m * 0 + t ′ 0 N [u(s), ρ(s), c(s)] ds .
For simplicity, we write N (s) := N [u(s), ρ(s), c(s)], and rewrite the previous difference as

M -1 [ρ(t) -ρ(t ′ )] m * 0 + t 0 N (s) ds + M -1 [ρ(t ′ )] m * 0 + t t ′ N (s) ds .
For the first term, we use (3.20) and the Hölder continuity of ρ given by Proposition 3.3. For the second term, we repeat the computations in the proof of Lemma 3.9. This ends the result.

We have the existence of a small interval [0, T (n)]. To iterate the procedure in order to prove that T (n) = T , it remains to find a bound on v n independent of T (n).

Lemma 3.11. v n is bounded in X n independently of T (n).
Proof. Note that we do not ask for a bound independent of n but only of T (n) since we use in the proof the fact that X n is finite-dimensional. The proof uses the energy structure of the equation. We differentiate Equation (3.12) in time and take η = v n as a test function. This yields

d dt Ω ρ |v n | 2 2 dx + 1 2 Ω (∂ t ρ + div(ρv n )) |v n | 2 dx - Ω p(ρ, c)div(v n ) dx - ε 2 Ω ∆ρ|v n | 2 dx + Ω T : ∇v n dx + γ Ω ( 1 2 |∇c| 2 I -(∇c ⊗ ∇c)) : ∇v n dx + Ω κ(ρ, c)|v n | 2 dx = 0. (3.22)
Here we used

Ω ∂ t (ρv n ) • v n = 1 2 d dt Ω ρ|v n | 2 dx + 1 2 Ω ∂ t ρ|v n | 2 dx, Ω div(ρv n ⊗ v n ) • v n dx = 1 2 Ω div(ρv n )|v n | 2 dx, ε Ω (∇v n ∇ρ) • v n dx = - ε 2 Ω ∆ρ|v n | 2 dx.
With (3.11), we see that (3.22) reads

d dt Ω ρ |v n | 2 2 dx - Ω p(ρ, c)div(v n ) dx + Ω T : ∇v n dx + γ Ω ( 1 2 |∇c| 2 I -(∇c ⊗ ∇c)) : ∇v n dx + Ω κ(ρ, c)|v n | 2 dx = 0. (3.23)
Now as in (3.9), we obtain with the artificial viscosity

∂ t (ρψ 0 ) + div(ρψ 0 v n ) + pdiv(v n ) -ψ 0 ε∆ρ -ερ ∂ψ 0 ∂ρ ∆ρ = div(b(c)∇µ)µ + div(∂ t c∇c) -∂ t |∇c| 2 2 + γ∆cv n • ∇c + µF c .
Integrating this equation in space, and summing with (3.23), we obtain

d dt Ω ρ |v n | 2 2 + ψ 0 + γ |∇c| 2 2 dx + ε Ω ∇ ψ 0 + ρ ∂ψ 0 ∂ρ • ∇ρ dx + Ω T : ∇v n dx + Ω b(c)|∇µ| 2 dx + Ω κ(ρ, c)|v n | 2 dx = Ω µF c dx. (3.24)
By definition of ψ 0 , we obtain

ε Ω ∇ ψ 0 + ρ ∂ψ 0 ∂ρ • ∇ρ dx = ε Ω ((a -1) + (a -1) 2 )ρ a-2 + H(c) ρ |∇ρ| 2 dx + ε Ω H ′ (c)(log(ρ) + 1) + Q ′ (c) ∇c • ∇ρ dx.
Therefore, the energy reads

d dt Ω ρ |v n | 2 2 + ψ 0 + γ |∇c| 2 2 dx + ε Ω ((a -1) + (a -1) 2 )ρ a-2 + H(c) ρ |∇ρ| 2 dx + Ω T : ∇v n dx + Ω b(c)|∇µ| 2 dx + Ω κ(ρ, c)|v n | 2 dx = Ω µF c dx -ε Ω H ′ (c)(log(ρ) + 1) + Q ′ (c) ∇c • ∇ρ dx. (3.25)
We need to prove that the right-hand side can be controlled in term of the left-hand side to obtain estimates. For the first term on the right-hand side, we treat it as in the proof of Proposition 3.1. For the second term, we know by assumption on H and Q, and the fact that (log(ρ) + 1) 2 is bounded by a constant times 1 ρ + (a -1)ρ a-2 that it can be bounded in terms of the left-hand side. Note that we used the hypothesis

|Q ′ (c)| ≤ C. This is based on the fact that Q is in fact Q ε Q so that we have |Q ′ (c)| ≤ C( 1 ε Q )
with a constant that blows up when ε Q is sent to 0. As we intend to send ε Q → 0 in the next step, it is important to notice that we can still manage to have this energy inequality since in fact the term

ε Ω Q ′ (c)∇c • ∇ρ dx can be estimated by ε 4 Ω ((a -1) + (a -1) 2 )ρ a-2 + H(c) ρ |∇ρ| 2 dx and Ω εC( 1 ε Q )|∇c| 2 dx
. Since ε will be sent to 0 before ε Q , the energy inequality will still hold independently of ε Q in the limit ε → 0. With Gronwall's lemma, and properties of the tensor T, we deduce that v n is bounded in L 2 (0, T (n); H 1 (Ω; R 3 )) independently of T (n). Since all the norms are equivalent, it is also bounded in L 1 (0, T (n); W 1,∞ (Ω, R 3 )). Therefore, we can apply the maximum principle stated in Lemma 3.4, and obtain that the density ρ is bounded from below by a constant independent of T (n). Then, using once again the energy inequality, we obtain that v n is bounded uniformly in time in L 2 (Ω; R 3 ). This procedure can be repeated for every final time T .

Finally, we are left with the following proposition Proposition 3.12. For any fixed n and T , there exists a solution (ρ, c, v n ) defined on (0, T ) (with appropriate regularity) to (3.11)-(3.13)-(3.12) subject to boundary conditions (3.14) and initial conditions (3.11). Moreover, this solution satisfies the energy dissipation inequality

E(t) + ε Ωt (a + a 2 )ρ a-1 + H(c) ρ |∇ρ| 2 dx dt + Ωt T : ∇v n dx dt + Ωt b(c)|∇µ| 2 dx dt + Ωt κ(ρ, c)|v n | 2 dx dt ≤ C + CE(0), (3.26)
where

E(t) = Ω ρ |v n | 2 2 + ψ 0 + γ |∇c| 2 2 dx,
and with a constant

C = C 1, ε ε Q .
Now, we need to find estimates, independent of n, to pass to the limit n → ∞. Since ρ and c depend on n, we write ρ n and c n from now on.

Integrating in space, using the boundary conditions, and Estimate (A7), the L 2 bound on {∇c n }, assumption 2.5 yields { Ω ρ n c n } is in L ∞ (0, T ). We deduce Estimate (A8). Estimate (A9) follows from the definition of ψ 0 and Estimate (A1). Estimate (A10) follows from Estimates (A5)-(A9) and Lemma 3.6. Estimate (A11) follows from Estimates (A2)-(A10). Estimate (A12) is a consequence of Equation (1.3), the previous estimates and interpolation. The two next estimates are a consequence of the other estimates and Sobolev embeddings. Finally, the last estimate on the pressure can be adapted from [START_REF] Cherfils | The compressible Navier-Stokes-Cahn-Hilliard equations with dynamic boundary conditions[END_REF]Subsection 2.5]. This estimate is useful when we obtain the convergence a.e. of ρ n and c n so we can obtain strong convergence of p(ρ n , c n ) in L 1 by Vitali's convergence theorem.

From [START_REF] Feireisl | Dynamics of viscous compressible fluids[END_REF], we also obtain the following Proposition Proposition 3.14. There exists r > 1 and p > 2 such that

∂ t ρ n , ∆ρ n are bounded in L r ((0, T ) × Ω), ∇ρ n is bounded in L p (0, T ; L 2 (Ω, R 3 )),
independently of n (but not independently of ε).

With all the previous bound, we can pass to the limit when n → ∞ and obtain the different equation and energy estimates in a weak formulation. Since the passage to the limit n → ∞ is simpler than the next passage ε → 0, we only detail the latter. Indeed, as n → ∞ we can obtain easily strong convergence of ρ which helps a lot in the different limits. So we assume that we can pass to the limit and that the bounds obtained in Proposition 3.13 still hold independently of ε. It remains now to send ε to 0. We recall the equations that we want to pass to the limit into:

∂ t ρ ε + div(ρ ε v ε ) = ε∆ρ ε , (3.27)
∂ t (ρ ε c ε ) + div(ρ ε c ε v ε ) = -εc ε ∆ρ ε + div(b(c ε )∇µ ε ) + F cε , (3.28)
and for every η (sufficiently regular)

Ω ρ ε v ε (t) • η dx - Ω m 0 • η dx - t 0 Ω ρ ε v ε ⊗ v ε : ∇ x η dx ds - t 0 Ω p(ρ ε , c ε )div(η) dx ds + ε t 0 Ω (∇v ε ∇ρ ε ) • η dx ds + t 0 Ω T ε : ∇η dx ds + γ t 0 Ω ( 1 2 |∇c ε | 2 I -(∇c ε ⊗ ∇c ε )) : ∇η dx ds + Ω t 0 κ(c ε )v ε • η dx ds = 0. (3.29)
Using Proposition 3.13, which yields uniform estimates in ε, we pass to the limit in the previous equations. The difficult terms are the one involving nonlinear combinations. Indeed, it is not clear that we can obtain strong convergence of ρ ε as we have no estimates on higher order derivatives. We use the following lemma, see [START_REF]Compressible models[END_REF].

Lemma 3.15. Let g n , h n converge weakly to g, h respectively in L p 1 (0, T ; L p 2 (Ω)), L q 1 (0, T ; L q 2 (Ω)) where 1 ≤ p 1 , p 2 ≤ +∞ and

1 p 1 + 1 q 1 = 1 p 2 + 1 q 2 = 1.
We assume in addition that

∂g n ∂t is bounded in L 1 (0, T ; W -m,1 (Ω)) for some m ≥ 0 independent of n, (3.30)
and ∥h n -h n (t, • + ξ)∥ L q 1 (0,T ;L q 2 (Ω)) → 0 as |ξ| → 0, uniformly in n.

(3.31)

Then, g n h n converges to gh in the sense of distributions.

Remark 3.16. This lemma admits many variants, and it is possible to identify the weak limit of the products with lower regularity, we refer for instance to [START_REF] Moussa | Some variants of the classical Aubin-Lions lemma[END_REF].

We want to apply the previous lemma to the terms

ρ ε v ε , ρ ε c ε , ρ ε µ ε , ρ ε c 2 ε , ρ ε v ε , ρ ε v ε ⊗v ε , ρ ε v ε c ε .
We admit that ∂ρε ∂t , ∂ρεvε ∂t and ∂ρεcε ∂t satisfy (3.30) by using Proposition 3.13 and Equations (3.27)-(3.28)- (3.29). The compactness in space required in (3.31) also uses Proposition 3.13. We refer also to [START_REF] Cherfils | The compressible Navier-Stokes-Cahn-Hilliard equations with dynamic boundary conditions[END_REF]Subsection 3.1] for similar results. The terms εc ε ∆ρ ε and ε t 0 Ω (∇v ε ∇ρ) • η dx ds converge to 0 (the first one in the distributional sense) thanks to estimates (A7)-(A8). It remains to pass to the limit in (i.e identifying the weak limits)

p(ρ ε , c ε ), 1 2 |∇c ε | 2 , ∇c ε ⊗ ∇c ε , (3.32) b(c ε )∇µ ε , F cε (ρ ε , c ε ), ρ ε ∂ c Ψ 0 . (3.33)
The convergence of the last term is used to identify ρµ. To prove the previous convergences, we need to prove strong compactness in L 2 of c ε , ∇c ε and convergence a.e. of ρ ε to use Vitali's convergence theorem. But they follow from the arguments in [START_REF] Abels | On a diffuse interface model for a two-phase flow of compressible viscous fluids[END_REF] and [START_REF] Cherfils | The compressible Navier-Stokes-Cahn-Hilliard equations with dynamic boundary conditions[END_REF]Section 3.3 and 3.4].

Altogether, we can pass to the limit in every term of the equations. This concludes the argument.

Sending ε Q → 0

The last step in our proof is to let ε Q vanishes and recover the existence of weak solutions for the double well potential

Q(c) = 1 4 c 2 (1 -c) 2 .
Since we have the energy estimates from before, the work is essentially the same but we have to be careful about two points. The first one is to indeed have an energy estimate independent of ε Q . We discussed this point after Equation (3.25) and, hence, we do not repeat it here. The second point are the estimates obtained in Proposition 3.13. However, the estimates are essentially the same, except for estimate (A9) which becomes

{ρ∂ c ψ 0 } in L ∞ (0, T ; L 2a a+2 (Ω)). (3.34)
This can be proved knowing that, when ε Q ≈ 0, we have that for c large ρQ ′ ε Q (c) ≈ ρc 3 , and we use estimates (A2)-(A8). Altogether, the reasoning to pass to the limit is the same and we conclude.

Numerical scheme for G-NSCH model

We recall that in the following, we use the assumptions on the functionals provided in Section 2 for the numerical part. These particular choices lead to stability problems and degeneracy of mobility, and friction in certain regions. Indeed, our model comprises a Navier-Stokes part that needs to be stabilized to be simulated and a degenerate Cahn-Hilliard system that is well-known to be difficult to simulate numerically (see our introductory section). We here propose a numerical scheme that combines recent advances in numerical analysis. We use the numerical scheme for the simplified variant of the compressible NSCH system in [START_REF] He | Numerical study of compressible Navier-Stokes-Cahn-Hilliard system[END_REF] and fast structure-preserving scheme for degenerate parabolic equations [START_REF] Huang | Bound/positivity preserving and energy stable scalar auxiliary variable schemes for dissipative systems: applications to Keller-Segel and Poisson-Nernst-Planck equations[END_REF][START_REF] Huang | Bound/positivity preserving and unconditionally stable schemes for a class of fourth order nonlinear equations[END_REF]. Namely, we adapt the relaxation [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] of the Navier-stokes part as used in [START_REF] He | Numerical study of compressible Navier-Stokes-Cahn-Hilliard system[END_REF]. The Cahn-Hilliard part is stabilized using the SAV method [START_REF]The scalar auxiliary variable (SAV) approach for gradient flows[END_REF]. More precisely, a variant designed for degenerate parabolic models that preserves the physical bounds of the solution [START_REF] Huang | Bound/positivity preserving and energy stable scalar auxiliary variable schemes for dissipative systems: applications to Keller-Segel and Poisson-Nernst-Planck equations[END_REF][START_REF] Huang | Bound/positivity preserving and unconditionally stable schemes for a class of fourth order nonlinear equations[END_REF]. Indeed, we expect that the volume fraction c remains within the physically (or biologically) relevant bounds c ∈ (0, 1) since we are here using a double-well logarithmic potential. Thus, following [START_REF] Huang | Bound/positivity preserving and energy stable scalar auxiliary variable schemes for dissipative systems: applications to Keller-Segel and Poisson-Nernst-Planck equations[END_REF][START_REF] Huang | Bound/positivity preserving and unconditionally stable schemes for a class of fourth order nonlinear equations[END_REF], we construct the invertible mapping T : R → (0, 1),

with c = T (v) transforming Equations (1.2)-(1.3) into ρ (∂ t v + (v • ∇)v) = 1 T ′ (v) (div(b(c)∇µ) + F c ) , ρµ = -γT ′ (v)∆v -γT ′′ (v)|∇v| 2 + ρ ∂ψ 0 ∂c . (4.1)
Following [START_REF] Huang | Bound/positivity preserving and energy stable scalar auxiliary variable schemes for dissipative systems: applications to Keller-Segel and Poisson-Nernst-Planck equations[END_REF] and [START_REF] Huang | Bound/positivity preserving and unconditionally stable schemes for a class of fourth order nonlinear equations[END_REF], we can choose

T (v) = 1 2 tanh(v) + 1 2 , or T (v) = 1 1 + exp(-v) ,
thus preserving the bounds c ∈ (0, 1).

The SAV method allows to solve efficiently (and also linearly) the nonlinear Cahn-Hilliard part while preserving the dissipation of a modified energy. In the following, we assume that it exists a positive constant C such that the energy associated with the Cahn-Hilliard part, i.e.

E[t](ρ, c) = Ω γ 2 |∇c| 2 + ρψ 0 (ρ, c) = E 0 [t] + E 1 [t],
with E 1 the nonlinear part of the energy, and E 0 the linear part, is bounded from below, i.e.

E 1 + C ≥ 1 We define r(t) = E(t) + C 0 , with C 0 = 2C + ∥E(ρ 0 , c 0 )∥ L ∞ (Ω) ,
and apply the SAV method. System (4.1) becomes

ρ (∂ t v + (v • ∇)v) = 1 T ′ (v) (div(b(c)∇µ) + F c ) , ρµ = -γT ′ (v)∆v -γT ′′ (v)|∇v| 2 + ρ ∂ψ 0 ∂c , dr dt = - r(t) E[t] + C 0 Ω b(c)|∇µ| 2 -µF c dx, (4.2)
One can easily see that the previous modifications do not change our system at the continuous level.

Remark 4.1. As stated in [START_REF] Huang | Bound/positivity preserving and unconditionally stable schemes for a class of fourth order nonlinear equations[END_REF], the transformation T (v) = 1 1+exp(-v) allows us to write

log(c) -log(1 -c) = v, T ′ (v) = (1 -c)c, T ′′ (v) T ′ (v) = 1 -2c.

One-dimensional scheme

We consider our problem in a one-dimensional domain Ω = (0, L). Even though v is now a scalar, we still denote it in bold font. As mentioned previously, we relax the Navier-Stokes part of our system. Namely, we introduce a relaxation parameter η ≥ 0 and U = (ρ, ρv). We rewrite Equation (1.4) as (up to rescaling ν(c) as 3 4 ν(c))

∂ t U + ∂ x V = G(U ), ∂ t V + A∂ x U = -1 η (V -F (U )), (4.3)
in which G(U ) = (0, -κv),

F (U ) = (ρv, ρv 2 + p -ν(c)∂ x v + γ 2 |∂ x c| 2 ) and A = diag(a 1 , a 2 ) satisfying Liu's subcharacteristic condition A ≥ F ′ (U ), ∀U.
In what follows, and following [START_REF] He | Numerical study of compressible Navier-Stokes-Cahn-Hilliard system[END_REF], we use

a 1 = a 2 = max sup v + ∂ ρ p 2 , sup v -∂ ρ p 2 .
We discretize the domain using a set of N x nodes located at the center of control volumes of size ∆x such that Ω = j=0,...,Nx [x j-1 2 , x j+ 1 2 ]. The time interval [0, T ] is also discretized using a uniform time step ∆t. Our scheme follows the discrete set of equations

U * j = U n j , (4.4)
V * j = V n j - ∆t η V * j -F (U * j ) , (4.5)
U n+1 j = U * j - ∆t ∆x V * j+ 1 2 -V * j-1 2 + ∆tG(U n+1 j ), (4.6)
V n+1 j = V * j - ∆t ∆x A U * j+ 1 2 -U * j-1 2 , (4.7)
v n+1 j -v n j ∆t + v n+1 j • (∇v n+1 ) j = g(c n , µ n+1 , ρ n+1 ) j , (4.8)
g(c n , µ n+1 , ρ n+1 ) j = 1 T ′ (v n j )ρ n+1 j 1 ∆x (b(c n )∇µ n+1 ) j+ 1 2 -(b(c n )∇µ n+1 ) j-1 2 + F c (ρ n j , c n j ) , (4.9)
µ n+1 j = 1 ρ n+1 j -γT ′ (v n j )(∆v n+1 ) j -γT ′′ (v n j )|(∇v n ) j | 2 + ∂ψ 0 ∂c n j , (4.10
)

Ω λT (v n+1 ) dx = Ω c n + ∆tF c dx, (4.11)
c n+1 j = λ j T (v n+1 j ), (4.12
)

1 ∆t r n+1 -r n = - r n+1 E(c n+1 ) + C 0 ∆x j b(c n+1 j )|(∇µ n+1 ) j | 2 + + r n E(c n+1 ) + C 0 ∆x j µ n+1 j F c (ρ n+1 j , c n+1 j ), (4.13)
ξ n+1 = r n+1 E(c n+1 ) + C 0 , (4.14)
c n+1 j = ν n+1 c n+1 j , with ν n+1 = 1 -(1 -ξ n+1 ) 2 , (4.15)
v n+1 j = ν n+1 v n+1 j . ( 4
, U * j-1 2 and V * j+ 1 2 , V * j-1 2
, we use the an upwind method. We also mention that similarly to [START_REF] He | Numerical study of compressible Navier-Stokes-Cahn-Hilliard system[END_REF], one can implement a MUSCL scheme (see e.g. [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]) to obtain a higher order reconstruction. The upwind method permits to rewrite Equations (4.6)-(4.7) as

U n+1 j = U * j - ∆t 2∆x (V * j+1 -V * j-1 ) + ∆t 2∆x √ a(δ 2 x U * j ) + ∆tG(U n+1 j ), (4.17)
V n+1 j = V * j - a∆t 2∆x (U * j+1 -U * j-1 ) + ∆t 2∆x √ a(δ 2 x V * j ), (4.18)
where we used the notation δ 2 x u = u j+1 -2u j + u j-1 . In Equations (4.17)-(4.18), we emphasize that U * = U n and

V * = V n -∆t ∆x (V * -F (U n )).
In the following, we use the notations, ⟨U, V ⟩ = ∆x j U j V j , and ∥U ∥2 = ⟨U, U ⟩.

We also use ∆ 0,x U := 1 2 (U j+1 -U j-1 ). Our numerical scheme possesses the following important properties: Proposition 4.4 (Energy stability, bounds and mass preserving). Assuming the CFL-like condition ∆t ∆x √ a 1 ≤ 1 and the condition

∆t < E(c n+1 ) + C 0 C Ω |µ n+1 | dx , (4.19)
our numerical scheme satisfies the energy dissipation-like inequality

∥ √ aU n+1 ∥ + ∥V n+1 ∥ + r n+1 ≤ ∥ √ aU n ∥ + ∥V ⋆ ∥ + C n+1 r n , (4.20)
where

C n+1 = 1 + ∆t E(c n+1 )+C 0 Ω µ n+1 F c (ρ n+1 , c n+1 ) dx 1 + ∆t E(c n+1 )+C 0 Ω b(c n+1 )|∇µ n+1 | 2 dx .
Furthermore, the numerical scheme preserves the physically relevant bounds of the mass fraction, i.e. 0 < c n+1 < 1.

Finally, the scheme is mass preserving, i.e.

x j ρ n+1 j =

x j ρ n j .

Remark 4.5. Note that the constant C n+1 is smaller than 1 whenever the nonnegative part of the dissipation of the energy is greater than the increase of energy induced by the source term F c . This of course satisfied when we have F c = 0 for instance. When the source term is dominant, we only know that C n+1 behaves like 1 + C∆t.

Proof. We start with Equation (4.17), and using the definition of the function G(U n+1 j ) as well as assuming κ(c) ≥ 0 (for c ∈ R), after taking the square on both sides, multiplying by ∆x and summing over the nodes j = 0, ..., N x , we have

∥U n+1 ∥ 2 ≤ ∥U n ∥ 2 + ∆t 2∆x 2 ∥∆ 0,x V ⋆ ∥ 2 + ∆t √ a 2∆x 2 ∥δ 2 x U n ∥ 2 - ∆t ∆x ⟨∆ 0,x V * , U n ⟩ + ∆t √ a ∆x ⟨U n , δ 2 x U n ⟩ - √ a∆t 2 2∆x 2 ⟨∆ 0,x V ⋆ , δ 2 x U n ⟩.
Repeating the same computations for Equation (4.18), we have

∥U n+1 ∥ 2 ≤ ∥U n ∥ 2 + a∆t 2∆x 2 ∥∆ 0,x U n ∥ 2 + ∆t √ a 2∆x 2 ∥δ 2 x V ⋆ ∥ 2 - a∆t ∆x ⟨∆ 0,x U n , V ⋆ ⟩ + ∆t √ a ∆x ⟨U ⋆ , δ 2 x V ⋆ ⟩ - a
At this point, the proof is similar to the proof of [START_REF] He | Numerical study of compressible Navier-Stokes-Cahn-Hilliard system[END_REF]Theorem 4.1] (these steps use the periodic boundary conditions and the summation by parts formula to cancel some terms when summing both of the previous equations together), to obtain for a constant C > 0,

∥ √ aU n+1 ∥ 2 + ∥V n+1 ∥ 2 ≤ C ∥ √ aU n ∥ 2 + ∥V n ∥ 2 .
Then, for the Cahn-Hilliard part, we easily obtain from Equation (4.13)

r n+1 = r n 1 + ∆t E(c n+1 )+C 0 Ω µ n+1 F c (ρ n+1 , c n+1 ) dx 1 + ∆t E(c n+1 )+C 0 Ω b(c n )|∇µ n+1 | 2 .

So as long as

1 + ∆t E(c n+1 ) + C 0 Ω µ n+1 F c (ρ n+1 , c n+1 ) dx ≥ 0, so does r n+1 . This implies that, assuming ∥F c ∥ L ∞ < C, ∆t < E(c n+1 ) + C 0 C Ω |µ n+1 | dx .
Under the previous condition, if r n ≥ 0 so does r n+1 , and (4.20) follows. Then, from the definition of ξ n+1 and C 0 , we have

0 < ξ n+1 < r 0 E(c n+1 ) + C 0 ≤ 2,
from the definition of the constant C 0 .

Remark 4.6. We emphasize that analytically, we can not verify the condition (4.19) as the solution can tend to 0 or 1 and consequently the integral of |µ| can blow up. However, we observe during numerical simulations that the condition (4.19) is obtained for reasonably small ∆t. We also note that if we do not consider any source term, i.e. F c = 0, the scheme satisfies the dissipation relation

∥ √ aU n+1 ∥ + ∥V n+1 ∥ + r n+1 ≤ ∥ √ aU n ∥ + ∥V ⋆ ∥ + r n ,
with the stability condition ∆t ∆x √ a ≤ 1.

Two-dimensional scheme

We describe the two-dimensional scheme. This scheme possesses the same properties as the onedimensional scheme. The proof of this is easily obtained from a very simple adaptation of the proof of Proposition 4.4. We write the velocity field v = (u x , u y ). System (1.1)-(1.4) with the transformation proposed at the beginning of this section, reads

∂ t ρ + ∂ x (ρu x ) + ∂ y (ρu y ) = 0, (4.21)
∂ t ρ u x u y + ∂ x (ρu 2 x + p) ∂ y (ρu 2 y + p) + ∂ y (ρu x u y ) ∂ x (ρu x u y ) = 2 ∂ x (ν(c)∂ x u x ) ∂ y (ν(c)∂ y u y ) + ∂ y (ν(c)(∂ y u x + ∂ x u y )) ∂ x (ν(c)(∂ y u x + ∂ x u y )) - 2 3
∂ x (ν(c)(∂ x u x + ∂ y u y )) ∂ y (ν(c)(∂ x u x + ∂ y u y )) - γ 2
∂ x ((∂ x c) 2 -(∂ y c) 2 ) ∂ y ((∂ y c) 2 -(∂ x c) 2 ) -γ ∂ y (∂ x c∂ y c) ∂ x (∂ x c∂ y c) -κ(ρ, c) u x u y , (4.22) ρ (∂ t v + u x ∂ x v + u y ∂ y v) = 1 T ′ (v) (∂ x (b(c)∂ x µ) + ∂ y (b(c)∂ y µ)) + 1 T ′ (v) F c , (4.23) ρµ = -γT ′ (v)(∂ xx c + ∂ yy c) -γT ′′ (v) (∂ x v) 2 + (∂ y v) 2 + ρ ∂ψ 0 ∂c , (4.24
)

dr dt = - r(t) E[t] + C 0 Ω b(c)|∇µ| 2 -µF c dx. (4.25)
We introduce the notations U = (ρ, ρu x , ρu y ), G(U ) = (0, -κu x , -κu y ) and

F (U ) = (ρu x , ρu 2 x + p -2ν(c)∂ x u x + 2 3 ν(c)(∂ x u x + ∂ y u y ) + 1 2 γ (∂ x c) 2 -(∂ y c) 2 , ρu x u y -ν(c) (∂ y u x + ∂ x u y ) + γ∂ x c∂ y c), K(U ) = (ρu y , ρu x u y -ν(c) (∂ y u x + ∂ x u y ) + γ∂ x c∂ y c, ρu 2 y + p -2ν(c)∂ y u y + 2 3 ν(c)(∂ x u x + ∂ y u y ) + 1 2 γ (∂ x c) 2 -(∂ y c) 2 ).
The stabilization (see [START_REF] He | Numerical study of compressible Navier-Stokes-Cahn-Hilliard system[END_REF][START_REF]The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF]) of the Navier-Stokes part of our system reads

     ∂ t U + ∂ x V + ∂ y W = G(U ), ∂ t V + A∂ x U = -1 η (V -F (U )), ∂ t W + B∂ y U = -1 η (W -K(U )), (4.26) in which A = diag(a 1 , a 2 , a 3 ) and B = diag(b 1 , b 2 , b 3 ).
In the following, we choose

a 1 = a 2 = a 3 = max{sup u x + ∂ ρ p 2 , sup u 2 x , sup u x -∂ ρ p 2 }, b 1 = b 2 = b 3 = max{sup u y + ∂ ρ p 2 , sup u 2 y , sup u y -∂ ρ p 2 }.
We assume that our two-dimensional domain is a square [0, L] × [0, L]. We discretize the domain using square control volumes of size ∆x × ∆y. The cell centers are located at positions (x j , y j ), and we approximate the value of a variable at the cell center by its mean, e.g.

ρ j,i = 1 ∆x∆y x j+ 1 2 x j-1 2 y j+ 1 2 y j-1 2 ρ(x, t) dx.
Simply employing a first-order time discretization, the numerical scheme becomes

U * j,i = U n j,i , (4.27) V * j,i = V n j,i - ∆t η V * j,i -F (U * j,i ) , (4.28) W * j,i = W n j,i - ∆t η W * j,i -K(U * j,i ) , (4.29) U n+1 j,i = U * j,i - ∆t ∆x V * j+ 1 2 ,i -V * j-1 2 ,i - ∆t ∆y W * j,i+ 1 2 -W * j,i-1 2 + ∆tG(U n+1 i,j ), (4.30) V n+1 j,i = V * j,i - ∆t ∆x A U * j+ 1 2 ,i -U * j-1 2 ,i , (4.31)
W n+1 j,i = W * j,i - ∆t ∆y B U * j,i+ 1 2 -U * j,i-1 2 , (4.32) v n+1 j,i -v n j,i ∆t + v n+1 j,i • (∇v n+1 ) j,i = g(c n , µ n+1 , ρ n+1 ) j,i , (4.33
)

g(c n , µ n+1 , ρ n+1 ) j,i = 1 T ′ (v n j,i )ρ n+1 j,i ∆x (b(c n )∇µ n+1 ) j+ 1 2 ,i -(b(c n )∇µ n+1 ) j-1 2 ,i + 1 T ′ (v n j,i )ρ n+1 j,i ∆y (b(c n )∇µ n+1 ) j,i+ 1 2 -(b(c n )∇µ n+1 ) j,i-1 2 + F c (ρ n j,i , c n j,i ) T ′ (v n j,i )ρ n+1 j,i , (4.34)
µ n+1 j,i = 1 ρ n+1 j,i -γT ′ (v n j,i )(∆v n+1 ) j,i -γT ′′ (v n j,i )|(∇v n ) j,i | 2 + ∂ψ 0 ∂c n j,i , (4.35
)

Ω λT (v n+1 ) dx = Ω c n + ∆tF c dx, (4.36)
c n+1 j,i = λT (v n+1 j,i ), (4.37)

1 ∆t r n+1 -r n = - r n+1 E(c n+1 ) + C 0 Ω b(c n+1 )|∇µ n+1 | 2 dx+ + r n E(c n+1 ) + C 0 Ω µ n+1 F c (ρ n+1 , c n+1 ) dx, (4.38)
ξ n+1 = r n+1 E(c n+1 ) + C 0 , (4.39) c n+1 j,i = ν n+1 c n+1 j,i , with ν n+1 = 1 -(1 -ξ n+1 ) 2 , (4.40) v n+1 j,i = ν n+1 v n+1 j,i . (4.41)

Numerical experiments

In the section, we begin using the one-dimensional scheme (4.4)-(4.16) with no source term and friction, i.e. κ(ρ, c) = 0 and F c = 0, and we verify that the scheme preserves all the properties stated in Proposition 4.4. Furthermore, we verify the spatial and temporal convergence orders for the scheme.

Remark 5.1 (Implementation details). All numerical schemes are implemented using Python 3 and the Numpy and Scipy modules. The linear systems present in the schemes are solved using the Generalized Minimal RESidual iteration (GMRES) indirect solver (function available in the scipy.sparse.linalg module).

One dimensional numerical test case 1: non-matching densities

We start with a one-dimensional test case to show the spatiotemporal evolution of the density, mass fraction, and velocity. We also verify numerically the properties stated in Proposition 4.4.

In this subsection, we use the double-well logarithmic potential

ψ mix = 1 2 (α 1 (1 -c) log(ρ(1 -c)) + α 2 c log(ρc)) - θ 2 (c - 1 2 ) 2 + k,
with k = 100, α 1 = 1.2, α 2 = 0.5, and θ = 4. This allows us to model a fluid for which the phase denoted by the index 1 is denser compared to the phase indicated by the index 2. Indeed, this can been seen on the effect of the values α 1 and α 2 on the potential. Taking α 2 < α 1 shifts the well corresponding to phase 1 very close to 1 compared to the other phase. This models the fact that the fluid 1 is in fact more compressible and thus aggregates of pure phase 1 appear denser. We use the computational domain Ω = (0, 1) discretized in N x = 128 cells. We take T = 0.3 (this has been chosen because the system reaches a meta-stable state by that time) and use the initial time step ∆t = 1 × 10 -6 (this time step size is adapted from the CFL-like condition stated in Proposition 4.4, however, the time step will never be larger than ∆t max = 1 × 10 -5 ). We choose the width of the diffuse interface to be γ = 1/500, the viscosity to be constant ν(c) = 1 × 10 -2 , the relaxation parameter to be η = 1 × 10 -3 , and the exponent for the barotropic pressure equals to a = 3. We choose constant initial conditions for the density and the pressure, i.e.

ρ 0 = 0.9, v 0 = 1.0.

The initial mass fraction is assumed to be a constant with a small random noise, i.e.

c 0 = c -0.05r, with c = 0.5 and r is a vector of random values between 0 and 1 given by the uniform distribution.

Figure 1 shows snapshots of (ρ, c, v) at different times (starting from the initial condition (i.e. Figure 1a) to the stable state depicted in Figure 1d). We observe that the numerical scheme catches well the spinodal decomposition of the binary fluid while it is transported to the right. Indeed, after an initial regularization of the initial condition (Fig. 1b), the separation of the two phases of the fluid occurs and small aggregates appear (Fig. 1b). Then, the coarsening of the small aggregates into larger ones occurs. We arrive at the end of the simulation to the solution depicted in Fig. 1d. During the last time steps, the evolution of the mass fraction was very slow. This slow evolution is due to the degeneracy of the potential. The solution in the longtime is expected to depict evenly spaced aggregates of equal size. We emphasize that even at the end of the simulation, the fluid continues to move to the right since the velocity v is positive. Due to the fact that α 1 > α 2 , fluid 1 is denser. This is observed in Fig. 1d as the total fluid density ρ and the pressure p are larger in the zones c ≈ 1 (corresponding to the pure phase of fluid 1). We also observe that the pressure in pure aggregates of fluid 1 is larger compared to the other zones. Figures 2 show that our numerical scheme preserved the properties presented in Proposition 4.4. Indeed, denoting by

dE dt = ∥ √ aU n+1 ∥ + ∥V n+1 ∥ + r n+1 -∥ √ aU n ∥ + ∥V n ∥ + C n+1 r n ,
we observe using Fig. 2a that Inequality (4.

Convergence tests

We study the numerical convergence of our one dimensional scheme (4.4)-(4.16) with κ(ρ, c) = 0 and F c = 0. The computational domain is Ω = (0, 1). The final time is T = 0.1 and choose α 1 = α 2 = 1. The other parameters γ, β, η, ν are chosen as for test case 1 (see Subsection 5.1). The initial condition for the mass fraction is given by c 0 = 0.5 + 0.01 cos(6πx).

The initial conditions for the velocity and the total density are chosen as in test case 1.

Convergence in space

We start by fixing the time step ∆t = 1 × 10 -5 and we vary the grid size. We vary the number of points N x = {64, 128, 256, 512, 1024, 2048}. For each ∆x, we compute the error at time t = T = 0.1 (we chose this value since we observed that the solution reached a stable state at that time) following

error = ∥ρ ∆x -ρ ∆x/2 ∥ + ∥c ∆x -c ∆x/2 ∥ + ∥v ∆x -v ∆x/2 ∥, (5.1)
in which each of the norms is computed following

∥ρ ∆x -ρ ∆x/2 ∥ =   ∆x Nx j=1 ρ ∆x (x 2(j-1)+1) ) -ρ ∆x/2 (x j ) 2 + ρ ∆x ((x 2(j-1)+2) ) -ρ ∆x/2 (x j )) 2 )   1/2
, with N x the number of points on the ∆x/2 grid, hence, the solution on the coarse grid (i.e. ∆x/2) is extended on the fine grid ∆x. We arrive at the results given in Figure 3. As expected by the upwind scheme, the spatial order of convergence is a little less than 1 for the total density ρ (see Figure 3a) and the velocity v (see Figure 3c).

Convergence in time

We here fix the grid size and select N x = 128 points. We choose ∆t = 1 × 10 -4 , and vary the time steps according to ∆t array = {∆t, ∆t 2 , ∆t 4 , ∆t 8 , ∆t 16 , ∆t 32 }. The other parameters and initial conditions are chosen as in the spatial convergence test (see Subsection 5.2.1). The error between the reference solution and a computed solution at time t = T = 0.2 is given by

∥ρ ∆t -ρ ∆t/2 ∥ =   ∆x Nx j=1 (ρ ∆t (x j ) -ρ ∆t/2 (x j )) 2   1/2 .
We obtain the results depicted in Figure 4 We observe that the order of convergence in time for our scheme is exactly 1 as expected.

Conclusion and simulation of a growing tumor in a healthy tissue

To conclude this article, we would like to go back to our initial interest, i.e. the modeling of tumor growth. We now present a simple numerical simulation for this application and introduce a forthcoming work focusing only on the modelling of tumor. Here, we assume that F c acts as a transfer of mass from the healthy tissue to the cancerous population (this can be viewed as tumor cells using material from the environment such as nutrients to divide and grow). Furthermore, we assume that the two cell populations adhere with different strengths to the ECM. Altogether, we choose

F c = 20ρc(1 -c/c * ), c * = 0.9, κ(ρ, c) = (κ 1 ρc + κ 2 ρ(1 -c)) , γ = 1 1000
, with c * representing the maximum saturation (decided to be c * = 1 in the following) and κ 1 = 0 and κ 2 = 20.

For this test case, we choose the initial conditions as follows:

c 0 = 0.008 + 0.6 exp -100(x -0.5) 2 -100(y -0.5) 2 , ρ 0 = c + 0.5 (1 -c),

v 0 = 0,
We chose a square domain Ω of length 1 and 64 × 64 cells for the spatial discretization. The other parameters that we did not specify are chosen as in Subsection 5.1. Figs 5 show the spatio-temporal evolution of the tumor with α 1 = α 2 = 1. We observe that as the tumor grows the shape remains radially symmetrical. We think that, due friction effects, instabilities may appear. The reason we are not able to see them in that case is double. Firstly, the scheme is less than first order in space and small structures may not be captured. Secondly, the regularization given by the double Laplacian from the Cahn-Hilliard equation may be too strong and the finger-like instabilities are not seen.

To verify the latter assertion, we propose another test case with the parameters

α 1 = 1.2, α 2 = 0.8, κ 2 = 100, γ = 1 1500
, and c 0 = 0.01 + 0.9 exp -100(x -0.5) 2 -100(y -0.5) 2 .

The other parameters and the initial conditions are chosen as the previous test case. Figs 6 show our numerical results for this "tumor growth" test case. We observe that as the tumor grows the shape does not remain perfectly symmetrical. This result emphasize that the behavior of the solution depends strongly on the regularization provided by the bi-laplacian term and the strength of the friction on the rigid fibers of the medium. Furthermore, one possible other explanation of the instabilities could be the difference of densities of between the two fluids. This is controlled by taking α 1 ̸ = α 2 (see Figs 1 to better see one-dimensional numerical results highlighting the effect of non-matching densities). Indeed, from Figs 1 we know that taking α 2 > α 2 leads to a heterogeneous pressure field which leads to pressure gradients. Hence, since the initial velocity field is zero in our case, the direction of the velocity is given by -∇p. Consequently, v tends to transport the cell densities to regions of less pressure, i.e. away from the tumor cell clusters. This movement in addition to heterogeneous friction effects is known to produce complex patterns (see e.g. [START_REF] Lorenzi | On interfaces between cell populations with different mobilities[END_REF]). We argue that, from these numerical simulations, our model seems capable to qualitatively represent patterns of invasive growth of tumors and could unravel the possible mechanical effects playing a role in the emergence of heterogeneous structures observed in tumor invasion of healthy tissue. However, we emphasize that to achieve these latter goals, we have to be able to capture accurately the possible fine structures emerging during the numerical simulations. Thus, our numerical scheme will be improved to increase the spatial and temporal orders of accuracy by taking advantage of the flexibility of the SAV method. In a forthcoming work, we will develop a highorder finite element scheme for the compressible SAV-NSCH system we proposed in the present work. This numerical scheme will allow efficient simulations of compressible diphasic fluids and will be used to simulate relevant test cases with applications in fluid mechanics such as rising bubbles and Rayleigh-Taylor instabilities.

In previously mentioned future work, we will study the possible mechanical effects at play during invasive tumor growth. Our strategy will be relatively simple: As presented in Appendix A, we will reduce the G-NSCH system to better represent tumor growth while removing non-necessary effects such as inertia. Our goal will be to present numerical simulations capturing Saffman-Taylor-like instabilities depicted by the protrusions of the tumor in the healthy tissue and commonly observed in the context of, e.g. , skin cancer [START_REF] Chatelain | Emergence of microstructural patterns in skin cancer: a phase separation analysis in a binary mixture[END_REF].

On the basis of the present work, numerous other work directions can be envisioned. Indeed, the analytical aspects of this work can be improved. This direction is challenging because as pointed out in the present work, necessary assumptions on the potential and the mobility term to perform the solutions' existence proof do not allow to take physically or biologically relevant forms. In fact, singular potentials, degenerate mobilities and viscosity functions are not allowed.

Working in this direction would require us to use recent results concerning the Cahn-Hilliard equation as well as the compressible Navier-Stokes model and adapting them to the context of the NSCH system. One possible solution that we will investigate is to derive a Bresch-Desjardins entropy estimate [START_REF] Bresch | Existence of global weak solutions for a 2D viscous shallow water equations and convergence to the quasi-geostrophic model[END_REF][START_REF]On the construction of approximate solutions for the 2D viscous shallow water model and for compressible Navier-Stokes models[END_REF] for the compressible NSCH as it has been done recently by Vasseur and Yu [START_REF] Vasseur | Existence of global weak solutions for 3D degenerate compressible Navier-Stokes equations[END_REF], and Bresch, Vasseur and Yu [START_REF] Bresch | Global existence of entropy-weak solutions to the compressible Navier-Stokes equations with non-linear density dependent viscosities[END_REF], for the compressible Navier-Stokes model with degenerate viscosities. We would like to conclude this research article by presenting a possible research perspective the authors will investigate. We believe that the coupling proposed by the G-NSCH model between fluid movement and tumor growth could have an original application to the mathematical representation of tumor-on-chips [START_REF] Liu | Tumor-on-a-chip: from bioinspired design to biomedical application[END_REF]. Indeed, these microfluidic devices present a growing interest in oncology since they can replicate in a very accurate manner the micro-environment of the tumor. Therefore, they present a lot of advantages for experiments and the G-NSCH model could be used to replicate such experiments in an in-silico manner. We believe this could be useful to optimize the devices and also help to answer questions in oncology.

A.2 Mass balance equations

We assume that each component has its own velocity and the component 1 is proliferating. The fluid being saturated, i.e. c 1 + c 2 = 1 Therefore, we have the mass balance equations for each component

∂ϕ 1 ∂t + div (ϕ 1 v 1 ) = F 1 (ρ, c 1 , c 2 ), ∂ϕ 2 ∂t + div (ϕ 2 v 2 ) = F 2 (ρ, c 1 , c 2 ).
(A.4)

In the previous system, the functions F i (ρ, c 1 , c 2 ) (i = 1, 2) act as source terms of mass. Summing the two equations, we obtain the continuity equation for the total density of the mixture, and using the mass fractions (denoting c 1 = c) and the relations (A.1), we obtain the balance equation for the density of the mixture

∂ρ ∂t + div (ρv) = F 1 + F 2 =: F ρ . (A.5)
To obtain a system analogous to (A.4), we rewrite the first equation of (A.4) using the definition of the mass fraction (A.1) to obtain

∂ρc ∂t + div (ρcv 1 ) = F 1 (ρ, c, 1 -c) =: F c . (A.6)
The mass of the component 1 is transported by the average velocity v and the remaining diffusive flux J 1 = ρc (v -v 1 ). Therefore, we can replace the previous equation by

∂ρc ∂t + div (ρcv) = div (J 1 ) + F c .
Then, using the definition of the material derivative (A.3) and the mass balance equation for the total mixture (A.5), the left-hand side of the previous equation reads

∂ρc ∂t + div (ρcv) = ρ Dc Dt + c ∂t + div (ρv) = ρ Dc Dt + cF ρ .
Altogether, we obtain the balance equation for the mass fraction of the component 1

ρ Dc Dt = div (J 1 ) + F c -cF ρ . (A.7)
Since c 2 = 1 -c, solving the equations (A.5) and (A.7) is equivalent to solving the system (A.4).

In the following, we refer to c as the order parameter (terminology often used in the framework of the Cahn-Hilliard model [START_REF] Cahn | On spinodal decomposition[END_REF][START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF]).

A.3 Balance of linear momentum

We write the balance of linear momentum [START_REF] Eck | Mathematical modeling[END_REF], which describes the evolution of the velocity v due to internal stresses. Indeed, we neglect the effect of any external forces, including gravity. Following continuum mechanics, the Cauchy stress tensor gives the stresses acting inside the mixture due to viscous and non-viscous effects. An additional stress must be taken into account to represent the effect of concentration gradients [START_REF] Ericksen | Liquid crystals with variable degree of orientation[END_REF]. Altogether, we assume that the stress tensor is a function of the total density ρ, the order parameter c (i.e. the mass fraction of population 1), its gradient ∇c, and the total velocity of the mixture v i.e.

σ σ σ = σ σ σ(ρ, c, ∇c, v).

The friction around the pores of the medium is modeled by a drag term in the balance equation [START_REF] Narsilio | Upscaling of navier-stokes equations in porous media: Theoretical, numerical and experimental approach[END_REF] with a permeability coefficient κ(ρ, c) = κ 1 (ρ, c) + κ 2 (ρ, c) (the sum of the two friction coefficients for each component of the mixture). The permeability coefficient relates the properties of the fluid and the porous medium.

For each dimension (for example if d = 3, then j = {x, y, z}), the balance of linear momentum reads [START_REF] Eck | Mathematical modeling[END_REF] ∂ρv

j ∂t + div (ρv j v) = div (σ σ σ) j -κ(ρ, c)v j + F v j ,
where F v j (v j , ρ) represents the gain or loss of velocity in the j-th direction from different effects such as external forces. Then, using the continuity equation (A.5), we can rearrange the left-hand side to obtain

∂ρv j ∂t + div (ρv j v) = ρ Dv j Dt + v j ∂ρ ∂t + div (ρv) = ρ Dv j Dt + v j F ρ + F v j .
Therefore, we have

ρ Dv j Dt = div(σ σ σ) j -(κ(ρ, c) + F ρ ) v j + F v j .
Then, we can rewrite the balance of linear momentum in a more compact form

ρ Dv Dt = div(σ σ σ) -(κ(ρ, c) + F ρ ) v + F v , (A.8)
where F v (v, ρ) is the vector of coordinates F v j .

A.4 Energy balance

The total energy of the mixture is the sum of the kinetic energy ρ 1 2 |v| 2 and of the internal energy ρu, where u = u(ρ, c, ∇c) is a specific internal energy. Compared to the classical conservation law for the total energy, we have an additional energy flux τ τ τ Dc Dt . Indeed, due to the interface region, surface effects must be taken into account. Following this direction, Gurtin [START_REF] Gurtin | On a nonequilibrium thermodynamics of capillarity and phase[END_REF] proposed to include in the second law of thermodynamics, the effect of an additional force called the microscopic-stress which is related to forces acting at the microscopic scale. We denote this supplementary stress by τ τ τ . Since we assume that the system is maintained in an isothermal state, the balance equation for the energy is given by [START_REF] Eck | Mathematical modeling[END_REF]

∂ ∂t ρ 1 2 |v| 2 + ρu + div ρ 1 2 |v| 2 + u v = div σ σ σ T v + div τ τ τ Dc Dt -div (q) + ρg + c ρ F ρ + c c F c + c v F v , (A.9)
where q is the heat flux and ρg is the density of heat sources to maintain the temperature constant. The last three terms in Equation (A.9) account for the energy supply coming from the mass and velocity sources (see e.g. [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF][START_REF] Lam | Thermodynamically consistent Navier-Stokes-Cahn-Hilliard models with mass transfer and chemotaxis[END_REF]). The prefactors c ρ , c c , c v will be determined later to satisfy the free energy imbalance. Then, repeating the same calculations on the left-hand side to use the balance of mass (A.5), we get

∂ ∂t ρ 1 2 |v| 2 + ρu + div ρ 1 2 |v| 2 + u v = ρ D Dt 1 2 |v| 2 + u + 1 2 |v| 2 + u F ρ .
Applying the chain rule to the kinetic part, we get

ρ D Dt 1 2 |v| 2 = ρv • Dv Dt ,
and using the balance of linear momentum (A.8), we obtain

ρv • Dv Dt = v • div(σ σ σ) -(κ(ρ, c) + F ρ ) |v| 2 + F v • v.
Using these previous equations inside (A.9), we obtain the balance equation for the internal energy

ρ Du Dt = div σ σ σ T v -v • div (σ σ σ) + div τ τ τ Dc Dt + (κ(ρ, c) + F ρ ) |v| 2 -F v v -div (q) + ρg - 1 2 |v| 2 + u F ρ + c ρ F ρ + c c F c + c v F v . However, since v • (div (σ σ σ)) -div σ σ σ T v = -σ σ σ . . ∇v,
where ∇v = ∂ x j v i i,j=1,...,d is the Jacobi matrix and, we have A . . B = i,j A ij B ij , for two matrices A, B. Altogether, we have the balance equation for the internal energy

ρ Du Dt = σ σ σ . . ∇v + div τ τ τ Dc Dt + (κ(ρ, c) + F ρ ) |v| 2 -F v v -div (q) + ρg - 1 2 |v| 2 + u F ρ + c ρ F ρ + c c F c + c v F v .
(A.10)

A.5 Entropy balance and Clausius-Duhem inequality

We aim to apply the second law of thermodynamics. To do so, we define the entropy s = s(ρ, c, ∇c) and the Helmholtz free energy F = F(ρ, c, ∇c), both related through the equation Then, using the definition of the Helmholtz free energy (A.11) and the balance of energy (A.10), we obtain

F = u -T s, (
ρ Ds Dt = - ρ T DF Dt + ρ T Du Dt = - ρ T DF Dt + 1 T σ σ σ . . ∇v + div τ τ τ Dc Dt + (κ(ρ, c) + F ρ ) |v| 2 -F v v -div (q) + ρg - 1 2 |v| 2 + u F ρ + c ρ F ρ + c c F c + c v F v , (A.13)
where we have replaced the material derivative of the internal energy using its balance equation (A.10).

The constitutive relations for the functions constituting the Navier-Stokes-Cahn-Hilliard model are often derived to satisfy the Clausius-Duhem inequality (Coleman-Noll Procedure) [START_REF] Eck | Mathematical modeling[END_REF]. Indeed, this inequality provides a set of restrictions for the dissipative mechanisms occurring in the system. However, in our case, due to the presence of source terms, we can not ensure that this inequality holds without some assumptions on the proliferation and friction of the fluid around the pores. Therefore, we use here a different method: the Lagrange multipliers method. Indeed, the Liu [START_REF] Liu | Method of Lagrange multipliers for exploitation of the entropy principle[END_REF] and Müller [START_REF] Müller | Thermodynamics[END_REF] method is based on using Lagrange multipliers to derive a set of restrictions on the constitutive relations that can be applied even in the presence of source terms. Following classical Thermodynamics [START_REF] Müller | Thermodynamics[END_REF], we state the second law as an entropy inequality, i.e., the Clausius-Duhem inequality in the local form [START_REF] Eck | Mathematical modeling[END_REF] ρ Ds Dt ≥ -div q T + ρg T + div (J ) , (A. [START_REF] Bresch | Global existence of entropy-weak solutions to the compressible Navier-Stokes equations with non-linear density dependent viscosities[END_REF] where J is the entropy flux. The inequality (A.14) results from the fact that the entropy of the mixture can only increase. Using the equation (A.13), we obtain

ρ T DF Dt - 1 T σ σ σ . . ∇v + div τ τ τ Dc Dt + (κ(ρ, c) + F ρ ) |v| 2 -F v v - 1 2 |v| 2 + u F ρ + c ρ F ρ + c c F c + c v F v + div (J ) ≤ 0.
ρ (-ρdiv(v) + F ρ ) ∂F ∂ρ + Dc Dt ∂F ∂c + ∇ Dc Dt -(∇v) T ∇c • ∂F ∂∇c -div τ τ τ Dc Dt -σ : ∇v -(κ(ρ, c) + F ρ ) |v| 2 -F v v - 1 2 |v| 2 + u F ρ + c ρ F ρ + c c F c + c v F v + T div (J ) ≤ 0.
(A.16) By the chain rule, we have

div τ τ τ Dc Dt = τ τ τ ∇ Dc Dt + Dc Dt div (τ τ τ ) .
Furthermore, we know that

-ρ 2 div (v) ∂F ∂ρ = -ρ 2 ∂F ∂ρ 1 . . ∇v, and
-ρ (∇v) T ∇c • ∂F ∂∇c = -ρ ∇c ⊗ ∂F ∂∇c . . ∇v.
Gathering the previous three relations and reorganizing the terms of (A.16), we obtain

-ρ 2 ∂F ∂ρ 1 -ρ∇c ⊗ ∂F ∂∇c -σ σ σ . . ∇v + ρ ∂F ∂c -div(τ τ τ ) Dc Dt + ρ ∂F ∂∇c -τ τ τ ∇ Dc Dt + T div (J ) -(κ(ρ, c) + F ρ ) |v| 2 -F v v - 1 2 |v| 2 + u -ρ ∂F ∂ρ F ρ + c ρ F ρ + c c F c + c v F v ≤ 0.
(A.17)

Then, we use Liu's Lagrange multipliers method [START_REF] Liu | Method of Lagrange multipliers for exploitation of the entropy principle[END_REF]. We denote by L c the Lagrange multiplier associated with the mass fraction equation (A.7). The method of Lagrange multipliers consists in setting the following local dissipation inequality that has to hold for arbitrary values of (ρ, c, ∇ρ, ∇c, v, p) Furthermore, using the same simplifying assumption as in Abels and Feireisl [START_REF] Abels | On a diffuse interface model for a two-phase flow of compressible viscous fluids[END_REF] to avoid vacuum zones, our final reduced system of equations is Problem 2: Biologically relevant variant of the system. For this variant of the system, we assume the production of mass and neglect certain effects. Namely, we neglect inertia effects, and the viscosity of the fluid, and assume no external source of velocity. This leads to the momentum equation

∇p + κ(ρ, c)v = -γdiv ∇c ⊗ ∇c - 1 2 |∇c| 2 1 -F ρ v.
Then, we assume one cell population proliferates while the other does not, leading to where C b is a positive constant.

Remark B.2. Note that we cannot prove the existence of solutions for this mobility because it is degenerate. Instead, in the analysis part of this paper, we use an approximate mobility bounded away from zero. The difficulty of the degenerate mobility comes when one tries to identify the chemical potential µ when the parameters in the approximating schemes are removed. Indeed the estimate of the dissipation of the energy Ω b(c)|∇µ| 2 dx does not provide estimates on ∇µ anymore.

We use for the pressure a power law such that p e (ρ) = 1 a -1 ρ a-1 . (B.13)

For H(c) and G(c), two choices can be considered depending on the behavior of the cells we want to represent. If the two cell populations exert attractive forces when they recognize cells of the same type and repulsion with the other type, the potential has to take a form of a double-well for which the two stable phases are located at the bottom of the two wells (see e.g. Figure 7a). This is a situation close to the phase separation in binary fluids. Thermodynamically consistent potentials are of Ginzburg-Landau type with the presence of logarithmic terms. Even though the double-well form of the potential is originally used for applications dealing with materials, it can also be motivated for biological purposes. Indeed, considering an application where the mixture is saturated with two cell types, a double-well potential is biologically relevant and reflects correctly the expected behavior of cells: they are attracted to each other respectively to their cell type at low densities and after a certain density they start to repel each other to avoid the creation of overcrowded zones. A typical example of biologically relevant double-well potential is given by thus giving

ψ mix = 1
H(c) = 1 2 (α 1 (1 -c) + α 2 c) , Q(c) = 1 2 (α 1 (1 -c) log(1 -c) + α 2 c log(c)) - θ 2 (c - 1 2 ) 2 + k,
where θ > 1 and k is an arbitrary constant.

To meet the phenomenological observations of the interaction between cells when the mixture is composed of only one cell population (and the other component of the mixture is supposed to be much more compressible), a single-well potential seems more appropriate [START_REF] Byrne | Modelling solid tumour growth using the theory of mixtures[END_REF][START_REF] Chatelain | Morphological changes in early melanoma development: influence of nutrients, growth inhibitors and cell-adhesion mechanisms[END_REF]. Indeed, when the distance between cells falls below a certain value (i.e. if the cell density is large enough), cells are attracted to each other. Then, it exists a threshold value called the mechanical equilibrium for which ρH(c e ) + Q(c e ) = 0 i.e. there is an equilibrium between attractive and repulsive forces. For larger cell densities, cells are packed too close to each other, they thus experience a repulsive force. When cells are so packed that they fill the whole control volume, then the repulsive force becomes infinite due to the pressure. The representation of such functional is depicted in Figure 7b. A typical example of single-well potential which has been used for the modeling of living tissue and cancer [START_REF] Agosti | A Cahn-Hilliard-type equation with application to tumor growth dynamics[END_REF][START_REF] Chatelain | Morphological changes in early melanoma development: influence of nutrients, growth inhibitors and cell-adhesion mechanisms[END_REF] is where k is an arbitrary constant.

ψ

Proposition 3 . 5 .

 35 Let ρ be given such that ρ ∈ C(0, T ; C 2 (Ω)) and ρ ≥ ρ > 0. Then Equation (3.13) with Neumann boundary conditions (3.14) admits a strong solution. Moreover, the mapping v n → c[v n ] takes bounded sets in the space C(0, T ; C 2 0 (Ω, R 3 )) into bounded sets in the space

  [START_REF] Cherfils | The compressible Navier-Stokes-Cahn-Hilliard equations with dynamic boundary conditions[END_REF]) is satisfied by our numerical solution, i.e. the derivative of the discrete energy remains negative indicating monotonic decay of the energy. The total mass Ω ρ dx is conserved throughout the simulation as observed in Fig 2b.The scalar variable ξ remains very close to 1 as depicted in Fig 2c.Finally, the bounds of the mass fraction are ensured as seen in Fig.2d

3 Figure 1 :

 31 Figure 1: Spatio-temporal evolution of density ρ, mass fraction c, velocity v and pressure p for test case 1.

Figure 2 :

 2 Figure 2: Temporal evolution of the dissipation of the energy dE dt , of the total mass Ω ρ dx, of the scalar variable ξ, and of the minimum and maximal values of the mass fraction c for test case 1.

Figure 3 :

 3 Figure 3: Convergence in space for the total density ρ, the mass fraction c and the velocity v. The orange dashed line represents the slope 1.

Figure 4 :

 4 Figure 4: Convergence in time for the total density ρ, the mass fraction c and the velocity v. The orange dashed line represents the slope 1.

  A.11) where T denotes the temperature. From the mass balance equation (A.5), we have the entropy balance equation

  ) T ∇c, Dρ Dt = -ρdiv(v) + F ρ ,in the entropy inequality (A.15), we obtain

-D iss := -ρ 2 -

 2 (κ(ρ, c) + F ρ ) |v| 2 -F v v -1 2 |v| 2 + u -ρ ∂F ∂ρ F ρ + c ρ F ρ + c c F c + c v F v -L c ρ Dc Dt -div (J 1 ) -F c -cF ρ ≤ 0.(A.[START_REF] Chatelain | Emergence of microstructural patterns in skin cancer: a phase separation analysis in a binary mixture[END_REF] with p defined in (A.24).B Model reduction, general assumptions and biologically relevant choices of the model's functionsB.1 Specific choices of functionals and model reductions Problem 1: General compressible NSCH with friction term and mass transfer. Assuming no creation of mass nor transfer of mass from the exterior of the system we haveF c = -F 1-c , (B.1)leading to mass conservationF ρ = 0. (B.2)Furthermore, we assume no external source of velocity and energy, leading toF v =0, and F u = 0. (B.3)

2 |∇c| 2 1 3 ∇(B. 7 )

 2137 ρcv) = div (b(c)∇µ) + F c , ρv ⊗ v) = -∇p + γdiv ∇c ⊗ ∇c -1 + div ν(c) ∇v + ∇v T -2 (ν(c) (div (v))) -κ(ρ, c)v, Remark B.1.In this article, we prove the existence of a weak solutions for Problem 1 and propose an efficient structure and bounds-preserving scheme.

F(B. 9 )B. 2

 92 c = F ρ = ρcP c (p), and F 1-c = 0, with a pressure-dependent proliferation rate P c (p) ≥ 0. The growth function P c (p) is used to represent the capacity of cells to divide accordingly to the pressure exerted on them. It is well known that cells are able to divide as long as the pressure is not too large. Once a certain pressure p max is reached cells enter a quiescent state. Therefore, we assume thatP ′ c (p) ≤ 0, and P c (p) = 0 for p > p max . (B.8)Combining these changes, the model becomes (ρv) = ρcP c (p),∂ρc ∂t + div (ρcv) = div (b(c)∇µ) + ρcP c (p), ρµ = -γ∆c + ρ ∂ψ 0 ∂c , ∇p + κ(ρ, c)v = -γdiv ∇c ⊗ ∇c -1 2 |∇c| 2 1 -ρcP c (p)v.Biologically consistent choices of functionsAs said in the derivation of the model, the free energy density F is the sum of two terms: γ 2 |∇c| 2 taking into account the surface tension effects existing between the phases of the mixture and the potential ψ 0 (ρ, c) representing the cell-cell interactions and pressure. Thus, we chooseψ 0 (ρ, c) = ψ e (ρ) + ψ mix (ρ, c), (B.10) with ψ mix (ρ, c) = H(c) log ρ + Q(c).Then, using the constitutive relation for the pressure we have p(ρ, c) = ρ 2 ∂ψ 0 ∂ρ = p e (ρ) + ρH(c). (B.11) The function b(ρ) is the active mobility of the cells of the growing population. Let us explain how the choices of functions for the free energy density and mobility are motivated by biological observations. To satisfy the conditions (2.7), we propose to choose b(c) = C b c(1 -c) α , α ≥ 1, (B.12)

2 (α 1

 21 (1 -c) log(ρ(1 -c)) + α 2 c log(ρc)) -

Figure 7 :

 7 Figure 7: For both figures ρ = 1. Double-well logarithmic potential (left) with α 1 = 1.2 and α 2 = 0.8 and single-well logarithmic potential (right)

3 -( 1 -c e ) c 2 2 -( 1 - 3 -( 1 -c e ) c 2 2 -( 1 -

 31213121 mix (ρ, c) = -(1 -c e ) log(ρ(1 -c)) -c 3 c e )c + k, (B.15)thus givingH(c) = -(1 -c e ), Q(c) = -(1 -c e ) log(1 -c) -c 3 c e )c + k, (B.16)

  terms are taken at the previous time step to linearize the equations). The solution v n+1 , µ n+1 of these equations, together with the array λ, is used in Equation (4.13) to find r n+1 and, then, in Equation (4.14), ξ n+1 . At this point, we solve Equation (4.15) and(4.16) from the solution of the previous steps.

	Remark 4.3. To obtain the interface values U * j+ 1 2

.16) Remark 4.2. We emphasize that the scheme (4.4)-(4.16) is linear. Indeed, Equations (4.4) to (4.7) are obviously linear. Then, the coupling between Equation (4.8) and Equation (4.10) is also linear (nonlinear

∆t 2 2∆x 2 ⟨∆ 0,x U n , δ 2 x V ⋆ ⟩.

Acknowledgements

The authors would like to thank Tommaso Lorenzi for his comments concerning the derivation of our G-NSCH model and for the very interesting discussions we had about the modelling of tumor growth. The authors would like also to thank Alain Miranville for fruitful discussions about the compressible NAvier-Stokes-Cahn-Hilliard model and diphasic fluid dynamics.

Proposition 3.13. We have the following estimates uniformly in n and ε:

(A8) {c n } in L ∞ (0, T ; H 1 (Ω)), (A9) {ρ n ∂ c ψ 0 } in L ∞ (0, T ; L r (Ω)) for r < 6a 6+a , (A10) {µ n } in L 2 (0, T ; H 1 (Ω)), (A11) {ρ n µ n } in L 2 (0, T ; L 6a/(6+a) ), (A12) {c n } in L 2 (0, T ; W 2,r (Ω)) ∩ L 2+ν (0, T ; W 1,2+ν ) for some ν > 0,

Proof. Estimates (A1)-(A2)-(A3)-(A4)-(A5) follow immediately from the energy equality (3.26). Estimate (A6) is the result of Lemma 3.7 and estimates (A2)-(A3)-(A4). To obtain estimate (A7), we multiply Equation (3.11) by ρ n , and using integration by parts, we obtain

Using (A2) and (A6), we deduce (A7). To prove Estimate (A8), we first notice that equality (3.26) provides the uniform bound on {∇c n } in L 2 (0, T ; L 2 (Ω)). To conclude with Lemma 3.6, we need to bound Ω ρ n c n . Combining Equations (3.11)-(3.13), we obtain

A Derivation of the model

In this Appendix, we present the rigourous derivation of our G-NSCH model.

A.1 Notation and definitions

We formulate our problem in Eulerian coordinates and in a smooth bounded domain Ω ⊂ R d (where d = {1, 2, 3} is the dimension). The balance laws derived in the following sections are in We have two cell populations in the model where ρ 1 , ρ 2 are the relative densities of respectively populations 1 and 2. Thus, ρ i represents the mass of the population M i per volume occupied by the i-th phase V i , i.e.

Then, we define the volume fractions φ 1 , φ 2 which are defined by the volume occupied by the i-th phase over the total volume of the mixture

Therefore, the mass density of population i which is the mass of population i in volume V is given by

We further assume that the fluid is saturated, i.e.

The total density of the mixture is then given by We also introduce the mass fractions c i = M i /M and we have the relations

We denote by p the pressure inside the mixture and v 1 , v 2 are the velocities of the different phases. We use a mass-average mixture velocity

We define the material derivative for a generic function g (scalar or vector-valued) by

and indicate the definition of the differential operator

In the following, we denote vectors with bold roman letters and we use bold Greek letters to denote second-order tensors.

Since,

we reorganize the terms of (A.18) to obtain

(A.19)

A.6 Constitutive assumptions and model equations

First of all, we assume that the free energy density F is of Ginzburg-Landau type and has the following form [START_REF] Cahn | On spinodal decomposition[END_REF][START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF] F(ρ, c, ∇c)

where ψ 0 is the homogeneous free energy accounting for the processes of phase separation and the gradient term γ 2 |∇c| 2 represents the surface tension between the two phases. This free energy is the basis of the Cahn-Hilliard model which describes the phase separation occurring in binary mixtures. Furthermore, as obtained in Wise et al. [START_REF] Wise | Three-dimensional multispecies nonlinear tumor growth-I: Model and numerical method[END_REF], the adhesion energy between different cell species is indeed well represented by such a choice of the free energy functional. To satisfy the inequality (A. [START_REF] Chatelain | Morphological changes in early melanoma development: influence of nutrients, growth inhibitors and cell-adhesion mechanisms[END_REF]), we first choose

Then, we define the chemical potential µ(ρ, c, ∇c) by

which in turn gives a condition for the Lagrange multiplier

Using these previous constitutive relations, we have already canceled some terms in the entropy inequality

Then, using classical results on isothermal diffusion [START_REF] Eck | Mathematical modeling[END_REF][START_REF] Lowengrub | Quasi-incompressible Cahn-Hilliard fluids and topological transitions[END_REF], we have

and using a generalized Fick's law, we have

where b(c) is a nonnegative mobility function that we will specify in the following. The two constitutive relations for the diffusive fluxes (A.22) and (A.23) together with (A.21), we obtain

Following [START_REF] Abels | On a diffuse interface model for a two-phase flow of compressible viscous fluids[END_REF][START_REF] Lowengrub | Quasi-incompressible Cahn-Hilliard fluids and topological transitions[END_REF], we define the pressure inside the mixture by

From standard rheology, we assume that the fluid meets Newton's rheological laws. The stress tensor is composed of two parts for the viscous P and non-viscous P contributions of stress σ σ σ := P + P, (A. [START_REF] Feireisl | Dynamics of viscous compressible fluids[END_REF] and we have by standard continuum mechanics (see e.g. [START_REF] Abels | On a diffuse interface model for a two-phase flow of compressible viscous fluids[END_REF][START_REF] Anderson | Diffuse-interface methods in fluid mechanics[END_REF][START_REF] Eck | Mathematical modeling[END_REF])

(A.26)

The second term in the non-viscous part of the stress (namely -γ (ρ∇c ⊗ ∇c)) is representing capillary stresses that act at the interface of the two populations. Furthermore, we assume that the bulk viscosity is zero and, consequently, we set λ(c) = -2 3 ν(c). This form for the stress tensor is also the same used for Navier-Stokes fluids [START_REF] Ericksen | Liquid crystals with variable degree of orientation[END_REF]. Using (A.26), we can cancel a new term in (A. [START_REF] Chatelain | Morphological changes in early melanoma development: influence of nutrients, growth inhibitors and cell-adhesion mechanisms[END_REF])

Therefore, the remaining terms of the entropy inequality are the ones associated with proliferation and friction. The last step to satisfy the entropy inequality is to choose arbitrarily a value for the Lagrange multiplier c ρ , such that

Reorganizing the terms we have

The obvious choices are of course

From the previous constitutive relations and choices of Lagrange multipliers, we have that the dissipation inequality (A.19) is satisfied.

A.7 Summary of the model's equations