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

Charles Elbar, Alexandre Poulain

► To cite this version:

Charles Elbar, Alexandre Poulain. Analysis and numerical simulation of a generalized compressible Cahn-Hilliard-Navier-Stokes model with friction effects. ESAIM: Mathematical Modelling and Numerical Analysis, 2024, 58 (5), pp.1989-2034. 10.1051/m2an/2024063. hal-04747238

HAL Id: hal-04747238

https://hal.science/hal-04747238v1

Submitted on 21 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ANALYSIS AND NUMERICAL SIMULATION OF A GENERALIZED COMPRESSIBLE CAHN–HILLIARD–NAVIER–STOKES MODEL WITH FRICTION EFFECTS

CHARLES ELBAR<fnref n="1" />,* AND ALEXANDRE POULAIN<fnref n="2" />

**Abstract.** We propose a new generalized compressible diphasic Navier–Stokes Cahn–Hilliard model that we name G-NSCH. This new G-NSCH model takes into account important properties of diphasic compressible fluids such as possible non-matching densities and contrast in mechanical properties (viscosity, friction) between the two phases of the fluid. The model also comprises a term to account for possible exchange of mass between the two phases. Our G-NSCH system is derived rigorously and satisfies basic mechanics of fluids and thermodynamics of particles. Under some simplifying assumptions, we prove the existence of global weak solutions. We also propose a structure preserving numerical scheme based on the scalar auxiliary variable method to simulate our system and present some numerical simulations validating the properties of the numerical scheme and illustrating the solutions of the G-NSCH model.

**Mathematics Subject Classification.** 35B40, 35B45, 35G20, 35Q35, 35Q92, 65M08.

Received May 11, 2023. Accepted July 30, 2024.

## 1. INTRODUCTION

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

$$ \frac{\partial \rho}{\partial t} + \text{div}(\rho \mathbf{v}) = 0, \qquad (1.1) $$

$$ \frac{\partial(\rho c)}{\partial t} + \text{div}(\rho c \mathbf{v}) = \text{div}(b(c)\nabla\mu) + F_c, \qquad (1.2) $$

$$ \rho\mu = -\gamma\Delta c + \rho\frac{\partial\psi_0}{\partial c}, \qquad (1.3) $$

$$ \frac{\partial(\rho\mathbf{v})}{\partial t} + \text{div}(\rho\mathbf{v} \otimes \mathbf{v}) = -\left[\nabla p + \gamma\text{div}\left(\nabla c \otimes \nabla c - \frac{1}{2}|\nabla c|^2\mathbf{I}\right)\right] + \text{div}(\nu(c)(\nabla\mathbf{v} + \nabla\mathbf{v}^T)) $$

*Keywords and phrases.* Cahn–Hilliard equation, Navier–Stokes equation, asymptotic analysis, mathematical modeling, numerical simulations, scalar auxiliary variable method.

<footnote n="1">¹ Sorbonne Université, CNRS, Université de Paris, Inria, Laboratoire Jacques-Louis Lions (LJLL), 75005 Paris, France.</footnote>

<footnote n="2">² Univ. Lille, CNRS, UMR 8524 - Laboratoire Paul Painlevé, 59000 Lille, France.</footnote>

<footnote>* Corresponding author: charles.elbar@sorbonne-universite.fr</footnote>
$$ -\frac{2}{3}\mathrm{div}(\nu(c)\mathrm{div}(\mathbf{v})\mathbb{I}) + \mathrm{div}(\eta(c)\mathrm{div}(\mathbf{v})\mathbb{I}) - \kappa(\rho, c)\mathbf{v}, \qquad (1.4) $$

stated in $(0, T) \times \Omega$, where $T > 0$ is finite time horizon, and $\Omega \subset \mathbb{R}^d$ ($d = 1, 2, 3$) is an open bounded domain with a smooth boundary $\partial\Omega$.

Interested by the modeling of invasive growth of tumors in healthy tissues, we motivate the different terms of the model with this biological application in mind. However, we emphasize that the model is a general compressible diphasic fluid model that could be used for other applications.

System (1.1)–(1.4) models the motion of a diphasic fluid composed of two immiscible components, i.e. two different cell types (e.g. tumor and healthy cells), and comprises viscosity effects, surface tension, and friction on rigid fibers representing the extracellular matrix (ECM in short). In System (1.1)–(1.4), $\rho$ 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), $\mathbf{v}$ is the mass averaged total velocity, $\mu$ is called the chemical potential, $p$ is the pressure. The coefficient $\gamma$ 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 $\kappa(\rho, c)$ is a non-negative function of the density and the mass fraction, and takes into account the possible difference of friction strength between the two populations. We use this friction term to model possible adhesive effects of the cells on the ECM. The coefficients $\nu(c)$ and $\eta(c)$ represents the viscosity coefficients (shear and dilatational, respectively) of the mixture. Possible differences in viscosities could be considered for the two populations. The function $\psi_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(\cdot)$ accounts for the possible proliferation and death of cells and these two effects are assumed to be modelled as an exchange of mass between the populations. The non-negative function $b(\cdot)$ models the mobility of cells. This function 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, as explained in Appendix B, under suitable choices of functionals, 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 and the rest of tissue (ECM, extracellular fluid, etc.). The proliferation of cells happens by consuming mass from the other phase (we are not injecting mass in the system). Both cell populations move in an ECM constituted of rigid fibers on which they can adhere. As we only focus on the mechanical effects generated by the properties of the cells that could play a role during invasion, we do not consider in the model other effects that are known to be important in tumor progression: e.g. angiogenesis, digestion of the ECM by proteolytic enzymes, role of helping cells located in the stroma.

We emphasize that this article only concerns the analysis and the numerical simulation of the G-NSCH model (1.1)–(1.4). This latter comprises effects that are negligible in biological situations, e.g. inertia effects. We propose here an analysis of the model and a structure preserving numerical scheme for the G-NSCH model. **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 Navier–Stokes–Cahn–Hilliard model (NSCH model in short). The well-known incompressible variant of the compressible NSCH model has been denominated model H (see e.g. [35, 38]). Model H has been proposed to represent viscous fluid flow in an incompressible binary mixture undergoing phase separation. This model assumes matching densities, i.e. $\rho_1 = \rho_2$ and, hence, constant total density $\rho$. To consider non-matching densities, Lowengrub and Truskinovsky [51] proposed the compressible NSCH model. Expanding the divergence term in the mass balance equation, the authors found a relation denoting the quasi-compressible nature of the fluid. Concomitantly, Anderson et al. [10] proposed a similar system. In the present work, we use a similar system. We also remark that a very recent work [63] 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 modelling framework by specific constitutive hypotheses.
Under some simplifying assumptions compared to the system proposed in [51] but being closer to the system in [10], the analysis of the compressible NSCH model with no-flux boundary conditions has been realized by Abels and Feireisl [3]. Their analysis requires to simplify the model proposed in [51] 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 [2]. Considering the same assumptions and dynamic boundary conditions, Cherfils et al. [21] 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 and a physically relevant choice of potential, the well-posedness and regularity analysis of model H has been performed by Abels [1] 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 [4]). For this reason, Abels, Depner and Garcke studied model H with a degenerate mobility [5]. 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 review partially the extensive literature about the Cahn-Hilliard equation 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 [17, 18]. 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 [53]. We are interested here in the application of the Cahn-Hilliard framework to tumor modelling (see e.g. [52]). Latter, different variants of the Cahn-Hilliard model appeared: e.g. (without giving a complete overview) its coupling to Darcy’s law [31], Brinkman’s law [22], chemotaxis [59]. Recently, a variant of the CH equation has been used to better represent the growth and organization of tumors. The main change is the use of a single-well logarithmic degenerate potential instead of a double-well potential [7, 19, 57]. This type of potential has been proposed in [9] to represent the action of the cells depending only on the local density, i.e. attraction at low cell density and repulsion for large cell density representing the tendency of cells to avoid overcrowding. The Cahn-Hilliard framework has also been utilized in systems representing invasive growth of tumors. The interested reader can find a lot of information about phase-field type systems modelling tumor growth and invasion in the very recent survey paper [30] and references therein.

We now review some of the literature about the numerical simulation of NSCH models. The numerical simulation of model H for binary fluids with non-matching densities has been the subject of numerous works (see e.g. [39] 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 [51] has been proposed in [33]. Around the same time, Giesselmann and Pryer [8, 32] designed a discontinuous Galerkin finite element scheme to simulate the quasi-incompressible NSCH system which preserves the total mass and the dissipation of energy. A numerical method has also been proposed in [37] in the case of constant mobility $b(c)$ and smooth polynomial potential $\psi(c)$. However, the system simulated in [37] is a simplification of the compressible NSCH system since the derivative pressure does not appear in the definition of the chemical potential $\mu$ 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 mobility combined with a logarithmic potential. We now review briefly some relevant discretization methods for the Cahn-Hilliard equation with degenerate mobility and singular
potentials. Considering a degenerate mobility and a double-well logarithmic potential, we mention the work of Barrett et al. [12]. 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. [7] proposed a similar finite element scheme for the single-well logarithmic potential case. The difficulty in this latter case lies in the fact that the degeneracy and the singularity sets do not coincide and, considering an order parameter that must remain within the bounds [0, 1), negative solutions can appear if a standard discretization method is used. The numerical scheme proposed in [7] solves this issue but does not preserve the mass. In a more recent work, Agosti [6] 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: they solve a strongly coupled nonlinear system and resort to the use of iterative algorithms.

Since the Cahn–Hilliard equation is a gradient flow (see e.g. [49]), a structure-preserving linear scheme can be constructed using the Scalar Auxiliary Variable (SAV in short) method [60]. The SAV method is a very powerful tool to design unconditionally energy-stable numerical schemes for models possessing a gradient-flow (see e.g. [62, 66] and references therein) or Hamiltonian structure (see e.g. [11, 58] and references therein). The SAV method has evolved during the past 6 years starting from the original SAV method [60, 61] to improved variants such as the generalized version GSAV (see e.g. [41, 66]) and the relaxed RSAV method [43]. In our work, we use the GSAV method that has been already used in [42] for the Cahn–Hilliard equation. 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 $\mathbb{R}$ to the physical relevant interval ((0, 1) in the case of a double-well potential). Hence, compared to other techniques, the SAV method has the advantage to allow for the design of a linear, efficient, structure-preserving scheme and can easily be used for our G-NSCH system. We also emphasize that the SAV method has been used for the simulation of the incompressible NSCH model with positive mobility and polynomial potential in [47]. In the present work, we use the GSAV method to design a numerical scheme for the G-NSCH model. Our numerical scheme allows to use degenerate mobility and singular potential functionals which are more physically relevant. To the best of our knowledge, our G-NSCH model is new because it comprises the friction force term and the exchange between the two phases of the fluid. Moreover, the use of the GSAV method for a compressible NSCH system is new, especially with our choice of functionals (i.e. degenerate mobility and singular potential).

**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 notations, functional spaces and assumptions we use in our work for the analytical part but also for the numerical part. Section 3 concerns 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 GSAV 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 the 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.

## 2. GENERAL ASSUMPTIONS, NOTATIONS AND FUNCTIONAL SETTING

The equations are set in a domain $\Omega_T = \Omega \times (0, T)$ with $\Omega$ an open and bounded subset of $\mathbb{R}^d$ ($d = 1, 2, 3$). We assume that the boundary $\partial\Omega$ is sufficiently smooth. We indicate the usual Lebesgue and Sobolev spaces by respectively $L^p(\Omega)$, $W^{m,p}(\Omega)$ with $H^m(\Omega) := W^{m,2}(\Omega)$, where $1 \le p \le +\infty$ and $m \in \mathbb{N}$. For $q \in [1, +\infty]$, 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.

## 2.1. 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 for 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 $\mu$ directly. The case of degenerate mobility, see for instance [24], 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 $\nu(c)$, $\eta(c)$ and permeability $\kappa(\rho, c)$ coefficients are smooth non-negative functions. The mobility is a non-negative function of the order parameter (mass fraction) $c$. Hence, we assume that

$$b \in C^1([0, 1]; \mathbb{R}^+), \quad \text{and} \quad b(c) \ge 0 \quad \text{for} \quad 0 \le c \le 1. \qquad (2.1)$$

In agreement with the literature (see e.g. [21]), the homogeneous free energy $\psi_0(\rho, c)$ is assumed to be of the form

$$\psi_0(\rho, c) = \psi_e(\rho) + \psi_{\text{mix}}(\rho, c), \qquad (2.2)$$

with $\psi_{\text{mix}}(\rho, c) = H(c) \log \rho + Q(c)$ and $Q(c)$ is a double-well (or single-well) potential. Then, using the constitutive relation for the pressure, we have

$$p(\rho, c) = \rho^2 \frac{\partial \psi_0}{\partial \rho} = p_e(\rho) + \rho H(c), \qquad (2.3)$$

where $p_e = \rho^2 \psi'_e(\rho)$ and is assumed to satisfy

$$p_1 \rho^{a-1} - p_2 \le p'_e(\rho) \le p_3(1 + \rho^{a-1}), \quad \text{for } a > 3/2, \quad p_1, p_2, p_3 > 0. \qquad (2.4)$$

We assume that the exchange term $F_c$ (that can depend on the mass fraction and the density) is bounded,

$$|F_c(\rho, c)| + \left| \frac{F_c(\rho, c)}{\rho} \right| \le C, \quad \forall (\rho, c) \in \mathbb{R}^2. \qquad (2.5)$$

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

$$\psi_{\text{mix}} = \frac{1}{2}(\alpha_1(1-c)\log(\rho(1-c)) + \alpha_2 c \log(\rho c)) - \frac{\theta}{2}\left(c - \frac{1}{2}\right)^2. \qquad (2.6)$$

This potential gives

$$H(c) = \frac{1}{2}(\alpha_1(1-c) + \alpha_2 c), \quad Q(c) = \frac{1}{2}(\alpha_1(1-c)\log(1-c) + \alpha_2 c \log(c)) - \frac{\theta}{2}\left(c - \frac{1}{2}\right)^2,$$

where $\theta > 1$.
**Additional assumptions for the existence of weak solutions and analysis of the numerical scheme.**

Concerning the existence of weak solutions and analysis of the numerical scheme, we need to strengthen our assumptions. The viscosity coefficients $\nu(c), \eta(c)$ are assumed to be bounded from below by a positive constant and the friction coefficient $\kappa(\rho, c)$ is assumed to be nonnegative. Moreover, $\nu(c), \eta(c)$ and $\kappa(\rho, c)$ are bounded functions. 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)^\alpha$. However, in the analysis, we consider a non-degenerate mobility by truncating the previous mobility. For instance, using a small parameter $0 < \varepsilon_b << 1$, we approximate the mobility $b(\cdot)$ by

$$b_{\varepsilon_b}(c) = \begin{cases} b(1 - \varepsilon_b), & \text{if } c \ge 1 - \varepsilon_b, \\ b(\varepsilon_b), & \text{if } c \le \varepsilon_b, \\ b(c), & \text{otherwise,} \end{cases}$$

and consider the case of a fixed $\varepsilon_b$. Dropping the $\varepsilon_b$ subscript, we obtain that

$$b \in C^1(\mathbb{R}; \mathbb{R}^+), \quad \text{and} \quad b(c) \ge C > 0 \quad \forall c \in \mathbb{R}. \qquad (2.7)$$

Concerning the functionals appearing in the definition of the free energy $\psi_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 \le H'(c), H(c) \le H_2, \quad c \in \mathbb{R}, \quad H_1, H_2 > 0, \\ Q(c) = \frac{1}{4}c^2(1-c)^2. \qquad (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 $\psi_{mix}$ as in the works [3, 21].

Also, to make the computations simpler, we assume that

- $a > 6$ where $a$ is the pressure exponent,

- $\psi_e(\rho) = \frac{\rho^{a-1}}{a-1}$ and therefore $p_e(\rho) = \rho^a$.

These two assumptions are not necessary and could be removed but simplify the analysis. We refer for instance to [3,27] 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 longer. Note that the assumptions on $\psi_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

$$\left| \rho \frac{\partial \psi_0}{\partial c} \right| \le C \rho \psi_0 + C.$$

Its proof uses the assumption on $H$ and the fact that for $c$ large, $Q'(c) \approx c^3 \le c^4 + 1 \approx Q(c) + 1$.

### 3. 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

$$\mathbf{v} = 0, \quad \frac{\partial c}{\partial \mathbf{n}} = b(c) \frac{\partial \mu}{\partial \mathbf{n}} = 0, \quad \text{on } \partial\Omega, \qquad (3.1)$$

and initial conditions

$$\rho(0, x) = \rho_0 \ge 0 \in L^a(\Omega), \quad c(0, x) = c_0 \in H^1(\Omega) \quad \rho_0 \mathbf{v}(0, x) = \mathbf{m}_0, \text{ with } \frac{|\mathbf{m}_0|^2}{\rho_0} \in L^1(\Omega). \qquad (3.2)$$
Also, we suppose $\rho_0 \neq 0$. In this section we take $d = 3$. The proof of the result is quite long and technical. Therefore, when possible 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 outline of 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 $\epsilon_Q$. Then, for this fixed truncation, we prove the existence of weak solutions using the ideas of [3, 21, 27, 48]. Then, we pass to the limit $\epsilon_Q \to 0$. Namely, recalling that $Q(c) = \frac{1}{4}c^2(1-c)^2$ we first consider $Q_{\epsilon_Q}(c)$ a smooth truncated approximation of $Q$ that satisfies

$$|Q_{\epsilon_Q}|, |Q'_{\epsilon_Q}|, |Q''_{\epsilon_Q}| \le C\left(\frac{1}{\epsilon_Q}\right). \quad (3.3)$$

In the first subsections, we drop the $\epsilon_Q$ notation and work with the regularized problem. We will use the $\epsilon_Q$ notation when we pass to the limit. For the moment, we benefit from the properties of the regularization.

## 3.1. 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

$$\frac{d}{dt}E + D = \int_{\Omega} \mu F_c \, dx, \quad (3.4)$$

where E is the energy, and D is the dissipation defined as

$$E = \int_{\Omega} \rho \frac{|\mathbf{v}|^2}{2} + \rho\psi_0 + \frac{\gamma}{2} |\nabla c|^2 \, d\mathbf{x}, \quad (3.5)$$

$$D = \int_{\Omega} \frac{\nu(c)}{2} \left| \nabla \mathbf{v} + \nabla \mathbf{v}^T - \frac{2}{3} \mathrm{div}(\mathbf{v}) \mathbf{I} \right|^2 + \eta(c) |\mathrm{div}(\mathbf{v}) \mathbf{I}|^2 + b(c) |\nabla \mu|^2 + \kappa(\rho, c) |\mathbf{v}|^2 \, d\mathbf{x}. \quad (3.6)$$

This yields a priori estimates on the solution i.e. there exists a positive constant C such that

$$E(t) + \int_{0}^{t} D(s) \, ds \le C + CE(0).$$

Note that the energy is bounded from below since $\rho \log \rho H(c)$ is bounded from below with (2.8). Also, the purpose of the assumptions $\nu(c)$, $\eta(c)$ and $b(c)$ bounded from below by a positive constant becomes clear, they are crucial to obtain estimates on the $H^1(\Omega)$ norm of $\mu$ and $\mathbf{v}$.

*Proof.* We recall the formula

$$\nabla c \Delta c = \mathrm{div}(\nabla c \otimes \nabla c) - \frac{1}{2} \nabla |\nabla c|^2. \quad (3.7)$$

We denote by $\mathbb{T}$ the tensor

$$\mathbb{T} = \nu(c) \left( \nabla \mathbf{v} + \nabla \mathbf{v}^T - \frac{2}{3} \mathrm{div}(\mathbf{v}) \mathbf{I} \right) + \eta(c) \mathrm{div}(\mathbf{v}) \mathbf{I}. \quad (3.8)$$

Then we multiply equation (1.1) by $|\mathbf{v}|^2$ and sum it with the scalar product of equation (1.4) with $\mathbf{v}$. We obtain
$$ \frac{\partial}{\partial t} \left( \rho \frac{|\mathbf{v}|^2}{2} \right) + \mathrm{div} \left( \frac{1}{2} \rho |\mathbf{v}|^2 \mathbf{v} + p(\rho, c) \mathbf{v} - \mathbb{T} \cdot \mathbf{v} \right) + \mathbb{T} : \nabla \mathbf{v} + \kappa(\rho, c) \mathbf{v}^2 = p(\rho, c) \mathrm{div}(\mathbf{v}) \\ + \gamma \mathrm{div} \left( \frac{1}{2} |\nabla c|^2 \mathbb{I} - (\nabla c \otimes \nabla c) \right) \cdot \mathbf{v}, $$

which is equivalent to

$$ \frac{\partial}{\partial t} \left( \rho \frac{|\mathbf{v}|^2}{2} \right) + \mathrm{div} \left( \frac{1}{2} \rho |\mathbf{v}|^2 \mathbf{v} + p(\rho, c) \mathbf{v} - \mathbb{T} \cdot \mathbf{v} \right) + \mathbb{T} : \nabla \mathbf{v} + \kappa(\rho, c) \mathbf{v}^2 = p(\rho, c) \mathrm{div} \mathbf{v} - \gamma \Delta c \nabla c \cdot \mathbf{v}. \quad (3.9) $$

Then, we multiply equation (1.2) by $\mu$ and obtain using also (1.1)

$$ \rho \mu(\partial_t c + \mathbf{v} \cdot \nabla c) = \mathrm{div}(b(c)\nabla\mu)\mu + \mu F_c. $$

And, using (1.3) we obtain

$$ \rho \frac{\partial \psi_0}{\partial c} (\partial_t c + \mathbf{v} \cdot \nabla c) = \mathrm{div}(b(c)\nabla\mu)\mu + \gamma \Delta c (\partial_t c + \mathbf{v} \cdot \nabla c) + \mu F_c. $$

The previous equation can be rewritten using the chain rule as

$$ \begin{aligned} \partial_t(\rho\psi_0) + \mathrm{div}(\rho\psi_0\mathbf{v}) - \psi_0(\partial_t\rho + \mathrm{div}(\rho\mathbf{v})) - \rho\frac{\partial\psi_0}{\partial\rho}(\partial_t\rho + \mathbf{v}\cdot\nabla\rho) & \\ = \mathrm{div}(b(c)\nabla\mu)\mu + \gamma\Delta c(\partial_t c + \mathbf{v}\cdot\nabla c) + \mu F_c. \end{aligned} $$

We have $\rho \frac{\partial \psi_0}{\partial \rho} (\partial_t \rho + \mathbf{v} \cdot \nabla \rho) = \rho \frac{\partial \psi_0}{\partial \rho} (-\rho \mathrm{div}(\mathbf{v})) = -p \mathrm{div}(\mathbf{v})$ (see Eq. (2.3) for the definition of the pressure).

Moreover, we know that $\Delta c \partial_t c = \mathrm{div}(\partial_t c \nabla c) - \partial_t \left(\frac{|\nabla c|^2}{2}\right)$ and, hence,

$$ \partial_t(\rho\psi_0) + \mathrm{div}(\rho\psi_0\mathbf{v}) + p\mathrm{div}(\mathbf{v}) = \mathrm{div}(b(c)\nabla\mu)\mu + \gamma \left[ \mathrm{div}(\partial_t c\nabla c) - \partial_t \left(\frac{|\nabla c|^2}{2}\right) + \Delta c \mathbf{v} \cdot \nabla c \right] + \mu F_c. \quad (3.10) $$

Summing (3.9) and (3.10) we obtain

$$ \frac{\partial}{\partial t} \left( \rho \frac{|\mathbf{v}|^2}{2} + \rho \psi_0 + \frac{\gamma}{2} |\nabla c|^2 \right) + \mathrm{div} \left( \rho \psi_0 \mathbf{v} + \frac{1}{2} \rho |\mathbf{v}|^2 \mathbf{v} + p(\rho, c) \mathbf{v} - \mathbb{T} : \mathbf{v} - \gamma \partial_t c \nabla c \right) - \mathrm{div}(b(c)\nabla\mu)\mu \\ + \mathbb{T} : \nabla \mathbf{v} + \kappa(\rho, c) |\mathbf{v}|^2 = \mu F_c. $$

Now we use the fact that

$$ \mathbb{T} : \nabla \mathbf{v} = \frac{\nu(c)}{2} \left| \nabla \mathbf{v} + \nabla \mathbf{v}^T - \frac{2}{3} \mathrm{div}(\mathbf{v}) \mathbb{I} \right|^2 + \eta(c) |\mathrm{div}(\mathbf{v}) \mathbb{I}|^2. \quad (3.11) $$

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

$$ \left| \int_{0}^{t} \int_{\Omega} \mu F_c \, dx \, dt \right| \le C \int_{0}^{t} \int_{\Omega} |\mu| \, dx \, dt. $$

We want to use Lemma 3.7 to control the $L^1$ norm of $\mu$. Integrating the equations on $\rho$ to obtain $\int_{\Omega} \rho \, dx = \int_{\Omega} \rho_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 $\rho$ to be in $L^{6/5}$ we have the inequality

$$ \left\| \mathbf{u} - \frac{1}{|\Omega|} \int_{\Omega} \rho \mathbf{u} \right\|_{L^2} \le C \| \nabla \mathbf{u} \|_{L^2} + \| \rho \|_{L^{6/5}}. $$
Using Young's inequality, the fact that in the energy $\rho\psi_0$ contains a term of the form $\rho^{a+1}$ we obtain for $\tilde{C}$ small enough

$$\int_0^t \int_{\Omega} |\mu| \, dx \le C + \tilde{C} \int_0^t \int_{\Omega} |\mu|^2 \, dx \le C + CE(t) + \frac{\inf_c b(c)}{2} \int_{\Omega} |\nabla \mu|^2 \, dx + C \left| \int_{\Omega} \rho \mu \, dx \right|.$$

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 $\rho\mu = \rho \frac{\partial\psi_0}{\partial c} - \gamma\Delta c$. Using the Neumann boundary conditions on $c$, it remains to control $|\int_{\Omega} \rho \frac{\partial\psi_0}{\partial c}|$. Using Lemma 2.2, we obtain

$$\left| \int_{\Omega} \rho \frac{\partial \psi_0}{\partial c} \, dx \right| \le C + CE(t).$$

We conclude using Gronwall's lemma.

□

## 3.2. Existence of weak solutions for fixed $\varepsilon_Q$

The weak solutions of system (1.1)–(1.4) are defined as follows

**Definition 3.2.** We say that $(\rho, \mathbf{v}, c, \mu)$ is a weak of system (1.1)–(1.4) provided:

- $\rho \ge 0$ and we have the regularity

$$\begin{align*}
\rho &\in L^\infty(0, T; L^a(\Omega)), \\
\mathbf{v} &\in L^2(0, T; H_0^1(\mathbb{R}^3)), \quad \sqrt{\rho} \mathbf{v} \in L^\infty(0, T; L^2(\Omega; \mathbb{R}^3)), \quad \mathcal{T} : \nabla \mathbf{v} \in L^1(0, T; L^1(\Omega)), \\
c &\in L^\infty(0, T; H^1(\Omega)), \\
\mu &\in L^2(0, T; H^1(\Omega)).
\end{align*}$$

- Equations (1.1)–(1.4) are satisfied in the distributional sense.

- The initial conditions (3.2) are satisfied a.e. in $\Omega$.

- The boundary conditions (3.1) are satisfied.

We state our main theorem about the existence of weak solutions

**Theorem 3.3 (Existence of weak solutions).** There exist $(\rho, \mathbf{v}, c, \mu)$ weak solutions of (1.1)–(1.4) in the sense of Definition 3.2.

In order to prove the existence of weak solutions, we use an approximating scheme with a small parameter $\varepsilon > 0$ borrowing the idea from [27,48]. More precisely, let $X_n = \operatorname{span}\{\eta_i\}_{i=1,...,n}$ be the set of the first $n$ vectors of a basis of $H_0^1(\Omega; \mathbb{R}^3)$ such that $X_n \subset C^2(\bar{\Omega}; \mathbb{R}^3)$. We consider the following problem for $(\rho, \mathbf{v}_n, c)$ with $\mathbf{v}_n \in X_n$ (with coordinates depending on time):

$$\partial_t \rho + \operatorname{div}(\rho \mathbf{v}_n) = \varepsilon \Delta \rho, \qquad (3.12)$$

and for every $\eta \in X_n$,

$$\begin{align}
\int_{\Omega} \rho \mathbf{v}_n(t) \cdot \eta \, dx - \int_{\Omega} \mathbf{m}_0 \cdot \eta \, dx - \int_0^t \int_{\Omega} \rho \mathbf{v}_n \otimes \mathbf{v}_n : \nabla \eta \, dx \, ds - \int_0^t \int_{\Omega} p(\rho, c) \operatorname{div}(\eta) \, dx \, ds \nonumber \\
+ \varepsilon \int_0^t \int_{\Omega} (\nabla \mathbf{v}_n \cdot \nabla \rho) \cdot \eta \, dx \, ds + \int_0^t \int_{\Omega} \mathcal{T} : \nabla \eta \, dx \, ds + \gamma \int_0^t \int_{\Omega} \left( \frac{1}{2} |\nabla c|^2 \mathbb{I} - (\nabla c \otimes \nabla c) \right) : \nabla \eta \, dx \, ds \nonumber \\
+ \int_{\Omega} \int_0^t \kappa(\rho, c) \mathbf{v}_n \cdot \eta \, dx \, ds = 0. \tag{3.13}
\end{align}$$
And for the equation on the mass fraction

$$ \partial_t c + \mathbf{v}_n \cdot \nabla c = \frac{1}{\rho} \mathrm{div}(b(c)\nabla\mu) + \frac{F_c}{\rho}, \quad \mu = \frac{\partial\psi_0}{\partial c} - \gamma \frac{\Delta c}{\rho}. \qquad (3.14) $$

We consider Neumann boundary conditions

$$ \nabla\rho \cdot \mathbf{n} = b(c)\nabla\mu \cdot \mathbf{n} = \nabla c \cdot \mathbf{n} = 0 \quad \text{on } \partial\Omega, \qquad (3.15) $$

and the Dirichlet boundary condition for $\mathbf{v}_n$ is included in the definition of $X_n$. Finally, we consider the initial conditions

$$ \rho(0, \cdot) = \rho_{0,\varepsilon} > 0, \quad c(0, \cdot) = c_{0,\varepsilon}, \quad \rho \mathbf{v}_n(0, \cdot) = \mathbf{m}_0, \qquad (3.16) $$

where $\rho_{0,\varepsilon}$, $c_{0,\varepsilon}$ satisfy the Neumann boundary conditions and they are smooth approximations of $\rho_0$, $c_0$ (when $\varepsilon \to 0$).

We now comment on the scheme used above and detail the strategy of the proof. We add the artificial diffusion in (3.12) with the parameter $\varepsilon > 0$. Here, $\mathbf{v}_n$ is fixed and we can conclude the global in time existence of classical solutions to (3.12) 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.14) which is in fact a fourth-order parabolic equation. Having obtained $c$, we focus on equation (3.13) and we prove existence for a small time with Schauder's fixed point theorem. Note the presence of the additional term $\varepsilon \int (\nabla \mathbf{v}_n \nabla \rho) \cdot \eta$ which is useful to cancel energy terms introduced by $\varepsilon \Delta \rho$ in (3.12). 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 \to \infty$. It remains to send $\varepsilon$ and $\varepsilon_Q$ to 0 and obtain solutions of system (1.1)–(1.4).

We first turn our attention to equation (3.12). From [27], we obtain the following proposition, and lemma

**Proposition 3.4.** Let $\Omega \subset \mathbb{R}^3$ be a bounded domain of class $C^{2+\beta}$ for some $\beta > 0$. For a fixed $\mathbf{v}_n \in X_n$, there exists a unique solution to equation (3.12) with Neumann boundary conditions (3.15) and initial data conditions (3.16). Furthermore, the mapping $\mathbf{v}_n \mapsto \rho[\mathbf{v}_n]$, that assigns to any $\mathbf{v}_n \in X_n$ the unique solution of (3.12), takes bounded sets in the space $C(0, T; C_0^2(\bar{\Omega}, \mathbb{R}^3))$ into bounded sets in the space

$$ V := \{ \partial_t \rho \in C(0, T; C^\beta(\bar{\Omega})), \rho \in C(0, T; C^{2+\beta}(\bar{\Omega})) \}. $$

**Lemma 3.5.** The solutions of (3.12) satisfy

$$ \left( \inf_{x \in \Omega} \rho(0, x) \right) \exp \left( - \int_0^t \| \mathrm{div} \, \mathbf{v}_n(s) \|_{L^\infty(\Omega)} \mathrm{d}s \right) \le \rho(t, x) \\ \le \left( \sup_{x \in \Omega} \rho(0, x) \right) \exp \left( - \int_0^t \| \mathrm{div} \, \mathbf{v}_n(s) \|_{L^\infty(\Omega)} \mathrm{d}s \right), $$

for all $t \in [0, T]$ and $x \in \Omega$.

Using the latter lemma, if the velocity field is in $W^{1,\infty}$, the density is bounded from below by a positive constant (provided the initial condition is positive). We now focus on equation (3.14).

**Proposition 3.6.** Let $\rho$ be given such that $\rho \in C(0, T; C^2(\bar{\Omega}))$ and $\rho \ge \underline{\rho} > 0$. Then equation (3.14) with Neumann boundary conditions (3.15) admits a strong solution. Moreover, the mapping $\mathbf{v}_n \mapsto c[\mathbf{v}_n]$ takes bounded sets in the space $C(0, T; C_0^2(\bar{\Omega}, \mathbb{R}^3))$ into bounded sets in the space

$$ W := \{c \in L^\infty(0, T; H^1(\Omega)) \cap L^2(0, T; H^3(\Omega))\}. \qquad (3.17) $$
The existence of a strong solution is based on the remark that the highest order term of this equation is $-\gamma \frac{b(c)}{\rho} \Delta^2 c$. Using $b(c), \rho \ge C > 0$ we obtain a fourth-order parabolic equation with smooth coefficients and with zero Neumann boundary conditions. Therefore, we can admit the global in time strong solution than can be achieved through a Galerkin scheme and we focus on the estimates (3.17). In the proof, we need the following two lemmas

**Lemma 3.7** ([28], Lem. 3.1). Let $\Omega \in \mathbb{R}^3$ be a bounded Lipschitz domain and let $M_0 > 0, K > 0$. Assume that $\rho$ is a nonnegative function such that

$$0 < M_0 \le \int_{\Omega} \rho dx, \int_{\Omega} \rho^a dx \le K, \quad \text{with } a > \frac{6}{5}.$$

Then, there exists a positive constant $C = C(M_0, K, a)$ such that the inequality

$$\left\| \mathbf{u} - \frac{1}{|\Omega|} \int_{\Omega} \rho \mathbf{u} \right\|_{L^2(\Omega; \mathbb{R}^{3 \times 3})} \le C \| \nabla \mathbf{u} \|_{L^2(\Omega; \mathbb{R}^{3 \times 3})},$$

holds for any $\mathbf{u} \in W^{1,2}(\Omega; \mathbb{R}^3)$.

**Lemma 3.8** ([29], Thm. 10.17). Let $\Omega \subset \mathbb{R}^3$ be a bounded Lipschitz domain, and let $1 < p < +\infty, M_0 > 0, K > 0, a > 1$. Then there exists a positive constant $C = C(p, M_0, K, a)$ such that the inequality

$$\| \mathbf{u} \|_{W^{1,p}(\Omega; \mathbb{R}^3)} \le C \left( \| \nabla \mathbf{u} + \nabla^T \mathbf{u} - \frac{2}{3} \mathrm{div} \mathbf{u} \mathbb{I} \|_{L^p(\Omega; \mathbb{R}^{3 \times 3})} + \int_{\Omega} \rho |\mathbf{u}| \, dx \right),$$

holds for any $\mathbf{u} \in W^{1,p}(\Omega; \mathbb{R}^3)$ and any non-negative function $\rho$ such that

$$0 < M_0 \le \int_{\Omega} \rho dx, \quad \int_{\Omega} \rho^a dx \le K.$$

*Proof of Proposition 3.6.* We admit the existence of solutions and focus on a priori estimates. We multiply equation (3.14) by $-\Delta c$. Using the boundary conditions and integrating in space yields

$$\begin{aligned} \partial_t \int_{\Omega} \frac{|\nabla c|^2}{2} dx + \gamma \int_{\Omega} b(c) \left| \nabla \left( \frac{\Delta c}{\rho} \right) \right|^2 dx \\ = \int_{\Omega} \frac{1}{2} \mathrm{div}(\mathbf{v}_n) |\nabla c|^2 - \nabla \mathbf{v}_n : \nabla c \otimes \nabla c \, dx + \int_{\Omega} b(c) \nabla \left( \frac{\partial \psi_0}{\partial c} \right) \cdot \nabla \left( \frac{\Delta c}{\rho} \right) dx - \int_{\Omega} \frac{F_c}{\rho} \Delta c. \end{aligned}$$

Here, we have also used the formula (3.7). We use the $L^\infty$ bounds on $\mathbf{v}_n$, $\mathrm{div}(\mathbf{v}_n)$, $b(c)$, $\rho$ the fact that $\frac{F_c}{\rho}$ is also bounded in $L^\infty$, properties on $\partial_c \psi_0$ (3.3), and obtain

$$\partial_t \int_{\Omega} \frac{|\nabla c|^2}{2} dx + \gamma \int_{\Omega} b(c) \left| \nabla \left( \frac{\Delta c}{\rho} \right) \right|^2 dx \le C \int_{\Omega} |\nabla c|^2 dx + C \int_{\Omega} \left| \nabla \frac{\Delta c}{\rho} \right| dx + C \int_{\Omega} |\Delta c|.$$

We want to control the last term on the right-hand side. We use Lemma 3.7 with $\mathbf{u} = \frac{\Delta c}{\rho}(1, 0, 0)^T$ and obtain, together with Neumann boundary conditions on $c$,

$$\left\| \frac{\Delta c}{\rho} \right\|_{L^2(\Omega)} \le C \left\| \nabla \left( \frac{\Delta c}{\rho} \right) \right\|_{L^2(\Omega; \mathbb{R}^3)}. \quad (3.18)$$

Then, writing $\Delta c = \rho \frac{\Delta c}{\rho}$ and using the $L^\infty$ bound on $\rho$,

$$\int_{\Omega} |\Delta c| \le C \left\| \nabla \left( \frac{\Delta c}{\rho} \right) \right\|_{L^2(\Omega; \mathbb{R}^3)}.$$
Finally, using Young's inequality and Gronwall's lemma, we obtain

$$ \sup_{t \in (0, T)} \int_{\Omega} |\nabla c|^2 dx + \gamma \int_{0}^{T} \int_{\Omega} b(c) \left| \nabla \left( \frac{\Delta c}{\rho} \right) \right|^2 dx \leq C. \quad (3.19) $$

With Lemma 3.7 (and integrating the equation on $\rho c$ using also the boundary conditions) we obtain the bound

$$ c \in L^{\infty}(0, T; H^1(\Omega)) \cap L^2(0, T; H^3(\Omega)). \quad (3.20) $$

□

Having defined $\rho$ and $c$, we now solve equation (3.13) with a fixed point argument. We define the operator

$$ \mathcal{M}[\rho] : X_n \to X_n^*, \quad \langle \mathcal{M}[\rho] \mathbf{v}, \mathbf{w} \rangle := \int_{\Omega} \rho \mathbf{v} \cdot \mathbf{w} \, dx, \quad \mathbf{v}, \mathbf{w} \in X_n. $$

This operator ([27]) $\mathcal{M}[\rho]$ is invertible, and

$$ \|\mathcal{M}^{-1}[\rho]\|_{\mathcal{L}(X_n^*; X_n)} \le \frac{1}{\inf_{\Omega} \rho}, \quad \|\mathcal{M}^{-1}[\rho_1] - \mathcal{M}^{-1}[\rho_2]\|_{\mathcal{L}(X_n^*; X_n)} \le C(n, \rho) \|\rho_1 - \rho_2\|_{L^1(\Omega)}, \quad (3.21) $$

for any $\rho_1, \rho_2 \ge \rho$. Finally, equation (3.13) can be reformulated as

$$ \mathbf{v}_n(t) = \mathcal{M}^{-1}[\rho(t)] \left( \mathbf{m}_0^* + \int_0^t \mathcal{N}[\mathbf{v}_n(s), \rho(s), c(s)] \, ds \right), \quad (3.22) $$

with

$$ \langle \mathbf{m}_0^*, \eta \rangle = \int_{\Omega} \mathbf{m}_0 \cdot \eta \, dx, $$

and

$$ \langle \mathcal{N}[\mathbf{v}_n, \rho, c], \eta \rangle = \int_{\Omega} \left( \rho \mathbf{v}_n \otimes \mathbf{v}_n - \mathbb{T} - \frac{\gamma}{2} |\nabla c|^2 \mathbb{I} + \gamma \nabla c \otimes \nabla c \right) : \nabla \eta + p(\rho, c) \mathrm{div}(\eta) \\ - (\varepsilon \nabla \mathbf{v}_n \nabla \rho + \kappa(\rho, c) \mathbf{v}_n) \cdot \eta \, dx. $$

To prove that equation (3.22) 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.9** (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 \to S$ has at least one fixed point.

With notation of the Lemma 3.9, we call $A$ the operator from equation (3.22) and $S = B(\mathbf{u}_{0,n})$ the unit ball with center $\mathbf{u}_{0,n}$ in $C([0, T]; X_n)$, $\mathbf{u}_{0,n}$ is defined by

$$ \int_{\Omega} \rho_0 \mathbf{u}_{0,n} \cdot \eta \, dx = \int_{\Omega} \mathbf{m}_0 \cdot \eta \, dx, \quad \forall \eta \in X_n. $$

More precisely, we consider

$$ A: S \to C([0, T]; X_n), \\ \mathbf{u} \mapsto \mathcal{M}^{-1}[\rho(t)] \left( \mathbf{m}_0^* + \int_0^t \mathcal{N}[\mathbf{u}(s), \rho(s), c(s)] \, ds \right). $$
**Lemma 3.10.** 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 $\mathbf{m}_0^*$, we need to prove that $\|\mathcal{M}^{-1}[\rho(t)]\int_0^t \mathcal{N}(s)\mathrm{d}s\|_{C(0,T;X_n)} \le 1$. With properties (3.21), it is sufficient to prove that there exists a final time $T$ small enough such that

$$ \left\| \int_0^t N(s) \mathrm{d}s \right\|_{C(0,T;X_n^*)} \le \inf_{\Omega_T} \rho. $$

Note that the infimum of $\rho$ needs to be taken over the set $\Omega_T = (0, T) \times \Omega$ as $\rho$ depends on time. But, since we only consider small times, using Lemma 3.5 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 \le T_0$, $\inf_{\Omega_T} \rho \ge C(T_0)$. We recall that $X_n \subset C^2(\Omega; \mathbb{R}^3)$ is finite-dimensional. With the definition of the tensor $\mathbb{T}$ and the pressure $p(\rho, c)$ given by (3.11)–(2.3) we estimate by Hölder's inequality:

$$ \int_0^t \int_{\Omega} \left( \rho \mathbf{u} \otimes \mathbf{u} - \mathbb{T} - \frac{\gamma}{2} |\nabla c|^2 \mathbb{I} + \gamma \nabla c \otimes \nabla c \right) : \nabla \eta + p(\rho, c) \mathrm{div}(\eta) - (\varepsilon \nabla \mathbf{u} \nabla \rho + \kappa(\rho, c) \mathbf{u}) \cdot \eta \, \mathrm{d}x \, \mathrm{d}s \\ \le C\left(\sqrt{T} + T\right) \left( \|\eta\|_{X_n} + \|\nabla\eta\|_{X_n} \right) \left( \|\rho\|_{L^\infty} \|\mathbf{u}\|_{L^\infty}^2 + C\|\nu(c)\|_{L^2} \|\nabla\mathbf{u}\|_{L^\infty} + C\|\nabla c\|_{L^4}^2 + \|\rho\|_{L^\infty}^a \right. \\ \left. + \|\rho\|_{L^\infty} \|H(c)\|_{L^\infty} + \varepsilon \|\mathbf{u}\|_{X_n} \|\nabla\rho\|_{L^\infty} + \|\mathbf{u}\|_{L^\infty} \|\kappa(\rho, c)\|_{L^2} \right). $$

Using assumptions of the Section 2.1 and Propositions 3.4–3.6, we prove that all the quantities on the right-hand side are bounded, with a bound that may depend on $n$, except $\|\nabla c\|_{L^4}$ which needs an argument. Note that from (3.17), we deduce $\nabla c$ is bounded in $L^2(0, T; H^2(\Omega)) \cap L^\infty(0, T; L^2(\Omega))$ (by a constant which depends on $\rho$, and also on $\|\mathbf{u}\|_{L^\infty}, \|\nabla \mathbf{u}\|_{L^\infty}$). By Sobolev embedding with $d=3$, $\nabla c$ is bounded in $L^2(0, T; L^\infty(\Omega)) \cap L^\infty(0, T; L^2(\Omega))$. By Hölder inequality (or interpolation), we obtain an $L^4(0, T; L^4(\Omega))$ bound: $\|\nabla c\|_{L^4 L^4}^4 \le \|\nabla c\|_{L^\infty L^2}^2 \|\nabla c\|_{L^2 L^\infty}^2$. With the previous estimates, and for $T$ small enough, we obtain the result. $\square$

**Lemma 3.11.** 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' \le t$ the $X_n$ norm of $\mathcal{M}^{-1}[\rho(t)](\mathbf{m}_0^* + \int_0^t \mathcal{N}[\mathbf{u}(s), \rho(s), c(s)]\mathrm{d}s) - \mathcal{M}^{-1}[\rho(t')](\mathbf{m}_0^* + \int_0^{t'} \mathcal{N}[\mathbf{u}(s), \rho(s), c(s)]\mathrm{d}s)$. For simplicity, we write $\mathcal{N}(s) := \mathcal{N}[\mathbf{u}(s), \rho(s), c(s)]$, and rewrite the previous difference as

$$ \mathcal{M}^{-1}[\rho(t) - \rho(t')] \left( m_0^* + \int_0^t \mathcal{N}(s) \mathrm{d}s \right) + \mathcal{M}^{-1}[\rho(t')] \left( m_0^* + \int_{t'}^t \mathcal{N}(s) \mathrm{d}s \right). $$

For the first term, we use (3.21) and the Hölder continuity of $\rho$ given by Proposition 3.4. For the second term, we repeat the computations in the proof of Lemma 3.10. This ends the result. $\square$

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 $\mathbf{v}_n$ independent of $T(n)$.

**Lemma 3.12.** $\mathbf{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.13) in time and take $\eta = \mathbf{v}_n$ as a test function. This yields

$$ \frac{\mathrm{d}}{\mathrm{d}t} \int_{\Omega} \rho \frac{|\mathbf{v}_n|^2}{2} \mathrm{d}x + \frac{1}{2} \int_{\Omega} (\partial_t \rho + \mathrm{div}(\rho \mathbf{v}_n)) |\mathbf{v}_n|^2 \mathrm{d}x - \int_{\Omega} p(\rho, c) \mathrm{div}(\mathbf{v}_n) \mathrm{d}x - \frac{\varepsilon}{2} \int_{\Omega} \Delta \rho |\mathbf{v}_n|^2 \mathrm{d}x $$
$$ + \int_{\Omega} \mathbb{T} : \nabla \mathbf{v}_n \, dx + \gamma \int_{\Omega} \left( \frac{1}{2} |\nabla c|^2 \mathbb{I} - (\nabla c \otimes \nabla c) \right) : \nabla \mathbf{v}_n \, dx + \int_{\Omega} \kappa(\rho, c) |\mathbf{v}_n|^2 \, dx = 0. \quad (3.23) $$

Here we used

$$ \begin{aligned} \int_{\Omega} \partial_t (\rho \mathbf{v}_n) \cdot \mathbf{v}_n &= \frac{1}{2} \frac{d}{dt} \int_{\Omega} \rho |\mathbf{v}_n|^2 \, dx + \frac{1}{2} \int_{\Omega} \partial_t \rho |\mathbf{v}_n|^2 \, dx, \\ \int_{\Omega} \mathrm{div}(\rho \mathbf{v}_n \otimes \mathbf{v}_n) \cdot \mathbf{v}_n \, dx &= \frac{1}{2} \int_{\Omega} \mathrm{div}(\rho \mathbf{v}_n) |\mathbf{v}_n|^2 \, dx, \\ \varepsilon \int_{\Omega} (\nabla \mathbf{v}_n \cdot \nabla \rho) \cdot \mathbf{v}_n \, dx &= -\frac{\varepsilon}{2} \int_{\Omega} \Delta \rho |\mathbf{v}_n|^2 \, dx. \end{aligned} $$

With (3.12), we see that (3.23) reads

$$ \begin{aligned} & \frac{d}{dt} \int_{\Omega} \rho \frac{|\mathbf{v}_n|^2}{2} \, dx - \int_{\Omega} p(\rho, c) \mathrm{div}(\mathbf{v}_n) \, dx + \int_{\Omega} \mathbb{T} : \nabla \mathbf{v}_n \, dx \\ & \qquad + \gamma \int_{\Omega} \left( \frac{1}{2} |\nabla c|^2 \mathbb{I} - (\nabla c \otimes \nabla c) \right) : \nabla \mathbf{v}_n \, dx + \int_{\Omega} \kappa(\rho, c) |\mathbf{v}_n|^2 \, dx = 0. \end{aligned} \quad (3.24) $$

Now as in (3.10), we obtain with the artificial viscosity

$$ \begin{aligned} & \partial_t (\rho \psi_0) + \mathrm{div}(\rho \psi_0 \mathbf{v}_n) + p \mathrm{div}(\mathbf{v}_n) - \psi_0 \varepsilon \Delta \rho - \varepsilon \rho \frac{\partial \psi_0}{\partial \rho} \Delta \rho = \mathrm{div}(b(c) \nabla \mu) \mu \\ & \qquad + \mathrm{div}(\partial_t c \nabla c) - \partial_t \left( \frac{|\nabla c|^2}{2} \right) + \gamma \Delta c \mathbf{v}_n \cdot \nabla c + \mu F_c. \end{aligned} $$

Integrating this equation in space, and summing with (3.24), we obtain

$$ \begin{aligned} & \frac{d}{dt} \int_{\Omega} \rho \left( \frac{|\mathbf{v}_n|^2}{2} + \psi_0 \right) + \gamma \frac{|\nabla c|^2}{2} \, dx + \varepsilon \int_{\Omega} \nabla \left( \psi_0 + \rho \frac{\partial \psi_0}{\partial \rho} \right) \cdot \nabla \rho \, dx \\ & \qquad + \int_{\Omega} \mathbb{T} : \nabla \mathbf{v}_n \, dx + \int_{\Omega} b(c) |\nabla \mu|^2 \, dx + \int_{\Omega} \kappa(\rho, c) |\mathbf{v}_n|^2 \, dx = \int_{\Omega} \mu F_c \, dx. \end{aligned} \quad (3.25) $$

By definition of $\psi_0$, we obtain

$$ \begin{aligned} \varepsilon \int_{\Omega} \nabla \left( \psi_0 + \rho \frac{\partial \psi_0}{\partial \rho} \right) \cdot \nabla \rho \, dx = \varepsilon \int_{\Omega} \left( \left( (a-1) + (a-1)^2 \right) \rho^{a-2} + \frac{H(c)}{\rho} \right) |\nabla \rho|^2 \, dx \\ + \varepsilon \int_{\Omega} \left( H'(c)(\log(\rho) + 1) + Q'(c) \right) \nabla c \cdot \nabla \rho \, dx. \end{aligned} $$

Therefore, the energy reads

$$ \begin{aligned} & \frac{d}{dt} \int_{\Omega} \rho \left( \frac{|\mathbf{v}_n|^2}{2} + \psi_0 \right) + \gamma \frac{|\nabla c|^2}{2} \, dx + \varepsilon \int_{\Omega} \left( \left( (a-1) + (a-1)^2 \right) \rho^{a-2} + \frac{H(c)}{\rho} \right) |\nabla \rho|^2 \, dx \\ & \qquad + \int_{\Omega} \mathbb{T} : \nabla \mathbf{v}_n \, dx + \int_{\Omega} b(c) |\nabla \mu|^2 \, dx + \int_{\Omega} \kappa(\rho, c) |\mathbf{v}_n|^2 \, dx = \int_{\Omega} \mu F_c \, dx \\ & \qquad - \varepsilon \int_{\Omega} \left( H'(c)(\log(\rho) + 1) + Q'(c) \right) \nabla c \cdot \nabla \rho \, dx. \end{aligned} \quad (3.26) $$

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(\rho) + 1)^2$ is bounded by a constant times $\frac{1}{\rho} + (a-1)\rho^{a-2}$ that it can be bounded in terms of the left-hand side. Note that we used the hypothesis $|Q'(c)| \le C$. This is based on the fact that $Q$ is in fact $Q_{\varepsilon_Q}$ so that we have $|Q'(c)| \le C(\frac{1}{\varepsilon_Q})$ with a constant that blows up when $\varepsilon_Q$ is sent to 0. As we intend to send $\varepsilon_Q \to 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 $\varepsilon \int_{\Omega} Q'(c)\nabla c \cdot \nabla \rho dx$ can be estimated by $\frac{\varepsilon}{4} \int_{\Omega} ((a-1) + (a-1)^2)\rho^{a-2} + \frac{H(c)}{\rho}) |\nabla \rho|^2 dx$ and $\int_{\Omega} \varepsilon C(\frac{1}{\varepsilon_Q}) |\nabla c|^2 dx$. Since $\varepsilon$ will be sent to 0 before $\varepsilon_Q$, the energy inequality will still hold independently of $\varepsilon_Q$ in the limit $\varepsilon \to 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(\Omega; \mathbb{R}^3))$ independently of $T(n)$. Also, the previous bounds do not depend on $n$. Since all the norms are equivalent, it is also bounded in $L^1(0, T(n); W^{1,\infty}(\Omega, \mathbb{R}^3))$. Therefore, we can apply the maximum principle stated in Lemma 3.5, and obtain that the density $\rho$ 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(\Omega; \mathbb{R}^3)$. This procedure can be repeated for every final time $T$. $\square$

Finally, we are left with the following proposition

**Proposition 3.13.** For any fixed $n$ and $T$, there exists a solution $(\rho, c, v_n)$ defined on $(0, T)$ (with appropriate regularity) to (3.12)–(3.14) subject to boundary conditions (3.15) and initial conditions (3.12). Moreover, this solution satisfies the energy dissipation inequality

$$ E(t) + \varepsilon \int_{\Omega_t} \left( (a+a^2)\rho^{a-1} + \frac{H(c)}{\rho} \right) |\nabla \rho|^2 dx dt \\ + \int_{\Omega_t} \mathbb{T} : \nabla v_n dx dt + \int_{\Omega_t} b(c) |\nabla \mu|^2 dx dt + \int_{\Omega_t} \kappa(\rho, c) |v_n|^2 dx dt \le C + CE(0), \quad (3.27) $$

where

$$ E(t) = \int_{\Omega} \rho \left( \frac{|v_n|^2}{2} + \psi_0 \right) + \gamma \frac{|\nabla c|^2}{2} dx, $$

and with a constant $C = C(1, \frac{\varepsilon}{\varepsilon_Q})$ that does not depend on $n$.

Now, we need to find estimates, independent of $n$, to pass to the limit $n \to \infty$. Since $\rho$ and $c$ depend on $n$, we write $\rho_n$ and $c_n$ from now on.

**Proposition 3.14.** We have the following estimates uniformly in $n$ and $\varepsilon$:

(A1) $\{\rho_n\}_{n \in \mathbb{N}}$ in $L^\infty(0, T; L^1(\Omega))$,

(A2) $\{\rho_n\}_{n \in \mathbb{N}}$ in $L^\infty(0, T; L^a(\Omega))$,

(A3) $\{\mathbb{T} : \nabla v_n\}_{n \in \mathbb{N}}$ in $L^1(0, T; L^1(\Omega))$,

(A4) $\{\sqrt{\rho_n} v_n\}_{n \in \mathbb{N}}$ in $L^\infty(0, T; L^2(\Omega; \mathbb{R}^3))$,

(A5) $\{\sqrt{b(c_n)} \nabla \mu_n\}_{n \in \mathbb{N}}$ in $L^2(0, T; L^2(\Omega; \mathbb{R}^3))$,

(A6) $\{v_n\}_{n \in \mathbb{N}}$ in $L^2(0, T; H_0^1(\Omega; \mathbb{R}^3))$,

(A7) $\{\sqrt{\varepsilon} \nabla \rho_n\}_{n \in \mathbb{N}}$ in $L^2(0, T; L^2(\Omega))$,

(A8) $\{c_n\}_{n \in \mathbb{N}}$ in $L^\infty(0, T; H^1(\Omega))$,

(A9) $\{\rho_n \partial_c \psi_0\}_{n \in \mathbb{N}}$ in $L^\infty(0, T; L^r(\Omega))$ for $r < \frac{6a}{6+a}$,

(A10) $\{\mu_n\}_{n \in \mathbb{N}}$ in $L^2(0, T; H^1(\Omega))$,

(A11) $\{\rho_n \mu_n\}_{n \in \mathbb{N}}$ in $L^2(0, T; L^{6a/(6+a)})$,

(A12) $\{c_n\}_{n \in \mathbb{N}}$ in $L^2(0, T; W^{2,r}(\Omega)) \cap L^{2+\nu}(0, T; W^{1,2+\nu}(\Omega))$ for some $\nu > 0$,

(A13) $\{\rho_n c_n\}_{n \in \mathbb{N}}$ in $L^\infty(0, T; L^{\frac{6a}{6+a}}(\Omega))$,

(A14) $\{\rho_n c_n v_n\}_{n \in \mathbb{N}}$ in $L^2(0, T; L^{\frac{6a}{3+4a}}(\Omega))$,
(A15) {$p(\rho_n, c_n)$} in $L^{1+\tilde{\nu}}((0,T) \times \Omega))$ for some $\tilde{\nu} > 0$.

*Proof.* Estimates (A1)-(A5) follow immediately from the energy equality (3.27). Estimate (A6) is the result of Lemma 3.8 and estimates (A2)-(A4). To obtain estimate (A7), we multiply equation (3.12) by $\rho_n$, and using integration by parts, we obtain

$$2\varepsilon \int_0^T \int_{\Omega} |\nabla \rho_n|^2 dx dt \le \|\rho_0\|_{L^2(\Omega)}^2 + \|\rho_n\|_{L^\infty(0,T;L^2(\Omega))}^2 + \|\rho_n\|_{L^2(0,T;L^4(\Omega))}^2 \|\nabla \mathbf{v}_n\|_{L^2(0,T;L^2(\Omega)^d)}.$$

Using (A2) and (A6), we deduce (A7). To prove estimate (A8), we first notice that equality (3.27) provides the uniform bound on {$\nabla c_n$} in $L^2(0,T;L^2(\Omega))$. To conclude with Lemma 3.7, we need to bound $\int_{\Omega} \rho_n c_n$. Combining equations (3.12)-(3.14), we obtain

$$\partial_t(\rho_n c_n) + \operatorname{div}(\rho_n c_n \mathbf{v}_n) = -\varepsilon c \Delta \rho + \operatorname{div}(b(c) \nabla \mu) + F_c.$$

Integrating in space, using the boundary conditions, and estimate (A7), the $L^2$ bound on {$\nabla c_n$}, assumption (2.5) yields {$\int_{\Omega} \rho_n c_n$} is in $L^\infty(0,T)$. We deduce estimate (A8). Estimate (A9) follows from the definition of $\psi_0$ and estimate (A1). Estimate (A10) follows from estimates (A5)-(A9) and Lemma 3.7. 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 Section 2.5 of [21]. This estimate is useful when we obtain the convergence a.e. of $\rho_n$ and $c_n$ so we can obtain strong convergence of $p(\rho_n, c_n)$ in $L^1$ by Vitali's convergence theorem. □

From [27], we also obtain the following proposition

**Proposition 3.15.** There exists $r > 1$ and $p > 2$ such that

$$\begin{array}{l} \partial_t \rho_n, \Delta \rho_n \text{ are bounded in } L^r((0,T) \times \Omega), \\ \nabla \rho_n \text{ is bounded in } L^p(0,T; L^2(\Omega, \mathbb{R}^3)), \end{array}$$

independently of $n$ (but not independently of $\varepsilon$).

With all the previous bound, we can pass to the limit when $n \to \infty$ and obtain the different equation and energy estimates in a weak formulation. Since the passage to the limit $n \to \infty$ is simpler than the next passage $\varepsilon \to 0$, we only detail the latter. Indeed, as $n \to \infty$ we can obtain easily strong convergence of $\rho$ 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.14 still hold independently of $\varepsilon$. It remains now to send $\varepsilon$ to 0.

We recall the equations that we want to pass to the limit into:

$$\partial_t \rho_{\varepsilon} + \operatorname{div}(\rho_{\varepsilon} \mathbf{v}_{\varepsilon}) = \varepsilon \Delta \rho_{\varepsilon}, \quad (3.28)$$

$$\partial_t(\rho_{\varepsilon} c_{\varepsilon}) + \operatorname{div}(\rho_{\varepsilon} c_{\varepsilon} \mathbf{v}_{\varepsilon}) = -\varepsilon c_{\varepsilon} \Delta \rho_{\varepsilon} + \operatorname{div}(b(c_{\varepsilon}) \nabla \mu_{\varepsilon}) + F_{c_{\varepsilon}}, \quad (3.29)$$

and for every $\eta$ (sufficiently regular)

$$\begin{aligned} & \int_{\Omega} \rho_{\varepsilon} \mathbf{v}_{\varepsilon}(t) \cdot \eta \,dx - \int_{\Omega} \mathbf{m}_0 \cdot \eta \,dx - \int_{0}^{t} \int_{\Omega} \rho_{\varepsilon} \mathbf{v}_{\varepsilon} \otimes \mathbf{v}_{\varepsilon} : \nabla_{x} \eta \,dx \,ds - \int_{0}^{t} \int_{\Omega} p(\rho_{\varepsilon}, c_{\varepsilon}) \operatorname{div}(\eta) \,dx \,ds \\ & + \varepsilon \int_{0}^{t} \int_{\Omega} (\nabla \mathbf{v}_{\varepsilon} \nabla \rho_{\varepsilon}) \cdot \eta \,dx \,ds + \int_{0}^{t} \int_{\Omega} \mathbb{T}_{\varepsilon} : \nabla \eta \,dx \,ds + \gamma \int_{0}^{t} \int_{\Omega} \left( \frac{1}{2} |\nabla c_{\varepsilon}|^2 \mathbb{I} - (\nabla c_{\varepsilon} \otimes \nabla c_{\varepsilon}) \right) : \nabla \eta \,dx \,ds \\ & + \int_{\Omega} \int_{0}^{t} \kappa(c_{\varepsilon}) \mathbf{v}_{\varepsilon} \cdot \eta \,dx \,ds = 0. \end{aligned} \quad (3.30)$$
Using Proposition 3.14, which yields uniform estimates in $\epsilon$, 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 $\rho_\epsilon$ as we have no estimates on higher order derivatives. We use the following lemma, see [48].

**Lemma 3.16.** Let $g_n$, $h_n$ converge weakly to $g$, $h$ respectively in $L^{p_1}(0, T; L^{p_2}(\Omega))$, $L^{q_1}(0, T; L^{q_2}(\Omega))$ where $1 \le p_1, p_2 \le +\infty$ and

$$ \frac{1}{p_1} + \frac{1}{q_1} = \frac{1}{p_2} + \frac{1}{q_2} = 1. $$

We assume in addition that

$$ \frac{\partial g_n}{\partial t} \quad \text{is bounded in } L^1(0, T; W^{-m,1}(\Omega)) \quad \text{for some } m \ge 0 \text{ independent of } n, \qquad (3.31) $$

and

$$ \|h_n - h_n(t, \cdot + \xi)\|_{L^{q_1}(0, T; L^{q_2}(\Omega))} \to 0 \quad \text{as } |\xi| \to 0, \text{ uniformly in } n. \qquad (3.32) $$

Then, $g_n h_n$ converges to $gh$ in the sense of distributions.

**Remark 3.17.** 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 [54].

We want to apply the previous lemma to the terms $\rho_\epsilon v_\epsilon$, $\rho_\epsilon c_\epsilon$, $\rho_\epsilon \mu_\epsilon$, $\rho_\epsilon c_\epsilon^2$, $\rho_\epsilon v_\epsilon$, $\rho_\epsilon v_\epsilon \otimes v_\epsilon$, $\rho_\epsilon v_\epsilon c_\epsilon$. We admit that $\frac{\partial \rho_\epsilon}{\partial t}$, $\frac{\partial \rho_\epsilon v_\epsilon}{\partial t}$ and $\frac{\partial \rho_\epsilon c_\epsilon}{\partial t}$ satisfy (3.31) by using Proposition 3.14 and equations (3.28)–(3.30). The compactness in space required in (3.32) also uses Proposition 3.14. We refer also to Section 3.1 from [21] for similar results. The terms $\epsilon c_\epsilon \Delta \rho_\epsilon$ and $\epsilon \int_0^t \int_\Omega (\nabla v_\epsilon \nabla \rho) \cdot \eta \, dx \, ds$ converge to 0 (the first one in the distributional sense) thanks to estimates (A7) and (A8).

It remains to pass to the limit in (i.e. identifying the weak limits)

$$ \begin{aligned} & p(\rho_\epsilon, c_\epsilon), \quad \frac{1}{2}|\nabla c_\epsilon|^2, \quad \nabla c_\epsilon \otimes \nabla c_\epsilon, \\ & b(c_\epsilon)\nabla \mu_\epsilon, \quad F_{c_\epsilon}(\rho_\epsilon, c_\epsilon), \quad \rho_\epsilon \partial_c \psi_0. \end{aligned} $$

The convergence of the last term is used to identify $\rho\mu$. To prove the previous convergences, we need to prove strong compactness in $L^2$ of $c_\epsilon$, $\nabla c_\epsilon$ and convergence a.e. of $\rho_\epsilon$ to use Vitali's convergence theorem. But they follow from the arguments in [3] and Sections 3.3 and 3.4 of [21]. We obtain

**Lemma 3.18.** Up to a subsequence (not relabeled),

$$ \rho_\epsilon \to \rho \quad \text{a.e.} \qquad (3.33) $$

$$ c_\epsilon \to c \quad \text{a.e. and strongly in } L^2(0, T; L^2(\Omega)) \qquad (3.34) $$

$$ \nabla c_\epsilon \to \nabla c \quad \text{a.e. and strongly in } L^2(0, T; L^2(\Omega)) \qquad (3.35) $$

Altogether, we can pass to the limit in every term of the equations:

- $p(\rho_\epsilon, c_\epsilon)$: we use (3.33), (3.34) and (A15)

- $\frac{1}{2}|\nabla c_\epsilon|^2$ and $\nabla c_\epsilon \otimes \nabla c_\epsilon$: (3.35)

- $b(c_\epsilon)\nabla \mu_\epsilon$: (3.34) and (A10)

- $F_{c_\epsilon}(\rho_\epsilon, c_\epsilon)$: (3.33), (3.34) and (2.5)

- $\rho_\epsilon \partial_c \psi_0$: (3.33), (3.34) and (A9).

This concludes the argument.
### 3.3. Sending $\varepsilon_Q \to 0$

The last step in our proof is to let $\varepsilon_Q$ vanishes and recover the existence of weak solutions for the double well potential $Q(c) = \frac{1}{4}c^2(1-c)^2$. Since we have the energy estimates from before, that still hold by properties of the weak convergence, 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 $\varepsilon_Q$. We discussed this point after equation (3.26) and, hence, we do not repeat it here. The second point are the estimates obtained in Proposition 3.14. However, the estimates are essentially the same, except for estimate (A9) (that is the only one containing $Q$) which becomes

$$ \{\rho \partial_c \psi_0\} \text{ in } L^\infty(0, T; L^{\frac{2\alpha}{\alpha+2}}(\Omega)). \quad (3.36) $$

This can be proved knowing that, when $\varepsilon_Q \approx 0$, we have that for $c$ large $\rho Q'_{\varepsilon_Q}(c) \approx \rho c^3$, and we use estimates (A2)–(A8). Altogether, the reasoning to pass to the limit is the same and we conclude.

## 4. NUMERICAL SCHEME FOR THE G-NSCH MODEL

We propose a numerical scheme for the G-NSCH model (1.1)–(1.4) subjected to periodic boundary conditions.

We combine ideas from the numerical scheme for the variant of the compressible NSCH system in [37] and fast structure-preserving scheme for degenerate parabolic equations [40, 42]. Namely, we adapt the relaxation [44] of the Navier-stokes part as used in [37]. The part of the scheme for the Cahn-Hilliard part of the system is designed using the GSAV method. More precisely, a variant used for degenerate parabolic models that preserves the physical bounds of the solution [40, 42].

Indeed, we expect that the volume fraction $c$ remains within the physically (or biologically) relevant bounds $c \in (0, 1)$. Thus, following [40, 42], we construct the invertible mapping $T: \mathbb{R} \to (0, 1)$, with $c = T(v)$, transforming equations (1.2)–(1.3) into

$$ \rho(\partial_t v + (\mathbf{v} \cdot \nabla)v) = \frac{1}{T'(v)}(\mathrm{div}(b(c)\nabla\mu) + F_c), \quad (4.1) \\ \rho\mu = -\gamma T'(v)\Delta v - \gamma T''(v)|\nabla v|^2 + \rho \frac{\partial \psi_0}{\partial c}. $$

Following [40, 42], we can choose

$$ T(v) = \frac{1}{2} \tanh(v) + \frac{1}{2}, \quad \text{or} \quad T(v) = \frac{1}{1 + \exp(-v)}, $$

thus, preserving the bounds $c \in (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 $\underline{C}$ such that the energy associated with the Cahn-Hilliard part, i.e.

$$ E[t](\rho, c) = \int_{\Omega} \frac{\gamma}{2} |\nabla c|^2 + \rho \psi_0(\rho, 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 + \underline{C} \ge 1$

We define

$$ r(t) = E(t) + C_0, \quad \text{with} \quad C_0 = 2\underline{C} + \|E(\rho^0, c^0)\|_{L^\infty(\Omega)}, $$

and apply the SAV method. System (4.1) becomes

$$ \rho(\partial_t v + (\mathbf{v} \cdot \nabla)v) = \frac{1}{T'(v)}(\mathrm{div}(b(c)\nabla\mu) + F_c), \quad (4.2) \\ \rho\mu = -\gamma T'(v)\Delta v - \gamma T''(v)|\nabla v|^2 + \rho \frac{\partial \psi_0}{\partial c}, \\ \frac{dr}{dt} = -\frac{r(t)}{E[t] + C_0} \int_{\Omega} b(c)|\nabla\mu|^2 - \mu F_c \,dx, $$
One can easily see that the previous modifications do not change our system at the continuous level.

## 4.1. One-dimensional scheme

We consider our problem in a one-dimensional domain $\Omega = (0, L)$. Even though $\mathbf{v}$ is now a scalar, we still denote it in bold font to not make the confusion with $v$ from the transformation $c = T(v)$. As mentioned previously, we relax the Navier-Stokes part of our system. Namely, we introduce a relaxation parameter $\iota \ge 0$ and write $U = (\rho, \rho\mathbf{v})$. We rewrite equation (1.4) as

$$ \left\{ \begin{array}{l} \partial_t U + \partial_x V = G(U), \\ \partial_t V + A \partial_x U = -\frac{1}{\iota} (V - F(U)), \end{array} \right. \qquad (4.3) $$

in which $G(U) = (0, -\kappa\mathbf{v})$, $F(U) = (\rho\mathbf{v}, \rho\mathbf{v}^2 + p - (\frac{4}{3}\nu(c) + \eta(c))\partial_x\mathbf{v} + \frac{\gamma}{2}|\partial_x c|^2)$ and $A = \text{diag}(a_1, a_2)$ satisfying Liu's subcharacteristic condition

$$ A \geq F'(U), \quad \forall U. $$

In what follows, and following [37], we use

$$ a_1 = a_2 = \max\left\{\sup(\mathbf{v} + \sqrt{\partial_{\rho}p})^2, \sup(\mathbf{v} - \sqrt{\partial_{\rho}p})^2\right\}. $$

We discretize the domain using a set of $N_x$ nodes located at the center of control volumes of size $\Delta x$ such that $\Omega = \cup_{j=0,...,N_x-1} [x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}]$.

Our scheme follows the discrete set of equations

$$ U_j^* = U_j^n, \qquad (4.4) $$

$$ V_j^* = V_j^n - \frac{\Delta t}{\iota} (V_j^* - F(U_j^*)), \qquad (4.5) $$

$$ U_j^{n+1} = U_j^* - \frac{\Delta t}{\Delta x} (V_{j+\frac{1}{2}}^* - V_{j-\frac{1}{2}}^*) + \Delta t G(U_j^{n+1}), \qquad (4.6) $$

$$ V_j^{n+1} = V_j^* - \frac{\Delta t}{\Delta x} A (U_{j+\frac{1}{2}}^* - U_{j-\frac{1}{2}}^*), \qquad (4.7) $$

$$ \rho_j^{n+1} T'(v_j^n) \left( \frac{\bar{v}_j^{n+1} - v_j^n}{\Delta t} + \mathbf{v}_j^{n+1} \cdot (\nabla \bar{v}^{n+1})_j \right) = g(c^n, \mu^{n+1}, \rho^{n+1})_j, \qquad (4.8) $$

$$ g(c^n, \mu^{n+1}, \rho^{n+1})_j = \left( \frac{1}{\Delta x} \left( (b(c^n)\nabla\mu^{n+1})_{j+\frac{1}{2}} - (b(c^n)\nabla\mu^{n+1})_{j-\frac{1}{2}} \right) \right) + F_c(\rho_j^n, c_j^n), \qquad (4.9) $$

$$ \rho_j^{n+1} \mu_j^{n+1} = \left( -\gamma T'(v_j^n) (\Delta \bar{v}^{n+1})_j - \gamma T''(v_j^n) (\nabla v^n)_j \cdot (\nabla \bar{v}^{n+1})_j \right) + \rho_j^{n+1} \left( \frac{\partial \psi_0}{\partial c} \right)_j^n, \qquad (4.10) $$

$$ \sum_j \Delta x T(\lambda \bar{v}_j^{n+1}) = \sum_j \Delta x c^0 + \sum_{r=1}^{n} \Delta t \sum_j \Delta x F_c(\rho_j^r, c_j^r), \qquad (4.11) $$

$$ \bar{c}^{n+1} = T(\lambda \bar{v}^{n+1}), \qquad (4.12) $$

$$ \frac{1}{\Delta t} (r^{n+1} - r^n) = - \frac{r^{n+1}}{E(\bar{c}^{n+1}) + C_0} \Delta x \sum_j b(\bar{c}_j^{n+1}) |(\nabla \mu^{n+1})_j|^2 \\ + \frac{r^{n+1}}{E(\bar{c}^{n+1}) + C_0} \Delta x \sum_j \mu_j^{n+1} F_c(\rho_j^{n+1}, \bar{c}_j^{n+1}), \qquad (4.13) $$

$$ \xi^{n+1} = \frac{r^{n+1}}{E(\bar{c}^{n+1}) + C_0}, \qquad (4.14) $$
$$c_j^{n+1} = \nu^{n+1} \bar{c}_j^{n+1}, \quad \text{with} \quad \nu^{n+1} = 1 - (1 - \xi^{n+1})^2, \qquad (4.15)$$

$$v_j^{n+1} = \lambda \nu^{n+1} \bar{v}_j^{n+1}. \qquad (4.16)$$

**Remark 4.1** (Computation of interface values). To obtain the interface values $U_{j+\frac{1}{2}}^*$, $U_{j-\frac{1}{2}}^*$ and $V_{j+\frac{1}{2}}^*$, $V_{j-\frac{1}{2}}^*$, we use the upwind method, i.e.

$$U_{j+\frac{1}{2}} = \frac{1}{2}(U_j + U_{j+1}) - \frac{\sqrt{a_1}}{2}(V_{j+1} - V_j), \quad V_{j+\frac{1}{2}} = \frac{1}{2}(V_j + V_{j+1}) - \frac{1}{2\sqrt{a_2}}(U_{j+1} - U_j).$$

We also mention that similarly to [37], one can implement a MUSCL scheme (see e.g. [46]) to obtain a higher order reconstruction. The upwind method permits to rewrite equations (4.6) and (4.7) as

$$U_j^{n+1} = U_j^* - \frac{\Delta t}{2\Delta x}(V_{j+1}^* - V_{j-1}^*) + \frac{\Delta t}{2\Delta x}\sqrt{a}(\delta_x^2 U_j^*) + \Delta t G(U_j^{n+1}), \qquad (4.17)$$

$$V_j^{n+1} = V_j^* - \frac{a \Delta t}{2 \Delta x} (U_{j+1}^* - U_{j-1}^*) + \frac{\Delta t}{2 \Delta x} \sqrt{a} (\delta_x^2 V_j^*), \qquad (4.18)$$

where we used the notation $\delta_x^2 U = U_{j+1} - 2U_j + U_{j-1}$. In equations (4.17) and (4.18), we emphasize that $U^* = U^n$ and $V^* = V^n - \frac{\Delta t}{\epsilon}(V^* - F(U^n))$.

**Remark 4.2** (Algorithm to compute the solution of the discrete equations' system). Equations (4.4)-(4.7) are solved from equations (4.17) and (4.18), hence, a solution $(U^{n+1}, V^{n+1})$ is computed just from vector computations. The coupling between equations (4.8) and (4.10) is also linear (nonlinear terms are taken at the previous time step to linearize the equations). We solve this coupled system using the GMRES algorithm but we emphasize that other iterative solver could work as long as they allow the matrix of the linear system to be non-symmetric. The coefficient $\lambda$ is computed using an iterative method. Then, the discrete solution $(\bar{v}^{n+1}, \mu^{n+1})$, together with the coefficient $\lambda$, is used in equation (4.13) to find $r^{n+1}$ and, in equation (4.14), $\xi^{n+1}$. At this point, we solve equations (4.15) and (4.16) from the previous steps.

In the following, we use the notations,

$$\langle U, V \rangle = \Delta x \sum_j U_j V_j, \quad \text{and} \quad \|U\|^2 = \langle U, U \rangle.$$

We also use $\Delta_{0,x} U := \frac{1}{2}(U_{j+1} - U_{j-1})$.

Our numerical scheme possesses the following important properties:

**Proposition 4.3** (Energy stability, bounds and mass preserving). Assuming the CFL-like condition $\frac{\Delta t}{\Delta x} \sqrt{a_1} \le 1$ and the condition

$$\Delta t \le C \frac{C_0}{E[\bar{c}^n]}, \qquad (4.19)$$

our numerical scheme satisfies the energy dissipation-like inequality

$$\|\sqrt{a}U^{n+1}\|^2 + \|V^{n+1}\|^2 + r^{n+1} \le \|\sqrt{a}U^n\|^2 + \|V^*\|^2 + C^{n+1}r^n, \qquad (4.20)$$

where $r^{n+1} \ge 0$ and

$$C^{n+1} = \frac{1}{1 + \frac{\Delta t}{E(\bar{c}^{n+1}) + C_0} \sum_{j=1}^{N_x} b(\bar{c}_j^{n+1}) |\nabla \mu_j^{n+1}|^2 - \mu_j^{n+1} F_c(\rho_j^{n+1}, \bar{c}_j^{n+1})}.$$

The previous constant can be estimated only in terms of $E[\bar{c}^n]$ and therefore do not depend on the step $n + 1$. Furthermore, the numerical scheme preserves the physically relevant bounds of the mass fraction, i.e.

$$0 < c^{n+1} < 1.$$
**Remark 4.4.** 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 is of course satisfied when we have $F_c = 0$ for instance.

*Proof.* We start with equation (4.17), and using the definition of the function $G(U_j^{n+1})$ as well as assuming $\kappa(c) \ge 0$ (for $c \in \mathbb{R}$), after taking the square on both sides, multiplying by $\Delta x$ and summing over the nodes $j = 0, \dots, N_x$, we have

$$ \|U^{n+1}\|^2 \le \|U^n\|^2 + \left(\frac{\Delta t}{2\Delta x}\right)^2 \|\Delta_{0,x}V^*\|^2 + \left(\frac{\Delta t\sqrt{a}}{2\Delta x}\right)^2 \|\delta_x^2 U^n\|^2 - \frac{\Delta t}{\Delta x} \langle \Delta_{0,x}V^*, U^n \rangle \\ + \frac{\Delta t\sqrt{a}}{\Delta x} \langle U^n, \delta_x^2 U^n \rangle - \frac{\sqrt{a}\Delta t^2}{2\Delta x^2} \langle \Delta_{0,x}V^*, \delta_x^2 U^n \rangle. $$

Repeating the same computations for equation (4.18), we have

$$ \|V^{n+1}\|^2 \le \|V^n\|^2 + \left(\frac{a\Delta t}{2\Delta x}\right)^2 \|\Delta_{0,x}U^n\|^2 + \left(\frac{\Delta t\sqrt{a}}{2\Delta x}\right)^2 \|\delta_x^2 V^*\|^2 - \frac{a\Delta t}{\Delta x} \langle \Delta_{0,x}U^n, V^* \rangle \\ + \frac{\Delta t\sqrt{a}}{\Delta x} \langle U^*, \delta_x^2 V^* \rangle - \frac{a^{\frac{3}{2}}\Delta t^2}{2\Delta x^2} \langle \Delta_{0,x}U^n, \delta_x^2 V^* \rangle. $$

At this point, the proof is similar to the proof of Theorem 4.1 from [37] (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$,

$$ \|\sqrt{a}U^{n+1}\|^2 + \|V^{n+1}\|^2 \le C(\|\sqrt{a}U^n\|^2 + \|V^*\|^2). $$

Then, for the Cahn-Hilliard part, we easily obtain from equation (4.13)

$$ r^{n+1} \left( 1 + \Delta t \frac{\Delta x \sum_j b(\bar{c}_j^{n+1}) |(\nabla \mu^{n+1})_j|^2 - \mu_j^{n+1} F_c(\rho_j^{n+1}, \bar{c}_j^{n+1})}{E[\bar{c}^{n+1}] + C_0} \right) = r^n. $$

Therefore, as long as

$$ E(\bar{c}^{n+1}) + C_0 + \Delta t \left( \Delta x \sum_j b(\bar{c}_j^{n+1}) |(\nabla \mu^{n+1})_j|^2 - \mu_j^{n+1} F_c(\rho_j^{n+1}, \bar{c}_j^{n+1}) \right) \ge 0, $$

so does $r^{n+1}$. Assuming $\|F_c\|_{L^\infty} < C$, it remains to control the discrete $L^1$ norm of $\mu^{n+1}$. Performing the same computations as in the proof of Proposition 3.1 in continuous case, it follows that

$$ \left| \Delta x \sum_j \mu_j^{n+1} F_c(\rho_j^{n+1}, \bar{c}_j^{n+1}) \right| \le C + C E[\bar{c}^n] + \frac{1}{2} \Delta x \sum_j b(\bar{c}_j^{n+1}) |(\nabla \mu^{n+1})_j|^2. $$

Of course, one first needs to prove that a discrete version of Lemma 3.7 holds. At the continuous level, this theorem is proved by contradiction using Rellich's theorem. Hence, a similar proof can be obtained at the discrete level, in the spirit of the Poincaré-Wirtinger inequality, see Lemma 3.8, Remark 3.16 from [26]. Based on these evidences we use a discrete version of Lemma 3.7 to conclude that there exists $C$ a universal constant such that provided $r^n \ge 0$ and

$$ \Delta t \le C \frac{C_0}{E[\bar{c}^n]}, $$
so does $r^{n+1} \ge 0$, and (4.20) follows.

Finally, from the definition of $\xi^{n+1}$ and $C_0$, we have

$$0 < \xi^{n+1} < \frac{r^0}{E(c^{n+1}) + C_0} \le 2.$$

The bounds for the mass fraction $c$ are ensured by the transformation $T(v)$. This finishes the proof. $\square$

**Remark 4.5.** We observe during numerical simulations that the condition (4.19) is obtained for reasonably small $\Delta t$. We also note that if we do not consider any source term, i.e. $F_c = 0$, the scheme satisfies the dissipation relation

$$\|\sqrt{a}U^{n+1}\|^2 + \|V^{n+1}\|^2 + r^{n+1} \le C(\|\sqrt{a}U^n\|^2 + \|V^\star\|^2) + r^n,$$

with the stability condition

$$\frac{\Delta t}{\Delta x} \sqrt{a_1} \le 1.$$

## 5. NUMERICAL EXPERIMENTS

In this section, we use the assumptions on the functionals stated in the “Framework for numerical simulations” paragraph in Section 2.1. Throughout this section we use the double-well logarithmic potential

$$\psi_{\text{mix}} = \frac{1}{2}(\alpha_1(1-c)\log(\rho(1-c)) + \alpha_2c\log(\rho c)) - \frac{\theta}{2}\left(c - \frac{1}{2}\right)^2 + k,$$

with $k = 100$, and $\theta = 4$ ($\alpha_1$ and $\alpha_2$ are specified later). We also use a degenerate mobility, i.e.

$$b(c) = c(1-c).$$

We start by using the one-dimensional scheme (4.4)-(4.16) with no exchange term and friction, i.e. $\kappa(\rho, c) = 0$ and $F_c(\rho, c) = 0$, and we verify that the scheme preserves all the properties stated in Proposition 4.3. We then use a non-zero exchange term and we compare the solution with same friction forces for the two phases or contrast of friction forces. Then, we perform two-dimensional simulations with friction forces contrast.

Finally, we verify the spatial and temporal convergence orders of the scheme.

**Remark 5.1 (Implementation details).** All numerical schemes are implemented using Python 3 and the Numpy and Scipy modules. The linear system for the Cahn-Hilliard part of the model is solved using the Generalized Minimal RESidual iteration (GMRES) iterative solver (function available in the `scipy sparse linalg` module). The tolerance on the convergence of the residual is indicated in each of the following subsections. To find the $\lambda$ that allows to compute the correct mass, we use the function `fsolve` of the `scipy_optimize` module which uses a modification of the Powell's conjugate direction method.

### 5.1. One dimensional numerical test cases

**Comparison between matching and non-matching densities.** We start with a one-dimensional test cases to show the spatiotemporal evolution of the density, mass fraction, and velocity. We also verify numerically the properties stated in Proposition 4.3. We compare numerical results for matching and non-matching densities for the phases of the fluid. For this comparison, we set $\kappa(\rho, c) = 0$ and $F_c(\rho, c) = 0$.

We use the computational domain $\Omega = (0, 1)$ discretized in $N_x = 128$ cells. We take $T = 0.5$ (this has been chosen because the system reaches a meta-stable state by that time) and use the initial time step $\Delta t = 1 \times 10^{-5}$ (this time step size is adapted from the CFL-like condition stated in Proposition 4.3).

We choose the width of the diffuse interface to be $\gamma = 1/600$, the viscosity to be constant $\nu(c) = 1 \times 10^{-2}$, and $\eta = 2 \times 10^{-2}$, the relaxation parameter to be $\iota = 1 \times 10^{-5}$, and the exponent for the barotropic pressure equals to $a = 3$.
To model matching densities for the two phases of the fluid, we choose $\alpha_1 = \alpha_2 = 1$. To represent non-matching densities for the two phases, we can choose $\alpha_1 \neq \alpha_2$. 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 $\alpha_1$ and $\alpha_2$ on the potential. Taking $\alpha_1 < \alpha_2$ shifts the well corresponding to phase 1 very close to 0 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 choose constant initial conditions for the density and the pressure, i.e.

$$ \rho_j^0 = 0.8, \quad \mathbf{v}_j^0 = 0.5, \quad j = 0, \dots, N_x - 1. $$

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

$$ c_j^0 = \underline{c} + 0.05r_j, \quad j = 0, \dots, N_x - 1, $$

with $\underline{c} = 0.5$ and $r$ is a vector of random values between 0 and 1 given by the uniform distribution.

We choose a tolerance for the convergence of the residual for the GMRES algorithm of rtol = $10^{-10}$.

Figure 1 compares the results obtained for matching and non-matching densities for the two phases of the fluid. For the two cases, we report the evolution of the density $\rho$, the mass fraction $c$, the velocity $\mathbf{v}$ and the pressure $p$ at different times. We observe that, for both cases, after an initial regularization of the initial condition, the separation of the two phases of the fluid occurs and small aggregates appear (see first and second columns of Fig. 1). Then, the coarsening of the small aggregates into larger ones occurs. We arrive at the end of the simulation to the solution depicted in the two figures on the last column of Figure 1. Hence, we can conclude that our numerical scheme catches well the spinodal decomposition of the binary fluid while it is transported to the right (since the velocity $\mathbf{v}$ is positive during these simulations).

A difference between the two simulations is observed on the densities and pressures. Indeed, for matching densities, we observe that $\rho$ organizes such that it is equal in aggregates of each phases and drops at the interfaces between the aggregates. We also observe a drop of pressure $p$ at the interface, probably explained by capillary effects. For non-matching densities, as expected, there is a density difference between aggregates of phase 1 and 2. Indeed, selecting $\alpha_1 < \alpha_2$ makes the aggregates of phase 1 denser compared to aggregates of phase 2. Our explanation is that, due to the fact that attractive effects are stronger in phase 1, more mass is allowed to move inside aggregates of phase 1. However, as the pressure function accounts for the difference $\alpha_1 \neq \alpha_2$, the pressure equilibrates to a field similar to the matching density case (i.e. $p$ varies from an equilibrium value and depicts a drop at the interface between the aggregates of the different phases).

Figure 2 shows that, for both cases (i.e. matching and non-matching densities), our numerical scheme preserved the properties presented in Proposition 4.3. We defined the discrete dissipation of energy

$$ \frac{dE}{dt} = \|\sqrt{a}U^{n+1}\|^2 + \|V^{n+1}\|^2 + r^{n+1} - \left[ \|\sqrt{a}U^n\|^2 + \|V^*\|^2 + C^{n+1}r^n \right]. \quad (5.1) $$

We emphasize that as no exchange term was present in the previous simulation $C^{n+1}$ is bounded from above by 1, hence we used $C^{n+1} = 1$ for the simulations in this paragraph. Figure 2 presents the temporal evolution of the dissipation $\frac{dE}{dt}$, the mass $\int_{\Omega} \rho c dx$, the minimum and maximum values of $c$, and the value of $\xi$. We observe for both cases that the dissipation (5.1) is strictly negative, as expected by Proposition 4.3. The mass of fluid 1 is preserved up to a small numerical error (we emphasize that the error on the initial mass at the end of the simulation is less than $10^{-9}$ for both simulations). This latter result is expected as we set $F_c(\rho, c) = 0$ for these simulations. We observe that the physical bounds of the mass fraction are ensured, i.e. maximum and minimum values for $c$ lie the interval $(0, 1)$. The scalar variable is very close to 1 (up to an error of order $10^{-5}$) as observed in Figure 2c. This verifies that the modified energy $r^{n+1}$ and the real energy of the Cahn-Hilliard part of the model

$$ E^{n+1} = \Delta x \sum_j \frac{\gamma}{2} |(\nabla c)_j|^2 + \rho^{n+1} \left( \frac{\partial \psi_0}{\partial c^{n+1}} \right)_j, $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Simulation of compressible Navier-Stokes-Cahn-Hilliard model (with $κ = F_c = 0$). Matching densities (top row) and non-matching densities (bottom row) for the two phases of the fluid. (a) $T = 0$. (b) $T = 0.03$. (c) $T = 0.1$. (d) $T = 0.5$. (e) $T = 0$. (f) $T = 0.03$. (g) $T = 0.1$. (h) $T = 0.5$.</figcaption></figure>

are close.

**Mass exchange and contrast of friction forces.** In this test case, we consider mass exchange between the two phases and friction effects. We choose

$$F_c(\rho, c) = r_{\text{trans}}\rho c(1 - c/c_{\text{max}}), \quad \kappa(\rho, c) = \rho c\kappa_1 + \rho(1 - c)\kappa_2,$$

where $0 < c_{\text{max}} < 1$ denotes the mass fraction at which we have an equilibrium for exchange of mass and $r_{\text{trans}}$ is the rate of mass exchange. In this test case, we use $c_{\text{max}} = 0.9$, and $r_{\text{trans}} = 1$. We compare the solution obtained with no contrast of friction effect, i.e. $\kappa_1 = \kappa_2 = 10$, and the solution obtained with $\kappa_1 = 0$, $\kappa_2 = 10$. To study the long time behavior of the numerical simulations, we set $T = 5$. The rest of the parameters and the initial conditions are chosen as for the previous non-matching densities test case.

Figure 3 compares the solutions for the two cases. In both cases, we observe that the separation of the two phases occurs and that the velocity decreases in time due to friction effects. Furthermore, as time passes, phase 1 of the fluid increases due to the exchange term $F_c(\rho, c) \neq 0$ and, hence, zones of mass fraction close to the value 1 enlarge. At the end of both simulations, there is one large aggregate of fluid 1. The difference between the two simulations appears clearly at time $t = 0.1$. When the friction forces are stronger in phase two compared to phase 1, i.e. $\kappa_2 > \kappa_1$, zones of larger density appear destabilized, i.e. the shape of the aggregates is not symmetric (compare the solution for Figs. 3b and 3f). For each aggregate of fluid 1, the density at the right of the aggregate is larger compared to the left. We conclude that the contrast in friction forces is captured well by the model as simulations depict a contrast of velocity for the two phases of the fluid and leading to less regular patterns for the densities. Furthermore, we emphasize that even with non-zero mass exchange and friction forces, the numerical scheme ensures the properties stated in Proposition 4.3 as observed in Figure 4. We emphasize that compared to the simulation without source in which $\xi$ seems to converge to a constant value,
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Temporal evolution of the dissipation of the energy $\frac{dE}{dt}$, mass of the fluid 1 given by $\int_{\Omega} \rho c \, dx$, scalar variable $\xi$, and of the minimum and maximal values of the mass fraction $c$ for matching densities (solid lines) and non-matching densities (dash-dotted lines). (a) Dissipation $\frac{dE}{dt}$. The value of the dissipation remains negative ensuring that Inequality (4.20) is satisfied at the discrete level. (b) Mass $\int_{\Omega} \rho c \, dx$ of fluid 1. Difference of initial mass is due to the random initial conditions for the two simulations. (c) Scalar variable $\xi$. The value remains close to 1 (error is of order $10^{-5}$), ensuring that the modified energy $\tau^{n+1}$ is close to the real energy $E^{n+1}$. (d) Minimum and maximum of $c$. The mass fraction lies in the physically relevant interval (0, 1) throughout the simulation.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Simulation of compressible Navier-Stokes-Cahn-Hilliard model nonmatching densities, exchange term ($F_c(\rho, c) \neq 0$), same friction effects for both fluids (top row) and contrast of friction forces (bottom row). (a) $t = 0$. (b) $t = 0.1$. (c) $t = 1$. (d) $t = 5$. (e) $t = 0$. (f) $t = 0.1$. (g) $t = 1$. (h) $t = 5$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Temporal evolution of the dissipation of the energy $\frac{dE}{dt}$, mass of the fluid 1 given by $\int_{\Omega} \rho c \, dx$, scalar variable $\xi$, and of the minimum and maximal values of the mass fraction $c$ for matching densities (solid lines) and non-matching densities (dash-dotted lines). (a) Dissipation $\frac{dE}{dt}$. (b) Total mass $\int_{\Omega} \rho \, dx$. (c) Scalar variable $\xi$. (d) min $c$ and max $c$.</figcaption></figure>

the variable $\xi$ increases slightly with time (compare Figs. 2c and 4c). A possible remedy to this issue is discussed in the conclusion of this article.
## 5.2. Two-dimensional numerical test cases

We now simulate the G-NSCH system (1.1)-(1.4) in two dimensions. Details about the two-dimensional numerical scheme can be found in Appendix C.

**Phase separation with non-matching densities and contrast of friction strengths.** We use $N_x = N_y = 64$ cells in each direction. We fix the final time at $T = 1$. We set up two simulations, both with no exchange terms $F_c(\rho, c) = 0$. For the first, we consider no contrast of friction effects, i.e. $\kappa_1 = \kappa_2 = 10$, $\nu_1 = \nu_2 = 0.01$, and $\eta_1 = \eta_2 = 0.02$. For the second simulation, we take $\kappa_1 = 0$, $\kappa_2 = 10$, $\nu_1 = \nu_2 = 0.01$, and $\eta_1 = \eta_2 = 0.02$. Hence, fluid 2 has stronger friction effects.

The other parameters are $\gamma = \frac{1}{800}$, $\theta = 4$, $\alpha_1 = 0.8$, $\alpha_2 = 1.2$, $\iota = 10^{-4}$, $C_0 = 100$, $a = 1.5$. The tolerance of GMRES solver is set to $tol = 10^{-10}$.

The initial velocities in both directions are constants in space $v_x^0 = 0.5$ and $v_y^0 = 0.5$. The initial density is also constant in space $\rho^0 = 0.8$. The initial mass fraction is set to a perturbed constant $c^0 = 0.3 - 0.05r$, with $r$ a random uniform number for each cell center.

Figure 5 compares the temporal evolution of the density of fluid 1 given by $\rho c$ for both cases. We observe that both solutions depict phase separation and progressive coarsening of small aggregates into larger ones. This phenomenon occurs as the fluid is transported to the top right corner (we recall that we implemented periodic boundary conditions). Careful inspection of the relative density $\rho c$ distribution inside each aggregates reveals the effect of the contrast of friction between the two solutions. Indeed, as the fluid 1 encounters a resistance when transported by the flow (because it pushes a fluid that experiences more friction), the mass of fluid 1 seems to concentrate in the top right corner of each aggregate. This can be observed inspecting the level lines depicted on Figure 6. Indeed on this figure, we see that the top of each aggregate is not localized in the center of the aggregates but is shifted to the top-right. This indicates that the 2D scheme captures correctly the effect seen with the 1D numerical scheme when contrast of friction between the two phases is considered.

## 5.3. Convergence tests

We study the numerical convergence of the one dimensional scheme (4.4)-(4.16) with $\kappa(\rho, c) = 0$ and $F_c(\rho, c) = 0$. The computational domain is $\Omega = (0, 1)$. The final time is $T = 0.05$. The other parameters $\gamma, \beta, \eta, \nu, a, \alpha_1, \alpha_2$ are chosen as for the previous 1D test case with non-matching densities (see Sect. 5.1). The initial condition for the mass fraction is given by

$$c^0 = 0.4 + 0.01 \cos(6\pi x).$$

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

We set the tolerance $rtol$ of the GMRES algorithm to $rtol = 10^{-6}$ for the spatial convergence test and $rtol = 10^{-10}$ for the temporal convergence test.

### 5.3.1. Convergence in space

We fix the time step to $\Delta t = 1 \times 10^{-5}$ and we vary the grid size. We choose an increasing number of cells $N_x = \{64, 128, 256, 512, 1024, 2048\}$. For each quantity $c, v, \rho$, we compute the discrete errors

$$ \begin{aligned} \text{error}(\rho_{\Delta x}, \rho_{\Delta x/2}) &= \|\rho_{\Delta x} - \rho_{\Delta x/2}\|_{L^{\infty}(0,T;L^a(\Omega))}, \\ \text{error}(c_{\Delta x}, c_{\Delta x/2}) &= \|c_{\Delta x} - c_{\Delta x/2}\|_{L^2(0,T;L^2(\Omega))}, \\ \text{error}(v_{\Delta x}, v_{\Delta x/2}) &= \|v_{\Delta x} - v_{\Delta x/2}\|_{L^2(0,T;L^2(\Omega))}, \end{aligned} \quad (5.2) $$

where $(\rho_{\Delta x/2}, c_{\Delta x/2}, v_{\Delta x/2})$ denotes the solution computed using twice the number of cells of the simulation that computes the solution $(\rho_{\Delta x}, c_{\Delta x}, v_{\Delta x})$.
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Two dimensional simulations of compressible Navier-Stokes-Cahn-Hilliard model with nonmatching densities, same friction effects for both fluids (top row) and contrast of friction (bottom row). (a) $t = 0$. (b) $t = 0.1$. (c) $t = 0.5$. (d) $t = 1$. (e) $t = 0$. (f) $t = 0.1$. (g) $t = 1$. (h) $t = 5$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Relative density of fluid 1 at time $t = 0.25$ considering a contrast of friction strengths between the two fluids. The black circles represent the level $\rho c = \max(\rho c) - 0.03$. This corresponds to the tops of each aggregate. (a) No contrast of friction strengths. (b) Contrast of friction strengths.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Convergence in space for the total density $\rho$, the mass fraction $c$ and the velocity $v$. The orange dashed line represents the slope 1. (a) Spatial convergence of density $\rho$. (b) Spatial convergence of mass fraction $c$. (c) Spatial convergence of velocity $v$.</figcaption></figure>

To compute the discrete norms, we save the solution every $\Delta t_{save} = 0.001$. The norms in (5.2) are computed following

$$ \begin{align*} \|\rho\Delta x - \rho\Delta x/2\|_{L^{\infty}(0,T;L^a(\Omega))} &= \max_{t_{save}} \left( \frac{\Delta x}{2} \sum_{j=1}^{N_x} (\rho\Delta x(x_j) - \rho\Delta x/2(x_j))^a \right)^{1/a}, \\ \|c\Delta x - c\Delta x/2\|_{L^2(0,T;L^2(\Omega))} &= \left( \sum_{t_{save}} \Delta t_{save} \left( \left( \frac{\Delta x}{2} \sum_{j=1}^{N_x} (c\Delta x(x_j) - c\Delta x/2(x_j))^2 \right)^{1/2} \right)^2 \right)^{1/2}, \\ \|v\Delta x - v\Delta x/2\|_{L^2(0,T;L^2(\Omega))} &= \left( \sum_{t_{save}} \Delta t_{save} \left( \left( \frac{\Delta x}{2} \sum_{j=1}^{N_x} (v\Delta x(x_j) - v\Delta x/2(x_j))^2 \right)^{1/2} \right)^2 \right)^{1/2}, \end{align*} $$

with $N_x$ the number of points on the $\Delta x/2$ grid, and $t_{save}$ the array of times at which snapshots of the solutions have been taken. Hence, the solution on the coarse grid $\Delta x$ is extended on the fine grid $\Delta x/2$ using the nearest solution from the coarse grid.

We arrive at the results given in Figure 7. As expected by the upwind scheme, the spatial order of convergence is a little less than 1 for the total density $\rho$ (see Fig. 7a) and the velocity $v$ (see Fig. 7c). We recover first order for the mass fraction.

### 5.3.2. Convergence in time

We here fix the grid size and select $N_x = 128$ points. We choose $\Delta t = 1 \times 10^{-4}$, and decrease the time steps according to $\Delta t_{array} = \{\Delta t, \frac{\Delta t}{2}, \frac{\Delta t}{4}, \frac{\Delta t}{8}, \frac{\Delta t}{16}, \frac{\Delta t}{32}, \frac{\Delta t}{64}\}$. We deactivate the time step adaptive strategy from the CFL condition. The other parameters and initial conditions are chosen as in the spatial convergence test (see Sect. 5.3.1). To check the convergence in time of our scheme, we compute the errors between two solutions computed with two time steps $\Delta t$ that differ only from a factor $\frac{1}{2}$. We denote these two different solutions by $(\rho\Delta t, c\Delta t, v\Delta t)$ and $(\rho\Delta t/2, c\Delta t/2, v\Delta t/2)$. We use the same method as for the spatial convergence computations,
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Convergence in time for the total density $\rho$, the mass fraction $c$ and the velocity $v$. The orange dashed line represents the slope 1. (a) Temporal convergence of density $\rho$. (b) Temporal convergence of mass fraction $c$. (c) Temporal convergence of velocity $v$.</figcaption></figure>

we save the solutions every $\Delta t_{\text{save}} = 0.001$. We compute the norms

$$ 
\begin{aligned}
\text{error}(\rho_{\Delta t}, \rho_{\Delta t/2}) &= \|\rho_{\Delta t} - \rho_{\Delta t/2}\|_{L^{\infty}(0,T;L^a(\Omega))}, \\
\text{error}(c_{\Delta t}, c_{\Delta t/2}) &= \|c_{\Delta t} - c_{\Delta t/2}\|_{L^2(0,T;L^2(\Omega))}, \\
\text{error}(\mathbf{v}_{\Delta t}, \mathbf{v}_{\Delta t/2}) &= \|\mathbf{v}_{\Delta t} - \mathbf{v}_{\Delta t/2}\|_{L^2(0,T;L^2(\Omega))}.
\end{aligned} 
\quad (5.3) 
$$

with

$$ 
\begin{aligned}
\|\rho_{\Delta t} - \rho_{\Delta t/2}\|_{L^{\infty}(0,T;L^a(\Omega))} &= \max_{t_{\text{save}}} \left( \Delta x \sum_{j=1}^{N_x} (\rho_{\Delta t}(x_j) - \rho_{\Delta t/2}(x_j))^a \right)^{1/a}, \\
\|c_{\Delta t} - c_{\Delta t/2}\|_{L^2(0,T;L^2(\Omega))} &= \left( \sum_{t_{\text{save}}} \Delta t_{\text{save}} \left( \left( \Delta x \sum_{j=1}^{N_x} (c_{\Delta t}(x_j) - c_{\Delta t/2}(x_j))^2 \right)^{1/2} \right)^2 \right)^{1/2}, \\
\|\mathbf{v}_{\Delta t} - \mathbf{v}_{\Delta t/2}\|_{L^2(0,T;L^2(\Omega))} &= \left( \sum_{t_{\text{save}}} \Delta t_{\text{save}} \left( \left( \Delta x \sum_{j=1}^{N_x} (\mathbf{v}_{\Delta t}(x_j) - \mathbf{v}_{\Delta t/2}(x_j))^2 \right)^{1/2} \right)^2 \right)^{1/2},
\end{aligned} 
$$

We obtain the results depicted in Figure 8. We observe that the order of convergence in time for our scheme is exactly 1 for the three quantities.

## 6. CONCLUSION AND PERSPECTIVES

We presented a generalized model of diphasic compressible fluid termed G-NSCH, that comprises possible mass transfer between the two phases and friction effects. Under simplifying assumptions, summarized in Section 2.1, we proved the existence of weak solutions of the G-NSCH system. We also proposed a numerical scheme and prove, under the same simplifying assumptions, that it is stable and structure preserving (i.e. it ensures the physically relevant bounds for the mass fraction $c$, and it satisfies an energy dissipation inequality). For the numerical simulations, we chose relevant functionals, thus, relaxing the simplifying assumptions that were necessary for the analysis. We presented numerical simulations showing that our numerical scheme possesses the robustness found analytically. The numerical simulations allowed us to show the ability of our
model to represent diphasic fluids with matching or non-matching densities for the two phases. Furthermore, we computed numerically the spatial and temporal convergence rates of our numerical scheme.

Our model and numerical scheme allow to use physically relevant choice of functional and to consider contrast of properties between the two phases of the fluid. Our aim is to perform efficient simulation of general compressible diphasic fluids while being able to capture instabilities that could emerge considering contrasts of properties such as Saffman-Taylor or Rayleigh-Taylor instabilities. However, we emphasize that to achieve this latter goal, we have to be able to capture accurately the possible fine structures appearing during the numerical simulations. We plan to improve our numerical scheme in several ways. First, we plan to adapt the Relaxed version of the Generalized SAV method. Indeed, we observed during our numerical simulation that the variable $\xi$ is between $10^{-5}$ and $10^{-3}$. However, with the transfer term $F_c(\rho, c) \neq 0$, as the simulation progresses, the gap between $r$ and the real energy increases, i.e. $\xi$ becomes larger. As shown in [66], this problem is solved with the relaxed G-SAV method and the use of this method is one of our further developments. In the same work [66], it is mentioned that this relaxed method works well even when an external force is comprised in the model. In our case, this external force will take the form of a mass transfer or mass source term. The second improvement concerns the accuracy of our scheme. As shown in our work, the temporal and spatial orders of our scheme do not exceed 1. Thus, we will aim to design a high-order finite element scheme for the generalized compressible NSCH system that will remain structure preserving taking advantage of the flexibility of the relaxed GSAV method. On another aspect, we plan to use the reduced version of the G-NSCH model presented in the Appendix B to represent tumor growth while removing non-necessary effects such as inertia. Our goal is to present a model and 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 [19]. Furthermore, analytical aspects of this work can also be improved. This direction is challenging because as pointed out in the present work, necessary tools to perform the solutions' existence proof do not work with physically or biologically relevant potentials or mobility functions. In fact, singular potentials, degenerate mobilities and degenerate viscosity functions are not allowed. One possible solution is to derive a Bresch-Desjardins entropy estimate [13, 14] for the compressible NSCH as it has been done recently by Vasseur and Yu [64], and Bresch et al. [15], for the compressible Navier-Stokes model with degenerate viscosities.

To conclude, we emphasize that the G-NSCH model is the basis of a reduced system that takes into account only biologically relevant physical effects that play a role in tumor evolution (presented in the present article in Appendix B as Problem 2). Therefore, this work has to be seen as the first part. In a subsequent work, relying heavily on the present one, we will focus on numerical simulations, and sensitivity analysis of the reduced model.

# APPENDIX A. DERIVATION OF THE MODEL

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

## A.1. Notation and definitions

We formulate our problem in Eulerian coordinates and in a smooth bounded domain $\Omega \subset \mathbb{R}^d$ (where $d = \{1, 2, 3\}$ is the dimension). The balance laws derived in the following sections are in local form.

We have two fluids in the model where $\rho_1, \rho_2$ are the relative densities of respectively fluid 1 and 2. Thus, $\rho_i$ represents the mass $M_i$ of the fluid per volume occupied by the $i$th phase $V_i$, i.e.

$$ \rho_i = \frac{M_i}{V_i}. $$

Then, we define the volume fractions $\varphi_1, \varphi_2$ which are defined by the volume occupied by the $i$th phase over the total volume of the mixture

$$ \varphi_i = \frac{V_i}{V}. $$
Therefore, the mass density of population i which is the mass of population i in volume V is given by

$$ \phi_i = \rho_i \phi_i. $$

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

$$ \varphi_1 + \varphi_2 = 1. $$

The total density of the mixture is then given by

$$ \rho = \phi_1 + \phi_2. $$

We also introduce the mass fractions $c_i = M_i/M$ and we have the relations

$$ \rho c_i = \phi_i, \quad \text{and} \quad c_1 = (1 - c_2). \qquad (A.1) $$

We denote by $p$ the pressure inside the mixture and $\mathbf{v}_1, \mathbf{v}_2$ are the velocities of the different phases. We use a mass-average mixture velocity

$$ \mathbf{v} = \frac{1}{\rho}(\phi_1 \mathbf{v}_1 + \phi_2 \mathbf{v}_2). \qquad (A.2) $$

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

$$ \frac{Dg}{Dt} = \frac{\partial g}{\partial t} + \mathbf{v} \cdot \nabla g, \qquad (A.3) $$

and indicate the definition of the differential operator

$$ \mathbf{v} \cdot \nabla g = \sum_{j=1}^{d} \mathbf{v}_j \frac{\partial g}{\partial x_j}. $$

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

## A.2. Mass balance equations

We have the mass balance equations

$$ \left\{ \begin{aligned} \frac{\partial \phi_1}{\partial t} + \text{div}(\phi_1 \mathbf{v}_1) &= F_1(\rho, c_1, c_2), \\ \frac{\partial \phi_2}{\partial t} + \text{div}(\phi_2 \mathbf{v}_2) &= F_2(\rho, c_1, c_2). \end{aligned} \right. \qquad (A.4) $$

The functions $F_i(\rho, c_1, c_2)$ ($i = 1, 2$) act as source or exchange terms of mass.

Summing the two equations, we obtain the continuity equation for the total density of the mixture, using $c = c_1$, and the relations (A.1), we obtain the balance equation for the total density of the mixture

$$ \frac{\partial \rho}{\partial t} + \text{div}(\rho \mathbf{v}) = F_1 + F_2 =: F_{\rho}. \qquad (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

$$ \frac{\partial \rho c}{\partial t} + \text{div}(\rho c \mathbf{v}_1) = F_1(\rho, c, 1-c) =: F_c. \qquad (A.6) $$

The mass of the component 1 is transported by the average velocity $\mathbf{v}$ and the remaining diffusive flux $\mathbf{J}_1 = \rho c(\mathbf{v} - \mathbf{v}_1)$. Therefore, we can replace the previous equation by

$$ \frac{\partial \rho c}{\partial t} + \text{div}(\rho c \mathbf{v}) = \text{div}(\mathbf{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

$$ \frac{\partial \rho c}{\partial t} + \text{div}(\rho c \mathbf{v}) = \rho \frac{Dc}{Dt} + c \left[ \frac{\partial \rho}{\partial t} + \text{div}(\rho \mathbf{v}) \right] = \rho \frac{Dc}{Dt} + c F_{\rho}. $$

Altogether, we obtain the balance equation for the mass fraction of the component 1

$$ \rho \frac{Dc}{Dt} = \text{div}(\mathbf{J}_1) + F_c - cF_{\rho}. \qquad (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 [17, 18]).

### A.3. Balance of linear momentum

We write the balance of linear momentum [23], which describes the evolution of the velocity $\mathbf{v}$ due to internal stresses and external forces. 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 [25]. Altogether, we assume that the stress tensor is a function of the total density $\rho$, the order parameter $c$ (i.e. the mass fraction of fluid 1), its gradient $\nabla c$, and the total velocity of the mixture $\mathbf{v}$, i.e.

$$ \boldsymbol{\sigma} = \boldsymbol{\sigma}(\rho, c, \nabla c, \mathbf{v}). $$

The friction around the pores of the medium is modeled by a drag term in the balance equation [56] with a permeability coefficient $\kappa(\rho, c) = \kappa_1(\rho, c) + \kappa_2(\rho, 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 [23]

$$ \frac{\partial \rho \mathbf{v}_j}{\partial t} + \text{div}(\rho \mathbf{v}_j \mathbf{v}) = \text{div}(\boldsymbol{\sigma})_j - \kappa(\rho, c) \mathbf{v}_j + F_{\mathbf{v}_j}, $$

where $F_{\mathbf{v}_j}(\mathbf{v}_j, \rho)$ 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

$$ \frac{\partial \rho \mathbf{v}_j}{\partial t} + \text{div}(\rho \mathbf{v}_j \mathbf{v}) = \rho \frac{D\mathbf{v}_j}{Dt} + \mathbf{v}_j \left[ \frac{\partial \rho}{\partial t} + \text{div}(\rho \mathbf{v}) \right] = \rho \frac{D\mathbf{v}_j}{Dt} + \mathbf{v}_j F_{\rho} + F_{\mathbf{v}_j}. $$

Therefore, we have

$$ \rho \frac{D\mathbf{v}_j}{Dt} = \text{div}(\boldsymbol{\sigma})_j - (\kappa(\rho, c) + F_{\rho})\mathbf{v}_j + F_{\mathbf{v}_j}. $$

We can rewrite the balance of linear momentum in a more compact form

$$ \rho \frac{D\mathbf{v}}{Dt} = \text{div}(\boldsymbol{\sigma}) - (\kappa(\rho, c) + F_{\rho})\mathbf{v} + F_{\mathbf{v}}, \qquad (A.8) $$

where $F_{\mathbf{v}}(\mathbf{v}, \rho)$ is the vector of coordinates $F_{\mathbf{v}_j}$.

### A.4. Energy balance

The total energy of the mixture is the sum of the kinetic energy $\rho \frac{1}{2}|\mathbf{v}|^2$ and of the internal energy $\rho u$, where $u = u(\rho, c, \nabla c)$ is a specific internal energy. Compared to the classical conservation law for the total energy, we have an additional energy flux $\tau \frac{Dc}{Dt}$. Indeed, due to the interface region, surface effects must be taken into account. Following this direction, Gurtin [34] 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 $\tau$.

Since we assume that the system is maintained in an isothermal state, the balance equation for the energy is given by Eck et al. [23]

$$ \begin{aligned} & \frac{\partial}{\partial t} \left( \rho \frac{1}{2} |\mathbf{v}|^2 + \rho u \right) + \operatorname{div} \left( \rho \left( \frac{1}{2} |\mathbf{v}|^2 + u \right) \mathbf{v} \right) \\ & = \operatorname{div}(\boldsymbol{\sigma}^T \mathbf{v}) + \operatorname{div} \left( \tau \frac{\mathrm{D}c}{\mathrm{D}t} \right) - \operatorname{div}(\mathbf{q}) + \rho g + c_{\rho} F_{\rho} + c_{c} F_{c} + c_{\mathbf{v}} F_{\mathbf{v}}, \end{aligned} \quad (A.9) $$

where $\mathbf{q}$ is the heat flux and $\rho 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. [36, 45]). The prefactors $c_{\rho}, c_{c}, c_{\mathbf{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 have

$$ \frac{\partial}{\partial t} \left( \rho \frac{1}{2} |\mathbf{v}|^2 + \rho u \right) + \operatorname{div} \left( \rho \left( \frac{1}{2} |\mathbf{v}|^2 + u \right) \mathbf{v} \right) = \rho \left[ \frac{\mathrm{D}}{\mathrm{D}t} \left( \frac{1}{2} |\mathbf{v}|^2 + u \right) \right] + \left( \frac{1}{2} |\mathbf{v}|^2 + u \right) F_{\rho}. $$

Applying the chain rule to the kinetic part, we obtain

$$ \rho \frac{\mathrm{D}}{\mathrm{D}t} \left( \frac{1}{2} |\mathbf{v}|^2 \right) = \rho \mathbf{v} \cdot \frac{\mathrm{D}\mathbf{v}}{\mathrm{D}t}, $$

and using the balance of linear momentum (A.8), we arrive to

$$ \rho \mathbf{v} \cdot \frac{\mathrm{D}\mathbf{v}}{\mathrm{D}t} = \mathbf{v} \cdot \operatorname{div}(\boldsymbol{\sigma}) - (\kappa(\rho, c) + F_{\rho}) |\mathbf{v}|^2 + F_{\mathbf{v}} \cdot \mathbf{v}. $$

Using these previous equations inside (A.9), we obtain the balance equation for the internal energy

$$ \begin{aligned} \rho \frac{\mathrm{D}u}{\mathrm{D}t} = \operatorname{div}(\boldsymbol{\sigma}^T \mathbf{v}) - \mathbf{v} \cdot \operatorname{div}(\boldsymbol{\sigma}) + \operatorname{div} \left( \tau \frac{\mathrm{D}c}{\mathrm{D}t} \right) &+ (\kappa(\rho, c) + F_{\rho}) |\mathbf{v}|^2 - F_{\mathbf{v}} \mathbf{v} \\ &- \operatorname{div}(\mathbf{q}) + \rho g - \left( \frac{1}{2} |\mathbf{v}|^2 + u \right) F_{\rho} + c_{\rho} F_{\rho} + c_{c} F_{c} + c_{\mathbf{v}} F_{\mathbf{v}}. \end{aligned} $$

However, since

$$ \mathbf{v} \cdot (\operatorname{div}(\boldsymbol{\sigma})) - \operatorname{div}(\boldsymbol{\sigma}^T \mathbf{v}) = -\boldsymbol{\sigma} : \nabla \mathbf{v}, $$

where $\nabla \mathbf{v} = (\partial_{x_j} \mathbf{v}_i)_{i,j=1,...,d}$ is the Jacobi matrix and, we have $A:B = \sum_{i,j} A_{ij} B_{ij}$, for two matrices $A, B$. Altogether, we have the balance equation for the internal energy

$$ \begin{aligned} \rho \frac{\mathrm{D}u}{\mathrm{D}t} = \boldsymbol{\sigma} : \nabla \mathbf{v} + \operatorname{div} \left( \tau \frac{\mathrm{D}c}{\mathrm{D}t} \right) &+ (\kappa(\rho, c) + F_{\rho}) |\mathbf{v}|^2 - F_{\mathbf{v}} \mathbf{v} \\ &- \operatorname{div}(\mathbf{q}) + \rho g - \left( \frac{1}{2} |\mathbf{v}|^2 + u \right) F_{\rho} + c_{\rho} F_{\rho} + c_{c} F_{c} + c_{\mathbf{v}} F_{\mathbf{v}}. \end{aligned} \quad (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(\rho, c, \nabla c)$ and the Helmholtz free energy $\mathcal{F} = \mathcal{F}(\rho, c, \nabla c)$, both related through the equation

$$ \mathcal{F} = u - Ts, \quad (A.11) $$
where $T$ denotes the temperature.

From the mass balance equation (A.5), we have the entropy balance equation

$$ \frac{\partial \rho s}{\partial t} + \mathrm{div}(s\rho \mathbf{v}) = \rho \frac{\mathrm{D}s}{\mathrm{D}t} + s \left[ \frac{\partial \rho}{\partial t} + \mathrm{div}(\rho \mathbf{v}) \right] = \rho \frac{\mathrm{D}s}{\mathrm{D}t} + sF_{\rho}. \quad (A.12) $$

Then, using the definition of the Helmholtz free energy (A.11) and the balance of energy (A.10), we obtain

$$ \begin{aligned} \rho \frac{\mathrm{D}s}{\mathrm{D}t} &= -\frac{\rho}{T} \frac{\mathrm{D}\mathcal{F}}{\mathrm{D}t} + \frac{\rho}{T} \frac{\mathrm{D}u}{\mathrm{D}t} \\ &= -\frac{\rho}{T} \frac{\mathrm{D}\mathcal{F}}{\mathrm{D}t} + \frac{1}{T} \left[ \boldsymbol{\sigma} : \nabla \mathbf{v} + \mathrm{div} \left( \boldsymbol{\tau} \frac{\mathrm{D}c}{\mathrm{D}t} \right) + (\kappa(\rho, c) + F_{\rho}) |\mathbf{v}|^2 - F_v \mathbf{v} \right. \\ & \qquad \left. - \mathrm{div}(\mathbf{q}) + \rho g - \left( \frac{1}{2} |\mathbf{v}|^2 + u \right) F_{\rho} + c_{\rho} F_{\rho} + c_c F_c + c_{\mathbf{v}} F_{\mathbf{v}} \right], \end{aligned} \quad (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) [23]. 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 [50] and Müller [55] 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 [55], we state the second law as an entropy inequality, i.e., the Clausius-Duhem inequality in the local form [23]

$$ \rho \frac{\mathrm{D}s}{\mathrm{D}t} \geq -\mathrm{div}\left(\frac{\mathbf{q}}{T}\right) + \frac{\rho g}{T} + \mathrm{div}(\mathcal{J}), \quad (A.14) $$

where $\mathcal{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

$$ \frac{\rho}{T} \frac{\mathrm{D}\mathcal{F}}{\mathrm{D}t} - \frac{1}{T} \left[ \boldsymbol{\sigma} : \nabla \mathbf{v} + \mathrm{div} \left( \boldsymbol{\tau} \frac{\mathrm{D}c}{\mathrm{D}t} \right) + (\kappa(\rho, c) + F_{\rho}) |\mathbf{v}|^2 - F_v \mathbf{v} - \left( \frac{1}{2} |\mathbf{v}|^2 + u \right) F_{\rho} + c_{\rho} F_{\rho} + c_c F_c + c_{\mathbf{v}} F_{\mathbf{v}} \right] + \mathrm{div}(\mathcal{J}) \leq 0. \quad (A.15) $$

Then, using the chain rule

$$ \frac{\mathrm{D}\mathcal{F}}{\mathrm{D}t} = \frac{D\rho}{Dt} \frac{\partial\mathcal{F}}{\partial\rho} + \frac{Dc}{Dt} \frac{\partial\mathcal{F}}{\partial c} + \frac{D\nabla c}{Dt} \cdot \frac{\partial\mathcal{F}}{\partial\nabla c}, $$

and

$$ \frac{D\nabla c}{Dt} = \nabla \left[ \frac{Dc}{Dt} \right] - (\nabla \mathbf{v})^T \nabla c, \quad \frac{D\rho}{Dt} = -\rho \mathrm{div}(\mathbf{v}) + F_{\rho}, $$

in the entropy inequality (A.15), we obtain

$$ \begin{aligned} & \rho \left[ \left( -\rho \mathrm{div}(\mathbf{v}) + F_{\rho} \right) \frac{\partial \mathcal{F}}{\partial \rho} + \frac{Dc}{Dt} \frac{\partial \mathcal{F}}{\partial c} + \left( \nabla \left[ \frac{Dc}{Dt} \right] - (\nabla \mathbf{v})^T \nabla c \right) \cdot \frac{\partial \mathcal{F}}{\partial \nabla c} \right] - \mathrm{div} \left( \boldsymbol{\tau} \frac{Dc}{Dt} \right) - \boldsymbol{\sigma} : \nabla \mathbf{v} \\ & - \left[ (\kappa(\rho, c) + F_{\rho}) |\mathbf{v}|^2 - F_v \mathbf{v} - \left( \frac{1}{2} |\mathbf{v}|^2 + u \right) F_{\rho} + c_{\rho} F_{\rho} + c_c F_c + c_{\mathbf{v}} F_{\mathbf{v}} \right] + T \mathrm{div}(\mathcal{J}) \le 0. \end{aligned} \quad (A.16) $$
By the chain rule, we have

$$ \mathrm{div}\left(\boldsymbol{\tau} \frac{\mathrm{D}c}{\mathrm{D}t}\right) = \boldsymbol{\tau} \nabla\left[\frac{\mathrm{D}c}{\mathrm{D}t}\right] + \frac{\mathrm{D}c}{\mathrm{D}t} \mathrm{div}(\boldsymbol{\tau}). $$

Furthermore, we know that

$$ -\rho^2 \mathrm{div}(\mathbf{v}) \frac{\partial \mathcal{F}}{\partial \rho} = -\rho^2 \frac{\partial \mathcal{F}}{\partial \rho} \mathbb{I} : \nabla \mathbf{v}, $$

and

$$ -\rho ((\nabla \mathbf{v})^T \nabla c) \cdot \frac{\partial \mathcal{F}}{\partial \nabla c} = -\rho (\nabla c \otimes \frac{\partial \mathcal{F}}{\partial \nabla c}) : \nabla \mathbf{v}. $$

Gathering the previous three relations and reorganizing the terms of (A.16), we obtain

$$ \begin{aligned} & \left(-\rho^2 \frac{\partial \mathcal{F}}{\partial \rho} \mathbb{I} - \rho \nabla c \otimes \frac{\partial \mathcal{F}}{\partial \nabla c} - \boldsymbol{\sigma}\right) : \nabla \mathbf{v} + \left(\rho \frac{\partial \mathcal{F}}{\partial c} - \mathrm{div}(\boldsymbol{\tau})\right) \frac{\mathrm{D}c}{\mathrm{D}t} \\ & \qquad + \left(\rho \frac{\partial \mathcal{F}}{\partial \nabla c} - \boldsymbol{\tau}\right) \nabla \left[\frac{\mathrm{D}c}{\mathrm{D}t}\right] + T \mathrm{div}(\mathcal{J}) \\ & \qquad - \left[\left(\kappa(\rho, c) + F_{\rho}\right) |\mathbf{v}|^2 - F_{v} \mathbf{v} - \left(\frac{1}{2} |\mathbf{v}|^2 + u - \rho \frac{\partial \mathcal{F}}{\partial \rho}\right) F_{\rho} + c_{\rho} F_{\rho} + c_{c} F_{c} + c_{v} F_{v}\right] \le 0. \end{aligned} \quad (A.17) $$

Then, we use Liu's Lagrange multipliers method [50]. 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 $(\rho, c, \nabla\rho, \nabla c, \mathbf{v}, p)$

$$ \begin{aligned} -D_{\text{iss}} := & \left(-\rho^2 \frac{\partial \mathcal{F}}{\partial \rho} \mathbb{I} - \rho \nabla c \otimes \frac{\partial \mathcal{F}}{\partial \nabla c} - \boldsymbol{\sigma}\right) : \nabla \mathbf{v} \\ & + \left(\rho \frac{\partial \mathcal{F}}{\partial c} - \mathrm{div}(\boldsymbol{\tau}) - \rho L_c\right) \frac{\mathrm{D}c}{\mathrm{D}t} + \left(\rho \frac{\partial \mathcal{F}}{\partial \nabla c} - \boldsymbol{\tau}\right) \nabla \left[\frac{\mathrm{D}c}{\mathrm{D}t}\right] + T \mathrm{div}(\mathcal{J}) \\ & - \left[\left(\kappa(\rho, c) + F_{\rho}\right) |\mathbf{v}|^2 - F_{v} \mathbf{v} - \left(\frac{1}{2} |\mathbf{v}|^2 + u - \rho \frac{\partial \mathcal{F}}{\partial \rho}\right) F_{\rho} + c_{\rho} F_{\rho} + c_{c} F_{c} + c_{v} F_{v}\right] \\ & - L_c \left(\rho \frac{\mathrm{D}c}{\mathrm{D}t} - \mathrm{div}(\mathbf{J}_1) - F_c - c F_{\rho}\right) \le 0. \end{aligned} \quad (A.18) $$

Since,

$$ \mathrm{div}(L_c \mathbf{J}_1) = L_c \mathrm{div}(\mathbf{J}_1) + \nabla L_c \cdot \mathbf{J}_1, $$

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

$$ \begin{aligned} -D_{\text{iss}} := & \left(-\rho^2 \frac{\partial \mathcal{F}}{\partial \rho} \mathbb{I} - \rho \nabla c \otimes \frac{\partial \mathcal{F}}{\partial \nabla c} - \boldsymbol{\sigma}\right) : \nabla \mathbf{v} \\ & + \left(\rho \frac{\partial \mathcal{F}}{\partial c} - \mathrm{div}(\boldsymbol{\tau}) - \rho L_c\right) \frac{\mathrm{D}c}{\mathrm{D}t} + \left(\rho \frac{\partial \mathcal{F}}{\partial \nabla c} - \boldsymbol{\tau}\right) \nabla \left[\frac{\mathrm{D}c}{\mathrm{D}t}\right] + \mathrm{div}(T\mathcal{J} + L_c \mathbf{J}_1) - \nabla L_c \cdot \mathbf{J}_1 \\ & - \left[\left(\kappa(\rho, c) + F_{\rho}\right) |\mathbf{v}|^2 - F_{v} \mathbf{v} - \left(\frac{1}{2} |\mathbf{v}|^2 + u - \rho \frac{\partial \mathcal{F}}{\partial \rho}\right) F_{\rho} + c_{\rho} F_{\rho} + c_{c} F_{c} + c_{v} F_{v} - L_c(F_c + cF_{\rho})\right] \le 0. \end{aligned} \quad (A.19) $$

## A.6. Constitutive assumptions and model equations

First of all, we assume that the free energy density $\mathcal{F}$ is of Ginzburg-Landau type and has the following form [17, 18]

$$ \mathcal{F}(\rho, c, \nabla c) := \psi_0(\rho, c) + \frac{\gamma}{2} |\nabla c|^2, \quad (A.20) $$
where $\psi_0$ is the homogeneous free energy accounting for the processes of phase separation and the gradient term $\frac{\gamma}{2}|\nabla 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. [65], 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.19), we first choose

$$ \boldsymbol{\tau} := \rho \frac{\partial \mathcal{F}}{\partial \nabla c} = \gamma \rho \nabla c. $$

Then, we define the chemical potential $\mu(\rho, c, \nabla c)$ by

$$ \mu := \frac{\partial \mathcal{F}}{\partial c} - \frac{1}{\rho} \mathrm{div}(\boldsymbol{\tau}) = \frac{\partial \mathcal{F}}{\partial c} - \frac{1}{\rho} \mathrm{div} \left( \rho \frac{\partial \mathcal{F}}{\partial \nabla c} \right) = \frac{\partial \psi_0}{\partial c} - \frac{\gamma}{\rho} \mathrm{div}(\rho \nabla c), $$

which in turn gives a condition for the Lagrange multiplier

$$ L_c = \mu. \qquad (A.21) $$

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

$$ \left( \rho \frac{\partial \mathcal{F}}{\partial c} - \mathrm{div}(\boldsymbol{\tau}) - \rho L_c \right) \frac{Dc}{Dt} + \left( \rho \frac{\partial \mathcal{F}}{\partial \nabla c} - \boldsymbol{\tau} \right) \nabla \left[ \frac{Dc}{Dt} \right] = 0. $$

Then, using classical results on isothermal diffusion [23, 51], we have

$$ \mathcal{J} := -\frac{\mu \mathbf{J}_1}{T}, \qquad (A.22) $$

and, using a generalized Fick's law, we have

$$ \mathbf{J}_1 := b(c) \nabla \mu, \qquad (A.23) $$

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

$$ \mathrm{div}(T\mathcal{J} + L_c\mathbf{J}_1) - \nabla L_c \cdot \mathbf{J}_1 = -b(c)|\nabla \mu|^2 \le 0. $$

Following [3, 51], we define the pressure inside the mixture

$$ p := \rho^2 \frac{\partial \psi_0}{\partial \rho}. \qquad (A.24) $$

From standard rheology, we assume that the fluid satisfies Newton's rheological laws. The stress tensor is composed of two parts for the viscous $\tilde{\mathbf{P}}$ and non-viscous $\mathbf{P}$ contributions of stress

$$ \boldsymbol{\sigma} := \mathbf{P} + \tilde{\mathbf{P}}, \qquad (A.25) $$

and we have by standard continuum mechanics (see e.g. [3, 10, 23])

$$ \left\{ \begin{array}{l} \mathbf{P} = -(p - \frac{\gamma}{2}|\nabla c|^2)\mathbb{I} - \gamma \rho \nabla c \otimes \nabla c, \\ \tilde{\mathbf{P}} = \nu(c)(\nabla \mathbf{v} + \nabla \mathbf{v}^T) + \lambda(c)\mathrm{div}(\mathbf{v})\mathbb{I}. \end{array} \right. \qquad (A.26) $$
In (A.26), $\nu(c)$ denotes the shear viscosity and $\lambda(c) = \eta(c) - \frac{2}{3}\nu(c)$ where $\eta(c)$ is the dilatational viscosity that encodes the response of the fluid to volume changes. The second term in the non-viscous part of the stress (namely $-\gamma(\rho\nabla c \otimes \nabla c)$) represents capillary stresses that act at the interface of the two populations.

Using (A.26), we can cancel terms in (A.19)

$$ \left( -\rho^2 \frac{\partial \mathcal{F}}{\partial \rho} \mathbb{I} - \rho \nabla c \otimes \frac{\partial \mathcal{F}}{\partial \nabla c} - \boldsymbol{\sigma} \right) : \nabla \mathbf{v} = 0. $$

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 $c_\rho$, such that

$$ - \left[ (\kappa(\rho, c) + F_{\rho}) |\mathbf{v}|^2 - F_{\nu} \mathbf{v} - \left( \frac{1}{2} |\mathbf{v}|^2 + u - \rho \frac{\partial \mathcal{F}}{\partial \rho} \right) F_{\rho} + c_{\rho} F_{\rho} + c_c F_c + c_{\nu} F_{\nu} - L_c(F_c + cF_{\rho}) \right] \le 0. $$

Reorganizing the terms we have

$$ -\kappa(\rho, c) |\mathbf{v}|^2 - F_{\rho} \left[ c_{\rho} + |\mathbf{v}|^2 - \left( \frac{1}{2} |\mathbf{v}|^2 + u - \rho \frac{\partial \mathcal{F}}{\partial \rho} \right) - \mu c \right] - F_{\nu} [c_{\nu} - \mathbf{v}] - F_c [c_c - \mu] \le 0. $$

The obvious choices are

$$ \begin{cases} c_{\rho} = -|\mathbf{v}|^2 + \left( \frac{1}{2} |\mathbf{v}|^2 + u - \rho \frac{\partial \mathcal{F}}{\partial \rho} \right) + \mu c, \\ c_{\nu} = \mathbf{v}, \\ c_c = \mu. \end{cases} $$

From the previous constitutive relations, we satisfy the dissipation inequality (A.19).

## A.7. Summary of the model's equations

Using the previous constitutive relations our general model is the following compressible Navier-Stokes-Cahn-Hilliard system

$$ \begin{align} \frac{\partial \rho}{\partial t} &= -\operatorname{div}(\rho \mathbf{v}) + F_{\rho}, \nonumber \\ \rho \frac{D c}{D t} &= \operatorname{div}(b(c) \nabla \mu) + F_c - c F_{\rho}, \nonumber \\ \rho \mu &= -\gamma \operatorname{div}(\rho \nabla c) + \rho \frac{\partial \psi_0}{\partial c}, \tag{A.27} \\ \rho \frac{D \mathbf{v}}{D t} &= -[\nabla p + \gamma \operatorname{div}(\rho \nabla c \otimes \nabla c)] + \operatorname{div}(\nu(c)(\nabla \mathbf{v} + \nabla \mathbf{v}^T)) \nonumber \\ &\quad - \frac{2}{3} \operatorname{div}(\nu(c) \operatorname{div}(\mathbf{v}) \mathbb{I}) + \operatorname{div}(\eta(c) \operatorname{div}(\mathbf{v}) \mathbb{I}) - (\kappa(\rho, c) + F_{\rho}) \mathbf{v} + F_{\nu}, \nonumber \end{align} $$

with $p$ defined in (A.24).

# APPENDIX B. MODEL REDUCTIONS, GENERAL ASSUMPTIONS AND BIOLOGICALLY RELEVANT CHOICES OF THE MODEL'S FUNCTIONS

## B.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 have

$$ F_c = -F_{1-c}, \tag{B.1} $$
leading to mass conservation

$$F_{\rho} = 0. \qquad (B.2)$$

Furthermore, we assume no external source of velocity and energy, leading to

$$F_{\mathbf{v}} = 0, \text{ and } F_{u} = 0. \qquad (B.3)$$

Furthermore, using the same simplifying assumption as in Abels and Feireisl [3] to avoid vacuum zones, our final reduced system of equations is

$$\frac{\partial \rho}{\partial t} + \text{div}(\rho \mathbf{v}) = 0, \qquad (B.4)$$

$$\frac{\partial \rho c}{\partial t} + \text{div}(\rho c \mathbf{v}) = \text{div}(b(c)\nabla\mu) + F_c, \qquad (B.5)$$

$$\rho\mu = -\gamma\Delta c + \rho\frac{\partial\psi_0}{\partial c}, \qquad (B.6)$$

$$\frac{\partial \rho \mathbf{v}}{\partial t} + \text{div}(\rho \mathbf{v} \otimes \mathbf{v}) = -\left[ \nabla p + \gamma \text{div} \left( \nabla c \otimes \nabla c - \frac{1}{2} |\nabla c|^2 \mathbf{I} \right) \right] + \text{div}(\nu(c)(\nabla \mathbf{v} + \nabla \mathbf{v}^T)) \\ - \frac{2}{3} \text{div}(\nu(c)\text{div}(\mathbf{v})\mathbf{I}) + \text{div}(\eta(c)\text{div}(\mathbf{v})\mathbf{I}) - \kappa(\rho, c)\mathbf{v}, \qquad (B.7)$$

**Problem 2: biologically relevant variant of the system.** We here refer to the two phases of the mixture as cell populations and not fluids. 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

$$\nabla p + \kappa(\rho, c)\mathbf{v} = -\gamma\text{div}\left(\nabla c \otimes \nabla c - \frac{1}{2}|\nabla c|^2\mathbf{I}\right) - F_{\rho}\mathbf{v}.$$

Assuming that one cell population proliferates while the other does not leads to

$$F_c = F_{\rho} = \rho c P_c(p), \quad \text{and} \quad F_{1-c} = 0,$$

with a pressure-dependent proliferation rate $P_c(p) \ge 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 that

$$P_c'(p) \le 0, \quad \text{and} \quad P_c(p) = 0 \quad \text{for} \quad p > p_{\text{max}}. \qquad (B.8)$$

Combining these changes, the model becomes

$$\left\{ \begin{array}{l} \frac{\partial \rho}{\partial t} + \text{div}(\rho \mathbf{v}) = \rho c P_c(p), \\ \frac{\partial \rho c}{\partial t} + \text{div}(\rho c \mathbf{v}) = \text{div}(b(c)\nabla\mu) + \rho c P_c(p), \\ \rho\mu = -\gamma\Delta c + \rho\frac{\partial\psi_0}{\partial c}, \\ \nabla p + \kappa(\rho, c)\mathbf{v} = -\gamma\text{div}\left(\nabla c \otimes \nabla c - \frac{1}{2}|\nabla c|^2\mathbf{I}\right) - \rho c P_c(p)\mathbf{v}. \end{array} \right. \qquad (B.9)$$

## B.2. Biologically consistent choices of functions

As said in the derivation of the model, the free energy density $\mathcal{F}$ is the sum of two terms: $\frac{\gamma}{2}|\nabla c|^2$ taking into account the surface tension effects existing between the phases of the mixture and the potential $\psi_0(\rho, c)$ representing the cell-cell interactions and pressure. Thus, we choose

$$\psi_0(\rho, c) = \psi_e(\rho) + \psi_{\text{mix}}(\rho, c), \qquad (B.10)$$
with $\psi_{\text{mix}}(\rho, c) = H(c) \log \rho + Q(c)$. Then, using the constitutive relation for the pressure we have

$$p(\rho, c) = \rho^2 \frac{\partial \psi_0}{\partial \rho} = p_e(\rho) + \rho H(c). \qquad (B.11)$$

The function $b(c)$ is the active mobility of the cells.

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)^{\alpha}, \quad \alpha \ge 1, \qquad (B.12)$$

where $C_b$ is a positive constant.

We use for the pressure a power law such that

$$p_e(\rho) = \frac{1}{a-1} \rho^{a-1}. \qquad (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. Fig. B.1a). 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. An example of double-well potential is given by

$$\psi_{\text{mix}} = \frac{1}{2}(\alpha_1(1-c)\log(\rho(1-c)) + \alpha_2 c \log(\rho c)) - \frac{\theta}{2}\left(c - \frac{1}{2}\right)^2 + k, \qquad (B.14)$$

thus giving

$$H(c) = \frac{1}{2}(\alpha_1(1-c) + \alpha_2 c), \quad Q(c) = \frac{1}{2}(\alpha_1(1-c)\log(1-c) + \alpha_2 c \log(c)) - \frac{\theta}{2}\left(c - \frac{1}{2}\right)^2 + k,$$

where $\theta > 1$, and $k, \alpha_1, \alpha_2 > 0$ are an arbitrary constants.

To meet the phenomenological observations of the interaction between cells when the mixture is composed of only one cell population, a single-well potential seems more appropriate [16, 20].

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 $\partial_c \psi_0 = 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 B.1b. A typical example of single-well potential which has been used for the modeling of living tissue and cancer [7, 20] is

$$\psi_{\text{mix}}(\rho, c) = -(1 - c_e) \log(\rho(1 - c)) - \frac{c^3}{3} - (1 - c_e) \frac{c^2}{2} - (1 - c_e)c + k, \qquad (B.15)$$

thus giving

$$H(c) = -(1 - c_e), \quad Q(c) = -(1 - c_e) \log(1 - c) - \frac{c^3}{3} - (1 - c_e) \frac{c^2}{2} - (1 - c_e)c + k, \qquad (B.16)$$

where $k$ is an arbitrary constant.
<figure><img src="image_10.png" /><figcaption>FIGURE B.1. For both figures $\rho = 1$. Double-well logarithmic potential (**left**) with $\alpha_1 = 1.2$ and $\alpha_2 = 0.8$ and single-well logarithmic potential (**right**). (a) Double-well potential. (b) Single-well potential.</figcaption></figure>

## APPENDIX C. DESCRIPTION OF THE TWO-DIMENSIONAL NUMERICAL SCHEME

We describe the two-dimensional scheme. This scheme possesses the same properties as the one-dimensional scheme.

We write the velocity field $\mathbf{v} = (u_x, u_y)$. System (1.1)–(1.4) with the transformation proposed at the beginning of this section, reads

$$ \partial_t \rho + \partial_x (\rho u_x) + \partial_y (\rho u_y) = 0, \qquad (C.1) $$

$$ \begin{aligned} & \partial_t \left( \rho \begin{bmatrix} u_x \\ u_y \end{bmatrix} \right) + \begin{bmatrix} \partial_x (\rho u_x^2 + p) \\ \partial_y (\rho u_y^2 + p) \end{bmatrix} + \begin{bmatrix} \partial_y (\rho u_x u_y) \\ \partial_x (\rho u_x u_y) \end{bmatrix} = 2 \begin{bmatrix} \partial_x (\nu(c) \partial_x u_x) \\ \partial_y (\nu(c) \partial_y u_y) \end{bmatrix} + \begin{bmatrix} \partial_y (\nu(c)(\partial_y u_x + \partial_x u_y)) \\ \partial_x (\nu(c)(\partial_y u_x + \partial_x u_y)) \end{bmatrix} \\ & \qquad + \begin{bmatrix} \partial_x ((\eta(c) - \frac{2}{3}\nu(c))(\partial_x u_x + \partial_y u_y)) \\ \partial_y ((\eta(c) - \frac{2}{3}\nu(c))(\partial_x u_x + \partial_y u_y)) \end{bmatrix} - \frac{\gamma}{2} \begin{bmatrix} \partial_y ((\partial_x c)^2 - (\partial_y c)^2) \\ \partial_y ((\partial_y c)^2 - (\partial_x c)^2) \end{bmatrix} - \gamma \begin{bmatrix} \partial_y (\partial_x c \partial_y c) \\ \partial_x (\partial_x c \partial_y c) \end{bmatrix} \\ & \qquad - \kappa(\rho, c) \begin{bmatrix} u_x \\ u_y \end{bmatrix}, \end{aligned} \qquad (C.2) $$

$$ \rho(\partial_t v + u_x \partial_x v + u_y \partial_y v) = \frac{1}{T'(v)} (\partial_x(b(c)\partial_x\mu) + \partial_y(b(c)\partial_y\mu)) + \frac{1}{T'(v)} F_c, \qquad (C.3) $$

$$ \rho\mu = -\gamma T'(v)(\partial_{xx}c + \partial_{yy}c) - \gamma T''(v)((\partial_x v)^2 + (\partial_y v)^2) + \rho \frac{\partial\psi_0}{\partial c}, \qquad (C.4) $$

$$ \frac{dr}{dt} = - \frac{r(t)}{E[t] + C_0} \int_{\Omega} b(c) |\nabla \mu|^2 - \mu F_c \, dx. \qquad (C.5) $$
We introduce the notations $U = (\rho, \rho u_x, \rho u_y)$, $G(U) = (0, -\kappa u_x, -\kappa u_y)$ and

$$
\begin{align*}
F(U) &= (\rho u_x, \rho u_x^2 + p - 2\nu(c)\partial_x u_x + \left(\frac{2}{3}\nu(c) - \eta(c)\right)(\partial_x u_x + \partial_y u_y) + \frac{1}{2}\gamma((\partial_x c)^2 - (\partial_y c)^2), \\
&\qquad \rho u_x u_y - \nu(c)(\partial_y u_x + \partial_x u_y) + \gamma\partial_x c\partial_y c), \\
K(U) &= (\rho u_y, \rho u_x u_y - \nu(c)(\partial_y u_x + \partial_x u_y) + \gamma\partial_x c\partial_y c), \\
&\qquad \rho u_y^2 + p - 2\nu(c)\partial_y u_y + \left(\frac{2}{3}\nu(c) - \eta(c)\right)(\partial_x u_x + \partial_y u_y) + \frac{1}{2}\gamma((\partial_y c)^2 - (\partial_x c)^2).
\end{align*}
$$

The stabilization (see [37,44]) of the Navier-Stokes part of our system reads, with $\iota > 0$ the relaxation parameter,

$$
\left\{
\begin{array}{l}
\partial_t U + \partial_x V + \partial_y W = G(U), \\
\partial_t V + A \partial_x U = -\frac{1}{\iota} (V - F(U)), \\
\partial_t W + B \partial_y U = -\frac{1}{\iota} (W - K(U)),
\end{array}
\right.
\qquad (C.6)
$$

in which $A = \operatorname{diag}(a_1, a_2, a_3)$ and $B = \operatorname{diag}(b_1, b_2, b_3)$. In the following, we choose

$$
\begin{align*}
a_1 &= a_2 = a_3 = \max \left\{ \sup \left( u_x + \sqrt{\partial_\rho p} \right)^2, \sup u_x^2, \sup \left( u_x - \sqrt{\partial_\rho p} \right)^2 \right\}, \\
b_1 &= b_2 = b_3 = \max \left\{ \sup \left( u_y + \sqrt{\partial_\rho p} \right)^2, \sup u_y^2, \sup \left( u_y - \sqrt{\partial_\rho p} \right)^2 \right\}.
\end{align*}
$$

We assume that our two-dimensional domain is a square $[0, L] \times [0, L]$. We discretize the domain using square control volumes of size $\Delta x \times \Delta 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.

$$
\rho_{j,i} = \frac{1}{\Delta x \Delta y} \int_{x_{j-1/2}}^{x_{j+1/2}} \int_{y_{j-1/2}}^{y_{j+1/2}} \rho(\mathbf{x}, t) \, d\mathbf{x}.
$$

Simply employing a first-order time discretization, the numerical scheme becomes

$$
U_{j,i}^* = U_{j,i}^n, \tag{C.7}
$$

$$
V_{j,i}^* = V_{j,i}^n - \frac{\Delta t}{\eta} (V_{j,i}^* - F(U_{j,i}^*)), \qquad (C.8)
$$

$$
W_{j,i}^* = W_{j,i}^n - \frac{\Delta t}{\eta} (W_{j,i}^* - K(U_{j,i}^*)), \qquad (C.9)
$$

$$
U_{j,i}^{n+1} = U_{j,i}^{*} - \frac{\Delta t}{\Delta x} (V_{j+\frac{1}{2},i}^{*} - V_{j-\frac{1}{2},i}^{*}) - \frac{\Delta t}{\Delta y} (W_{j,i+\frac{1}{2}}^{*} - W_{j,i-\frac{1}{2}}^{*}) + \Delta t G(U_{i,j}^{n+1}), \quad (C.10)
$$

$$
V_{j,i}^{n+1} = V_{j,i}^{*} - \frac{\Delta t}{\Delta x} A(U_{j+\frac{1}{2},i}^{*} - U_{j-\frac{1}{2},i}^{*}), \qquad (C.11)
$$

$$
W_{j,i}^{n+1} = W_{j,i}^{*} - \frac{\Delta t}{\Delta y} B(U_{j,i+\frac{1}{2}}^{*} - U_{j,i-\frac{1}{2}}^{*}), \qquad (C.12)
$$

$$
\frac{\bar{v}_{j,i}^{n+1} - v_{j,i}^n}{\Delta t} + \mathbf{v}_{j,i}^{n+1} \cdot (\nabla \bar{v}^{n+1})_{j,i} = g(c^n, \mu^{n+1}, \rho^{n+1})_{j,i}, \qquad (C.13)
$$

$$
g(c^n, \mu^{n+1}, \rho^{n+1})_{j,i} = \frac{1}{T'(v_{j,i}^n) \rho_{j,i}^{n+1} \Delta x} \left( (b(c^n) \nabla \mu^{n+1})_{j+\frac{1}{2},i} - (b(c^n) \nabla \mu^{n+1})_{j-\frac{1}{2},i} \right) \\
+ \frac{1}{T'(v_{j,i}^n) \rho_{j,i}^{n+1} \Delta y} \left( (b(c^n) \nabla \mu^{n+1})_{j,i+\frac{1}{2}} - (b(c^n) \nabla \mu^{n+1})_{j,i-\frac{1}{2}} \right) + \frac{F_c(\rho_{j,i}^n, c_{j,i}^n)}{T'(v_{j,i}^n) \rho_{j,i}^{n+1}}, \quad (C.14)
$$
$$ \mu_{j,i}^{n+1} = \frac{1}{\rho_{j,i}^n} \left( -\gamma T'(v_{j,i}^n) (\Delta \bar{v}^{n+1})_{j,i} - \gamma T''(v_{j,i}^n) |(\nabla v^n)_{j,i}|^2 \right) + \left( \frac{\partial \psi_0}{\partial c} \right)_{j,i}^n, \quad (C.15) $$

$$ \int_{\Omega} T(\lambda \bar{v}^{n+1}) dx = \int_{\Omega} c^n + \Delta t F_c dx, \quad (C.16) $$

$$ \bar{c}_{j,i}^{n+1} = T(\lambda_j \bar{v}_{j,i}^{n+1}), \quad (C.17) $$

$$ \frac{1}{\Delta t} (r^{n+1} - r^n) = - \frac{r^{n+1}}{E(\bar{c}^{n+1}) + C_0} \int_{\Omega} b(\bar{c}^{n+1}) |\nabla \mu^{n+1}|^2 d\mathbf{x} + \frac{r^{n+1}}{E(\bar{c}^{n+1}) + C_0} \int_{\Omega} \mu^{n+1} F_c(\rho^{n+1}, \bar{c}^{n+1}) d\mathbf{x}, \quad (C.18) $$

$$ \xi^{n+1} = \frac{r^{n+1}}{E(\bar{c}^{n+1}) + C_0}, \quad (C.19) $$

$$ c_{j,i}^{n+1} = \nu^{n+1} \bar{c}_{j,i}^{n+1}, \quad \text{with} \quad \nu^{n+1} = 1 - (1 - \xi^{n+1})^2, \quad (C.20) $$

$$ v_{j,i}^{n+1} = \nu^{n+1} \bar{v}_{j,i}^{n+1}. \quad (C.21) $$

## 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. A.P. acknowledges support from the Labex CEMPI (ANR-11-LABX-0007-01).

## REFERENCES

[1] H. Abels, On a diffuse interface model for two-phase flows of viscous, incompressible fluids with matched densities. *Arch. Ration. Mech. Anal.* **194** (2009) 463–506.

[2] H. Abels, Strong well-posedness of a diffuse interface model for a viscous, quasi-incompressible two-phase flow. *SIAM J. Math. Anal.* **44** (2012) 316–340.

[3] H. Abels and E. Feireisl, On a diffuse interface model for a two-phase flow of compressible viscous fluids. *Indiana Univ. Math. J.* **57** (2008) 659–698.

[4] H. Abels, H. Garcke and G. Grün, Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities. *Math. Models Methods Appl. Sci.* **22** (2012) 40.

[5] H. Abels, D. Depner and H. Garcke, On an incompressible Navier-Stokes/Cahn-Hilliard system with degenerate mobility. *Ann. Inst. H. Poincaré C Anal. Non Linéaire* **30** (2013) 1175–1190.

[6] A. Agosti, Discontinuous Galerkin finite element discretization of a degenerate Cahn-Hilliard equation with a single-well potential. *Calcolo* **56** (2019) 47.

[7] A. Agosti, P.F. Antonietti, P. Ciarletta, M. Grasselli and M. Verani, A Cahn-Hilliard-type equation with application to tumor growth dynamics. *Math. Methods Appl. Sci.* **40** (2017) 7598–7626.

[8] G.L. Aki, W. Dreyer, J. Giesselmann and C. Kraus, A quasi-incompressible diffuse interface model with phase transition. *Math. Models Methods Appl. Sci.* **24** (2014) 827–861.

[9] D. Ambrosi and L. Preziosi, On the closure of mass balance models for tumor growth. *Math. Models Methods Appl. Sci.* **12** (2002) 737–754.

[10] D.M. Anderson, G.B. McFadden and AA. Wheeler, Diffuse-interface methods in fluid mechanics. *Annu. Rev. Fluid Mech.* **30** (1998) 139–165.

[11] X. Antoine, J. Shen and Q. Tang, Scalar auxiliary variable/Lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear Schrödinger/Gross-Pitaevskii equations. *J. Comput. Phys.* **437** (2021) 19.

[12] J.W. Barrett, J.F. Blowey and H. Garcke, Finite element approximation of the Cahn-Hilliard equation with degenerate mobility. *SIAM J. Numer. Anal.* **37** (1999) 286–318.

[13] D. Bresch and B. Desjardins, Existence of global weak solutions for a 2D viscous shallow water equations and convergence to the quasi-geostrophic model. *Commun. Math. Phys.* **238** (2003) 211–223.
[14] D. Bresch and B. Desjardins, On the construction of approximate solutions for the 2D viscous shallow water model and for compressible Navier-Stokes models. *J. Math. Pures Appl.* **86** (2006) 362-368.

[15] D. Bresch, A.F. Vasseur and C. Yu, Global existence of entropy-weak solutions to the compressible Navier-Stokes equations with non-linear density dependent viscosities. *J. Eur. Math. Soc. (JEMS)* **24** (2022) 1791–1837.

[16] H. Byrne and L. Preziosi, Modelling solid tumour growth using the theory of mixtures. *Math. Med. Biol.* **20** (2003) 341–366.

[17] J.W. Cahn, On spinodal decomposition. *Acta metall.* **9** (1961) 795–801.

[18] J.W. Cahn and J.E. Hilliard, Free energy of a nonuniform system. I. interfacial free energy. *J. Chem. Phys.* **28** (1958) 258–267.

[19] C. Chatelain, T. Balois, P. Ciarletta and M.B. Amar, Emergence of microstructural patterns in skin cancer: a phase separation analysis in a binary mixture. *New J. Phys.* **13** (2011) 115013.

[20] C. Chatelain, P. Ciarletta and M. Ben Amar, Morphological changes in early melanoma development: influence of nutrients, growth inhibitors and cell-adhesion mechanisms. *J. Theoret. Biol.* **290** (2011) 46–59.

[21] L. Cherfils, E. Feireisl, M. Michálek, A. Miranville, M. Petcu and D. Pražák, The compressible Navier-Stokes-Cahn-Hilliard equations with dynamic boundary conditions. *Math. Models Methods Appl. Sci.* **29** (2019) 2557–2584.

[22] M. Ebenbeck, H. Garcke and R. Nürnberg, Cahn-Hilliard-Brinkman systems for tumour growth. *Discrete Cont. Dyn. S* **14** (2021) 3989–4033.

[23] C. Eck, H. Garcke and P. Knabner, Mathematical Modeling. *Undergraduate Mathematics Series*. Springer, Cham (2017).

[24] C.M. Elliott and H. Garcke, On the Cahn-Hilliard equation with degenerate mobility. *SIAM J. Math. Anal.* **27** (1996) 404–423.

[25] J.L. Ericksen, Liquid crystals with variable degree of orientation. *Arch. Ration. Mech. Anal.* **113** (1990) 97–120.

[26] R. Eymard, T. Gallouët and R. Herbin, Finite volume methods, in Handbook of Numerical Analysis. Vol. VII. North-Holland, Amsterdam (2000) 713–1020.

[27] E. Feireisl, Dynamics of Viscous Compressible Fluids. *Oxford Lecture Series in Mathematics and its Applications*. Vol. 26. Oxford University Press, Oxford (2004).

[28] E. Feireisl and P. Laurençot, Non-isothermal Smoluchowski-Poisson equations as a singular limit of the Navier-Stokes-Fourier-Poisson system. *J. Math. Pures Appl.* **88** (2007) 325–349.

[29] E. Feireisl and A. Novotný, Singular limits in thermodynamics of viscous fluids, in Advances in Mathematical Fluid Mechanics. Birkhäuser Verlag, Basel (2009).

[30] M. Fritz, Tumor evolution models of phase-field type with nonlocal effects and angiogenesis. *Bull. Math. Biol.* **85** (2023) 34.

[31] H. Garcke, K.F. Lam, E. Sitka and V. Styles, A Cahn-Hilliard-Darcy model for tumour growth with chemotaxis and active transport. *Math. Models Methods Appl. Sci.* **26** (2016) 1095–1148.

[32] J. Giesselmann and T. Pryer, Energy consistent discontinuous Galerkin methods for a quasi-incompressible diffuse two phase flow model. *ESAIM Math. Model. Numer. Anal.* **49** (2015) 275–301.

[33] Z. Guo, P. Lin and J.S. Lowengrub, A numerical method for the quasi-incompressible Cahn-Hilliard-Navier-Stokes equations for variable density flows with a discrete energy law. *J. Comput. Phys.* **276** (2014) 486–507.

[34] M.E. Gurtin, On a nonequilibrium thermodynamics of capillarity and phase. *Q. Appl. Math.* **47** (1989) 129–145.

[35] M.E. Gurtin, D. Polignone and J. Viñals, Two-phase binary fluids and immiscible fluids described by an order parameter. *Math. Models Methods Appl. Sci.* **6** (1996) 815–831.

[36] M.E. Gurtin, E. Fried and L. Anand, The Mechanics and Thermodynamics of Continua. Cambridge University Press, Cambridge (2010).

[37] Q. He and X. Shi, Numerical study of compressible Navier-Stokes-Cahn-Hilliard system. *Commun. Math. Sci.* **18** (2020) 571–591.

[38] P.C. Hohenberg and B.I. Halperin, Theory of dynamic critical phenomena. *Rev. Mod. Phys.* **49** (1977) 435–479.

[39] B.S. Hosseini, S. Turek, M. Möller and C. Palmes, Isogeometric analysis of the Navier-Stokes-Cahn-Hilliard equations with application to incompressible two-phase flows. *J. Comput. Phys.* **348** (2017) 171–194.

[40] F. Huang and J. Shen, Bound/positivity preserving and energy stable scalar auxiliary variable schemes for dissipative systems: applications to Keller-Segel and Poisson-Nernst-Planck equations. *SIAM J. Sci. Comput.* **43** (2021) A1832–A1857.

[41] F. Huang and J. Shen, A new class of implicit-explicit BDFk SAV schemes for general dissipative systems and their error analysis. *Comput. Methods Appl. Mech. Eng.* **392** (2022) 25.
[42] F. Huang, J. Shen and K. Wu, Bound/positivity preserving and unconditionally stable schemes for a class of fourth order nonlinear equations. *J. Comput. Phys.* **460** (2022) 16.

[43] M. Jiang, Z. Zhang and J. Zhao, Improving the accuracy and consistency of the scalar auxiliary variable (SAV) method with relaxation. *J. Comput. Phys.* **456** (2022) 20.

[44] S. Jin and Z.P. Xin, The relaxation schemes for systems of conservation laws in arbitrary space dimensions. *Commun. Pure Appl. Math.* **48** (1995) 235-276.

[45] K.F. Lam and H. Wu, Thermodynamically consistent Navier-Stokes-Cahn-Hilliard models with mass transfer and chemotaxis. *Eur. J. Appl. Math.* **29** (2018) 595-644.

[46] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems. Vol. 31. Cambridge University Press, Cambridge (2002).

[47] X. Li and J. Shen, On a SAV-MAC scheme for the Cahn-Hilliard-Navier-Stokes phase-field model and its error analysis for the corresponding Cahn-Hilliard-Stokes case. *Math. Models Methods Appl. Sci.* **30** (2020) 2263-2297.

[48] P.-L. Lions, Mathematical Topics in Fluid Mechanics: Vol. 2: Compressible Models. Vol. 10 of *Oxford Lecture Series in Mathematics and its Applications*. The Clarendon Press, Oxford University Press, New York (1998).

[49] S. Lisini, D. Matthes and G. Savaré, Cahn-Hilliard and thin film equations with nonlinear mobility as gradient flows in weighted-Wasserstein metrics. *J. Differ. Equ.* **253** (2012) 814-850.

[50] I.S. Liu, Method of Lagrange multipliers for exploitation of the entropy principle. *Arch. Rational Mech. Anal.* **46** (1972) 131-148.

[51] J. Lowengrub and L. Truskinovsky, Quasi-incompressible Cahn-Hilliard fluids and topological transitions. *R. Soc. Lond. Proc. Ser. A Math. Phys. Eng. Sci.* **454** (1998) 2617-2654.

[52] J. Lowengrub, E. Titi and K. Zhao, Analysis of a mixture model of tumor growth. *Eur. J. Appl. Math.* **24** (2013) 691-734.

[53] A. Miranville, The Cahn-Hilliard Equation: Recent Advances and Applications. Vol. 95 of *CBMS-NSF Regional Conference Series in Applied Mathematics*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia (2019).

[54] A. Moussa, Some variants of the classical Aubin-Lions lemma. *J. Evol. Equ.* **16** (2016) 65-93.

[55] I. Müller, Thermodynamics. Pitman, Boston (1985).

[56] G.A. Narsilio, O. Buzzi, S. Fityus, T.S. Yun and D.W. Smith, Upscaling of Navier-Stokes equations in porous media: theoretical, numerical and experimental approach. *Comput. Geotechn.* **36** (2009) 1200-1206.

[57] B. Perthame and A. Poulain, Relaxation of the Cahn-Hilliard equation with singular single-well potential and degenerate mobility. *Eur. J. Appl. Math.* **32** (2021) 89-112.

[58] A. Poulain and K. Schratz, Convergence, error analysis and longtime behavior of the scalar auxiliary variable method for the nonlinear Schrödinger equation. *IMA J. Numer. Anal.* **42** (2022) 2853-2883.

[59] E. Rocca, G. Schimperna and A. Signori, On a Cahn-Hilliard-Keller-Segel model with generalized logistic source describing tumor growth. *J. Differ. Equ.* **343** (2023) 530-578.

[60] J. Shen, J. Xu and J. Yang, The scalar auxiliary variable (SAV) approach for gradient flows. *J. Comput. Phys.* **353** (2018) 407-416.

[61] J. Shen, J. Xu and J. Yang, A new class of efficient and robust energy stable schemes for gradient flows. *SIAM Rev.* **61** (2019) 474-506.

[62] J. Shen and X. Yang, The IEQ and SAV approaches and their extensions for a class of highly nonlinear gradient flow systems. *Contemp. Math.* **754** (2020) 217-245.

[63] M.F.P. ten Eikelder, K.G. van der Zee, I. Akkerman and D. Schillinger, A unified framework for Navier-Stokes Cahn-Hilliard models with non-matching densities. *Math. Models Methods Appl. Sci.* **33** (2023) 175-221.

[64] A.F. Vasseur and C. Yu, Existence of global weak solutions for 3D degenerate compressible Navier-Stokes equations. *Invent. Math.* **206** (2016) 935-974.

[65] S.M. Wise, J.S. Lowengrub, H.B. Frieboes and V. Cristini, Three-dimensional multispecies nonlinear tumor growth - I: model and numerical method. *J. Theor. Biol.* **253** (2008) 524-543.

[66] Y. Zhang and J. Shen, A generalized SAV approach with relaxation for dissipative systems. *J. Comput. Phys.* **464** (2022) 23.
<figure><img src="image_11.png" /></figure>

Please help to maintain this journal in open access!

This journal is currently published in open access under the [Subscribe to Open model](https://www.eds Sciences.org/). We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting [subscribers@edpsciences.org](mailto:subscribers@edpsciences.org).

More information, including a list of supporters and financial transparency reports, is available at [https://edpsciences.org/en/subscribe-to-open-s2o](https://edpsciences.org/en/subscribe-to-open-s2o).