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

HAL
open science

# On the Incompressible Limit for a Tumour Growth Model incorporating Convective Effects

Noemi David, Markus Schmidtchen

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

## To cite this version:

Noemi David, Markus Schmidtchen. On the Incompressible Limit for a Tumour Growth Model incorporating Convective Effects. Communications on Pure and Applied Mathematics, 2023. hal-03162169v2

HAL Id: hal-03162169

https://hal.science/hal-03162169v2

Submitted on 28 Dec 2022

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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ON THE INCOMPRESSIBLE LIMIT FOR A TUMOUR GROWTH MODEL
INCORPORATING CONVECTIVE EFFECTS

NOEMI DAVID<fnref n="1" />,<fnref n="2" /> AND MARKUS SCHMIDTCHEN<fnref n="3" />

ABSTRACT. In this work we study a tissue growth model with applications to tumour growth. The model is based on that of Perthame, Quirós, and Vázquez proposed in 2014 but incorporates the advective effects caused, for instance, by the presence of nutrients, oxygen, or, possibly, as a result of self-propulsion. The main result of this work is the incompressible limit of this model which builds a bridge between the density-based model and a geometry free-boundary problem by passing to a singular limit in the pressure law. The limiting objects are then proven to be unique.

2010 Mathematics Subject Classification. 35B45; 35K57; 35K65; 35Q92; 76N10; 76T99;
Keywords and phrases. Porous medium equation; Tumour growth; Aronson-Bénilan estimate;
Incompressible limit; Free boundary; Hele-Shaw problem.

# 1 INTRODUCTION

Modelling living tissue poses a whole range of challenges. On the one hand, it is important to identify the biomedical drivers that should be incorporated in the model, while, on the other hand there are certain modelling choices that need to be discussed. One of these choices that, in a way, separates the community is the type of model used to describe tissue growth. Roughly speaking we identify the following two types of models: those that describe the tissue as an evolving distribution in space and those that describe the tissue as an evolving domain in space. While the first type is mostly based on a partial differential equation description, the latter is known as a free-boundary or evolving boundary model.

The goal of this paper is to build a bridge between the two types of models by passing to the so-called stiff limit in the population-based model to obtain a free-boundary description. The model we propose here describes the evolution of the tissue density, $n_γ = n_γ(x, t)$, and is given by

$$ (1) \qquad \frac{\partial n_\gamma}{\partial t} - \nabla \cdot (n_\gamma \nabla p_\gamma) - \nabla \cdot (n_\gamma \nabla \Phi) = n_\gamma G(p_\gamma). $$

on $\mathbb{R}^d$ and for $t > 0$. It is equipped with some non-negative initial data $n_γ(0, x) = n_γ^0(x) \in L_+^1(\mathbb{R}^d)$. Here $p_γ = n_γ^γ$ denotes the pressure, $G = G(p_γ)$ models the cell proliferation (resp. cell death), and $\Phi = \Phi(x, t)$ denotes a chemical concentration. In order to pass to the incompressible limit $\gamma \to \infty$ we need to study the equation satisfied by the pressure, i.e., in the equation

$$ (2) \qquad \frac{\partial p_\gamma}{\partial t} = \gamma p_\gamma (\Delta p_\gamma + \Delta \Phi + G(p_\gamma)) + \nabla p_\gamma \cdot \nabla (p_\gamma + \Phi). $$

<footnote n="1">¹ SORBONNE UNIVERSITÉ, INRIA, UNIVERSITÉ DE PARIS, LABORATOIRE JACQUES-LOUIS LIONS UMR 7598, PARIS 75005, FRANCE (NOEMI.DAVID@LJLL.MATHUPMC.FR).</footnote>

<footnote n="2">² DIPARTIMENTO DI MATEMATICA, UNIVERSITÀ DI BOLOGNA, ITALY</footnote>

<footnote n="3">³ INSTITUTE FOR SCIENTIFIC COMPUTING, TECHNISCHE UNIVERSITÄT DRESDEN, ZELLESCHER WEG 12-14, 01069 DRESDEN, GERMANY. (MARKUS.SCHMIDTCHEN@TU-DRESDEN.DE).</footnote>
While it is intuitive to expect

$$p_{\infty}(\Delta p_{\infty} + \Delta\Phi + G(p_{\infty})) = 0, \quad \text{as well as} \quad p_{\infty}(n_{\infty} - 1) = 0,$$

in the limit, there are technical subtleties, obtaining strong compactness of the pressure gradient to be precise, that need to be overcome. We are by no means the first to ask this question. As a matter of fact, there are already some promising results towards this rigorous limit. However, all of them are borderline and just not good enough to obtain the strong compactness of the pressure gradient. A blend of two techniques finally allows us to settle this open question. The rest of the introduction is dedicated to presenting a historical view on this type of model as well as variations thereof. We will also use this as an opportunity to introduce the tools necessary for the limit passage in a brief, explanatory way.

## 1.1 Historical Notes – The Origin of Incompressible Limits & the Mesa Problem

The question of passing to the incompressible limit has a rich history and several variations of it have been studied in the literature. Historically, the problem has its early foundation in the work of Bénilan and Crandall on the continuous dependence on $\varphi$ of solutions to the filtration equation $\partial_t n = \Delta\varphi(n)$ in 1981, cf. [4], not too long after the first wellposedness results for the filtration equation around 1960, cf. [48, 54]. The continuous dependence of [4] is established using non-linear m-accretive semi-group theory, notably maximal monotone operators, enabling them to allow for cases of $\varphi$ being a monotone graph. As a matter of fact, it already covers the first result on incompressible limits by choosing $\varphi(z) = z^{\gamma}$ and assuming non-negative initial data bounded from above by unity. Henceforth the problem has been attracting a lot of attention.

In [30] the authors show the formation of a plateau-like region, which they refer to as ‘mesa’, of nearly constant density $n_{\gamma}$, for $\gamma \approx \infty$, using a formal asymptotic expansions and working with radial solutions. In [13], too, the authors consider the limit of the density of the porous equation but they can weaken the assumption on the initial data thus extending the results of [4]. Moreover, they are able to show that the limit density, $n_{\infty}$, is independent of time and bounded $0 \le n_{\infty} \le 1$. This ‘stationarity’ result on the limit density, $n_{\infty}$, is obtained upon combining three tools. First, the uniform essential bounds on the compactly supported densities, $n_{\gamma}$ imply the weak-star convergence of a subsequence. Second, by the classical Aronson-Bénilan estimate (see [2] for the original article as well as [7] and references therein for a survey), it can be inferred that $\partial_t n_{\infty} \ge 0$, and therefore $n_{\infty}(x, t) \ge n_{\infty}(s, x)$ for almost every $x \in \mathbb{R}^d$, $s < t$, and all $\gamma > 1$. Finally, the conservation of mass implies that, in fact, $n_{\infty}(x, t) = n_{\infty}(s, x)$, which shows that $n_{\infty}$ is independent of time, cf. [13] for the full argument.

Later, in 2001, Gil and Quirós revisit the study of the incompressible limit of the solution of the porous medium equation defined in $[0, +\infty) \times \Omega$. In their paper they prove that the solution of the porous medium equation converges to that of the Hele-Shaw problem in the sense of Elliot and Janovsky, i.e., in the form of a variational formulation whenever the boundary data $g = g(x)$ is independent of time and the initial data is the indicator function of some bounded set $\Omega_0 \subset \Omega$. In this case, the weak formulation and the variational formulation coincide, cf. [31, Corollary 4.5]). In their study, cf. [31], $\Omega$ is assumed to be a compact subset of $\mathbb{R}^d$ which is equipped with Dirichlet data on the pressure, $p_{\gamma}(x, t) = g(x, t)$ on $\partial\Omega$, for some $g(x, t) \ge 0$. Let us point out that, given a set $\Omega$ large enough, the case $g \equiv 0$ coincides with the problem studied by Caffarelli and Friedman in [13], and, again, the limit is independent of time. Indeed, Gil and Quirós are able to recover the same result from a different perspective, focusing on the role of the pressure rather than the density itself. In the absence of Dirichlet boundary data, i.e., $g \equiv 0$, the limit solution solves a Hele-Shaw problem where the free boundary is actually motionless since the limit pressure vanishes almost everywhere. This can be easily seen by passing to the limit $\gamma \to \infty$ in the porous medium pressure equation, Eq. (2), where, of course, the growth term and the migration term are absent. In conjunction with the uniform essential bounds this immediately yields $\|\nabla p_{\infty}\|_{L^2(\Omega \times (0,T))} = 0$.

On the other hand, in the case non-vanishing $g \ge 0$ on $\partial\Omega$, the pressure is “forced” to be positive near to the boundary, and then, since the pressure gradient is no longer zero, the motion of the
free boundary $\partial\{p_{\infty} > 0\}$ is governed by Darcy's law

$$V = -\partial_{\nu}p_{\infty},$$

where $\nu$ denotes the outward normal on the free boundary. In [32] the authors generalise there result towards a broader class of initial data give a description of the positivity set of the densities, $n_{\gamma}$, to that of the limit.

Let us also stress that the conservation of mass no longer holds since there is a source term on the boundary of $\Omega$. Therefore, the proof of the stationarity of $n_{\infty}$ using the Aronson-Bénilan estimate fails. Similarly, the proof of $\|\nabla p_{\infty}\|_{L^2} = 0$, no longer holds true due to the fact that the boundary terms arising from integration by parts no longer vanish.

It is also worthwhile noticing that $p_{\gamma} \approx n_{\gamma}p_{\gamma}$, for $\gamma \gg 1$, which leads to the relation

$$p_{\infty}(1 - n_{\infty}) = 0.$$

Hence, we infer the inclusion $\{p_{\infty} > 0\} \subset \{n_{\infty} = 1\}$, but we also stress that the two sets need not coincide. In fact, in the case $g = 0$, or equivalently the porous medium equation on $\mathbb{R}^d$ with compactly supported initial data, as mentioned above, the limiting pressure vanishes, $p_{\infty} = 0$, almost everywhere and the limit density is stationary, $n_{\infty}(x, t) = n^0(x)$, where $0 \le n^0(x) \le 1$. This means that, even if there are saturation zones, $\{n_{\infty} = 1\}$, the pressure does not become positive. This situation changes drastically if the model includes a positive growth term of the form

$$\frac{\partial n_{\gamma}}{\partial t} - \nabla \cdot (n_{\gamma} \nabla p_{\gamma}) = n_{\gamma} G(p_{\gamma}),$$

as was proposed in [50]. In this case it can be shown that the two sets coincide, i.e., $\{p_{\infty} > 0\} = \{n_{\infty} = 1\}$, and, what is more, the problem is no longer stationary!

## 1.2 Contemporary Advances – Generalisations of the Model

Emanating from the early works on the mesa problem for the porous medium equation, research began branching out in different directions. In this section we aim at giving a brief overview of different extensions of the porous medium equation, applications of the models obtained this way, as well as techniques used to study their respective incompressible limits analytically.

The first generalisation concerns the inclusion of a pressure-dependent growth term proposed in the work of [50]. Here the authors propose a tissue-growth model where cells move according to a population pressure generated by the total density of the form $p(n) = n^{\gamma}$. In conjunction with Darcy's law they recover the porous-medium type degenerate diffusion. In addition, they include a proliferation term, $nG(p)$, which models cells divisions with a pressure depending rate. Thus the proliferation rate, $G$, is assumed to be a decreasing function accounting for the fact that cells are less 'willing' to divide in packed regimes, cf. Section 1.2.1.

The model was then extended by a nutrient distribution, $c(x, t)$, which is assumed to diffuse in the domain and released (resp. decayed) by general $L^2$-processes, cf. Section 1.2.2. Most recently, the inclusion of migratory processes, i.e., drift terms given by a velocity field, $v(x, t)$, as a model extension received a lot of attention, cf. Section 1.2.3. This is also where our contribution to the current discourse enters, namely the first rigorous derivation of the complementarity relation, that is, an equation governing the pressure distribution inside of the moving boundary problem. Before we begin discussing our main result we shall also point out recent advances in the area of stiff-limits in the context of pressure laws that are different from Darcy's law, cf. Section 1.2.4. We conclude our short survey of the literature by mentioning some multi-phase results, where, instead of one equation, two interacting species are considered, cf. Section 1.2.5.

### 1.2.1 A Model including Proliferation

In [50], Perthame, Quirós, and Vázquez propose the model

$$(3) \qquad \frac{\partial n_{\gamma}}{\partial t} - \nabla \cdot (n_{\gamma} \nabla p_{\gamma}) = n_{\gamma} G(p_{\gamma}).$$

Their paper is seminal in that they were the first to perform the rigorous stiff pressure limit in the presence of growth terms. While strong compactness of the pressure is absolutely sufficient
for the Hele-Shaw limit itself, obtaining the so-called complementarity relation which provides an equation for the pressure in the limit is much more involved. In fact, in order to obtain it strong compactness of the pressure gradient is indispensable. To this purpose, using the comparison principle, they show that the Laplacian of the pressure satisfies an Aronson-Bénilan type estimate, namely $\Delta p + G(p) \gtrsim -C/\gamma t$.

In [39] the authors study the same model through a viscosity solution approach. They are able to show that the density converges locally uniformly away from the free boundary $\partial\{p_{\infty} > 0\}$. Moreover, they prove locally uniform convergence of the pressure (as long as the limit is continuous) and that $p_{\infty}$ is the viscosity solution of the Hele-Shaw problem

$$ (4) \qquad \left\{ \begin{array}{ll} -\Delta p_{\infty} = G(p_{\infty}), & \text{in } \{p_{\infty} > 0\}, \\ V = -\frac{|\nabla p_{\infty}|}{1 - \min(1, n_{\infty}^{E})}, & \text{on } \partial\{p_{\infty} > 0\}, \end{array} \right. $$

where the normal velocity law was only formally presumed in [50], but not rigorously proven. Here, $n_{\infty}^E$ denotes the trace of the “external” limit density on the free boundary, namely the trace of $n_{\infty}$ from the set $\{n_{\infty} < 1\}$.

Let us stress the fact that, as the velocity law suggests, the density shows jump discontinuities at the free boundary. Moreover, the velocity blows up when the density reaches value 1, therefore, when new *mesas* appear outside of $\{p_{\infty} > 0\}$, the pressure becomes instantaneously positive in the new nucleated regions, hence exhibiting time discontinuities.

The free boundary problem, Eq. (4), was further studied in [47], where the authors prove that the velocity law of the free boundary holds both in a weak (distributional) and in a measure theoretical sense. In the same paper, they also provide an $L^4$-bound of the pressure gradient that relies on the Aronson-Bénilan estimate, which we extend to our model, Eq. (1), through a self-contained proof in Lemma 3.2, independently of any estimate on $\Delta p_{\gamma}$. A different approach for the incompressible limit for Eq. (3) was taken in [16], where a transport-growth distance is introduced such that Eq. (3) can be understood as a gradient flow with respect to said metric.

### 1.2.2 A Model including Nutrients

In [50], the authors also study an extension of the model including the effect of a nutrient with concentration $c = c(x, t)$ in the growth term

$$ \frac{\partial n_{\gamma}}{\partial t} - \nabla \cdot (n_{\gamma} \nabla p_{\gamma}) = n_{\gamma} G(p_{\gamma}, c_{\gamma}). $$

While they were able to prove the strong convergence of $n_{\gamma}$ and $c_{\gamma}$ as $\gamma \to \infty$, they leave open the question of how to recover the $L^2$-strong compactness of the pressure gradient needed to pass to the limit in the pressure equation and obtain the complementarity relation.

This problem was addressed in [22], where the authors combine a weak version of the Aronson-Bénilan estimate in $L^3$ with a uniform bound of the pressure gradient in $L^4$ to infer strong compactness. In fact, the $L^{\infty}$-Aronson-Bénilan estimate does not hold in the nutrient case, since $G(p, c)$ can be negative and then the comparison principle used in [50] fails. Travelling waves solutions of the Hele-Shaw problem that arises in the stiff limit have been studied in [20, 52]. Besides, explicit solutions to the limit problem are presented in [44] for initial data of the form of an indicator of a bounded set. Recently, interesting progress have been made in [34] where the authors are able to establish the incompressible limit and the complementarity relation without relying on any Aronson-Bénilan-type estimates. Instead, their approach is based on viscosity solutions and establishing the equivalence between the complementarity relation and an obstacle problem.

### 1.2.3 Models including local and non-local Drifts

In 2010, Kim and Lei introduced the notion of viscosity solution for the porous medium equation with drift

$$ \frac{\partial n_{\gamma}}{\partial t} = \Delta n_{\gamma} + \nabla \cdot (n_{\gamma} \nabla \Phi), $$
and they prove that it coincides with the weak solution in the distributional sense, cf. [38]. Using the same viscosity approach, in [1] the authors study the link between the Hele-Shaw model with drift

$$ \begin{cases} -\Delta p = \Delta \Phi, & \text{in } \{p > 0\}, \\ V = -(\nabla p + \nabla \Phi) \cdot \nu, & \text{on } \partial\{p > 0\}, \end{cases} $$

and the congested crowd motion model

$$ \partial_t n + \nabla \cdot (n \nabla \Phi) = 0, $$

if $n < 1$, with $n \le 1$, where the latter constraint comes from the singular limit in the nonlinear diffusion term. To prove the equivalence of the two models, they study the asymptotics of the porous medium equation with drift as $\gamma \to \infty$. They show that the viscosity solution converges locally uniformly to a solution of the Hele-Shaw model. At the same time, using the metric setting of the 2-Wasserstein space, they infer the convergence to the aforementioned congested crowd motion model. To this purpose, they assume the potential $\Phi$ to be sub-harmonic, i.e., $\Delta\Phi > 0$. While the convergence in the 2-Wasserstein distance holds for general initial data $0 \le n_0 \le 1$, the locally uniform limit holds only for patches, namely $n^0 = 1_{\Omega_0}$, with $\Omega_0$ a compact set in $\mathbb{R}^d$. This result was extended in 2016, by Craig, Kim, and Yao, cf. [17] to a model with non-local Newtonian potential, $\mathcal{N}$,

$$ \frac{\partial n_{\gamma}}{\partial t} = \Delta n_{\gamma} + \nabla \cdot (n_{\gamma} \nabla \mathcal{N} \star n_{\gamma}). $$

The main novelty they introduce is that they are able to study the incompressible limit despite the lack of convexity. In fact, unlike the congested drift equation studied in [1], the energy related to the aggregation equation through the 2-Wasserstein gradient flow structure is not semi-convex, cf. [17]. Even more recently, the $\Gamma$-limit to obtain the incompressible limit has been studied in [18] for a wider class of interaction potentials, and in [15] at the level of the stationary states, cf. also [14] and references therein.

The question of how to pass to the limit $\gamma \to \infty$ in the porous medium equation with a drift and a non-trivial source term has been addressed in [40]. The authors propose a model with a generic vector field $v: \mathbb{R}^d \times \mathbb{R}^+ \to \mathbb{R}^d$ as drift term, i.e.,

$$ \frac{\partial n_{\gamma}}{\partial t} - \Delta n_{\gamma} + \nabla \cdot (n_{\gamma} v) = n_{\gamma} G, $$

with a growth rate $G = G(x, t)$. Through viscosity solutions methods, they prove that as $\gamma \to \infty$ the model converges to a free boundary model of Hele-Shaw type. Their work improves the results previously achieved in [1], extending the class of initial data from patches to any continuous and compactly supported function bounded between zero and one. In the absence of any growth dynamics, the rate of convergence as $\gamma \to \infty$ in the Wasserstein distance was obtained in [1] and was recently improved (in an $H^{-1}$ sense) by [21] who also allow for growth dynamics.

### 1.2.4 Different Pressure Laws and Relations

As foreshadowed above, in certain contexts Darcy's law may not be the appropriate relation that links the velocity field to the mechanical pressure. Depending on the modelling context and the model complexity, the pressure is incorporated in the fluid velocity through Stokes flow, Brinkman's law or Navier-Stokes' law, rather than Darcy's law. We briefly present recent works of incompressible limits for different pressure laws and relations.

#### Singular Pressure

Parallel to the advances in the context of incompressible limits with power-law pressures it has been observed that another pressure law of the form

$$ (5) \qquad p_{\epsilon}(n) = \epsilon \frac{n}{1-n}, $$
can be used to model living tissue, cf. [36]. Using this singular pressure law already introduces an incompressibility condition in the sense that the pressure blows up when the cell density reaches the saturated regime, $n = 1$. Thus, singular pressure laws of this kind are encountered in scenarios when non-overlap conditions are enforced already at a population-level, cf. [27, 49] in the context of congestive collective crowd motion, [5, 6] in the context of traffic flow modelling. In [36] the authors are able to show that the pressure in Eq. (5) is suitable to pass to the incompressible limit using a generalisation of the Aronson-Bénilan argument by Crandall and Pierre, cf. [19].

## Brinkman Law Pressure

Unlike Darcy's law using the Brinkman law,

$$-\nu\Delta W + W = p(n),$$

accounts for visco-elastic effects, [12]. Based on this observation, in [53] the authors propose a modification of the above model, Eq. (3), incorporating the Brinkman law, i.e.,

$$\frac{\partial n_{\gamma}}{\partial t} - \nabla \cdot (n_{\gamma} \nabla W_{\gamma}) = n_{\gamma} G(p_{\gamma}).$$

Different from the Darcy law setting the authors are forced to use a different set of techniques since the problem is no longer degenerate parabolic but, instead, of transport nature. While, at first glance, the Brinkman law has a regularising effect on the velocity field it makes obtaining compactness of the pressure a hard endeavour. Using a kinetic reformulation and controlling oscillations in the pressure finally yields the required compactness to pass to the incompressible limit and obtain a visco-elastic version of the complementarity relation, cf. [53, Theorem 1.1]. For pressure laws of the form $p_{\epsilon}(n) = \epsilon \mathbb{1}_{n \ge 1} \log(n)$, quite recently, explicit travelling wave profiles we obtained by [45]. Moreover, the authors provide an apt numerical scheme to track the moving front accurately.

## Stokes Flow

It is important to stress that both Darcy's law and Brinkman's law are, at least, formally related to the Navier-Stokes law which can therefore be seen as the most general relation between the fluid velocity and the mechanical pressure. In [56] the authors prove the incompressible limit for a proliferating species whose velocity is linked to the pressure through the Navier-Stokes law thus generalising the case without birth and death processes of [43]. The authors use the fact that the growth rate is linear in the pressure such that weak compactness of the pressure suffices in order to pass to the limit, so long as the density itself is strongly compact. While the weak compactness of the pressure follows from a renormalisation argument the strong compactness of the density is based on a compactness-propagation argument introduced (and later refined) in [3, 9, 10].

## Active Motion

In [51] the authors extend the model of [50] by an additional active motion term in form of a linear diffusion term. They are able to rigorously perform the incompressible limit, in fact they obtain the same complementarity relation as in the absence of active motion, for certain initial data not relying on the Aronson-Bénilan for certain initial data. Nonetheless, the restriction on the initial data can be dropped by employing the argument of Crandall and Pierre, in [19]. In [55] the authors propose a very similar model based on Brinkman's law (unlike [51]) including a linear diffusion term. They observe that travelling waves exist and analyse their profile.

## Fractional Diffusion

In 2015, J.-L. Vázquez opened another both fascinating and challenging research direction by addressing the mesa problem in the fractional pressure case, cf. [57]. More precisely, he studies
the incompressible limit, $\gamma \to \infty$, in the fractional porous medium equation,

$$ \frac{\partial n_{\gamma}}{\partial t} + (-\Delta)^{-s}(n_{\gamma})^{\gamma} = 0, $$

for $s \in (0, 1)$. Unlike the case of classical porous medium type diffusion, the limiting profile exhibits tails and does not remain compactly supported. The analysis is of orders of magnitude harder since the classical theory discussed in Section 1.1 relies on comparison principles and the fact that it is known what happens to the Barenblatt profiles in the incompressible limit. In the fractional setting the explicit source solutions are not known explicitly. None the less, they are the starting point of the analysis of [57]. Many questions remain open, in particular the inclusion of other processes such as reactions and drifts.

### 1.2.5 Multi-Species System

Recently, there has growing interest in multi-phase extensions of the above model. Instead of merely modelling the evolution of a single species, say, cancer tissue, other phases such as interstitial fluid, healthy tissue, dead tissue, ..., are incorporated into the model. The extension to multiple interacting species not only leads to interesting behaviours such as phase separation but also raises novel mathematical challenges such as the loss of regularity at so-called internal layers, i.e., regions where two or more phases get in contact. Recently, [11] have established the rigorous incompressible limit for a two-species model consisting of normal and abnormal tissue, respectively for a Darcy law type pressure. Unlike in the single-species case, the pressure is now generated by the joint population in form of a power law. However, the lack of regularity is such that only a one dimensional result could be obtained and the general case was successfully addressed only recently, cf. [46]. In a similar fashion, a one-dimensional result could be obtained, see [26], when the pressure is given by the singular law, Eq. (5) using the generalisation of the Aronson-Bénilan estimate introduced in [19]. A minute study of the interface of the two species in two dimensions was carried out in [41].

A more complete picture is available if the cells do not avoid overcrowding due to Darcy's law but if they move according to Brinkman's law. Coupling the cell's 'velocity' to the pressure accounts for visco-elastic effects, cf. [24, 25]. A coupling through the more general Stoke's flow remains a challenging open problem. Recently, [29] proposed a two-cell-type model coupled with nutrients to study the effect of autophagy on tumour growth. In their work they, too, consider an incompressible limit, however the results remain formal due to difficulties similar to that of the system without nutrients treated by [11, 26].

## 1.3 Our Contribution

As set out in the introduction, there have been several promising steps towards establishing the incompressible limit and the complementarity relation for reaction-diffusion models incorporating convective effects. As a matter of fact, just like the authors of [40], we address the problem of passing to the incompressible limit in a porous medium equation with both a drift and a source term. While their approach is based on a viscosity solution approach, we use a weak (distributional) interpretation. By employing a blend of recently developed tools, i.e., an $L^p$-version of the celebrated Aronson-Bénilan estimate, cf. [2], along with the optimal $L^4$-regularity of the pressure gradient observed in [22], we can obtain strong compactness of the pressure gradient and proceed to passing to the incompressible limit and obtain the complementarity relation in the same vein as [11]. To summarise:

* We obtain an $L^3$-space-time estimate on the negative part of the Laplacian of the pressure which ultimately helps us obtain strong compactness of the pressure gradient. We note that an $L^\infty$-version has been obtained recently in [42, Theorem 3.1]. However, the lower bound on the Laplacian of the pressure that they infer, $\Delta p \ge -C/t - C$, does not go to zero as $\gamma \to \infty$, as in the classical Aronson-Bénilan estimate. Nonetheless, this result in conjunction with our uniform $L^4$-estimate on the pressure gradient would already be sufficient to obtain the complementarity relation rigorously, following [11, 22, 47].

* Here, we choose a different route by only striving for the much weaker $L^3$-estimate on the negative part of the Laplacian of the pressure. This, in turn, allows us to drastically
relax the $C_{x,t}^{3,1}$-regularity of the velocity field, $\nabla\Phi$, required by [42]. In fact, our assumptions on the drift, cf. Eq. (A1-$\Phi$) and Eq. (A2-$\Phi$), in a way boil down to controlling certain third derivatives in $L_{\text{loc}}^{12/5}(Q_T)$.

* Finally, to the best of our knowledge, we are the first to prove the uniqueness of the solution, $(n_{\infty}, p_{\infty})$, to the limit problem

$$ \frac{\partial n_{\infty}}{\partial t} = \Delta p_{\infty} + n_{\infty}G(p_{\infty}) + \nabla \cdot (n_{\infty}\nabla\Phi). $$

This result is only possible since we work with weak solutions in the classical sense which ultimately allows us to apply a variation of Hilbert's duality method. The only related results in this direction in the literature are given by [1] where the uniqueness of so-called patch solutions is shown in the drift-diffusion model with $\Delta\Phi > 0$ in the absence of growth dynamics and the very recent preprint [37] where uniqueness of the limit equation is shown for signed solutions, linear drifts, and general growth dynamics. In the absence of drifts uniqueness was known since [50] and for a special type of growth term it can also be obtained from $\lambda$-contractivity of metric gradient flows, cf. [16, 28].

Moreover, our approach provides an answer to several open problems proposed in [40]:

* The first question the authors raise concerns the monotonicity assumption on $G(p) + \Delta\Phi > 0$, which in our case is not necessary. An improvement in this direction has also been obtained very recently, [34]. We stress that the growth rate in [40] does not depend on the pressure but on space and time, only.

* The next question concerns the class of initial data. In [40], the authors write "A more interesting question arises with the initial data that is larger than 1 at some points. In such cases there is a jump in the solution at $t = 0$ in the limit $\gamma \to \infty$" which adds another challenge in the analysis."<fnref n="1" /> This effect has already been observed at the early stages of this singular limit problem. The parts of the density that are larger than 1 are known to "collapse" immediately and a mesa-structure is obtained instantaneously, for instance, cf. [13]. Following our approach, we can allow for the larger class of non-negative $L^1(\mathbb{R}^d) \cap L^\infty(\mathbb{R}^d)$ functions with compact support as initial data.<fnref n="2" />

* Finally, in [40], the authors postulate BV-regularity of the limiting density, also suggested by [23] based on the "five-gradients-estimate" using tools from optimal transportation. Even though our arguments do not borrow techniques from optimal transport but, instead, rely on Sobolev compactness theory, we are able to improve the regularity result in that we obtain the BV-regularity of the limit density for any initial data. What is more, we additionally have an $L^4$-regularity of the limit pressure gradient, which, to the best of our knowledge, is novel.

## 1.4 Problem Setting and Main Results

Before we present the main results of our paper let us introduce some notation used throughout this work. Henceforth, we call $Q_T := \mathbb{R}^d \times (0, T)$ the truncated space-time cylinder and drop the subscript $T$ to denote the entire cylinder, i.e., $Q := \mathbb{R}^d \times (0, \infty)$. Besides, for the sake of readability, we shall employ the short-hand notation

$$ n_{\gamma} := n_{\gamma}(t) := n_{\gamma}(x, t), $$

and, similarly,

$$ p_{\gamma} = p_{\gamma}(t) := p_{\gamma}(x, t). $$

Moreover, throughout, $C > 0$ denotes a generic positive constant independent of $\gamma$ that may change from line to line.

<footnote n="1">¹This quote is directly taken from [40] where we only adapted the notation to that of our paper.</footnote>

<footnote n="2">²While $L^\infty$-data with compact support immediately implies integrability, we trust that the assumption on the support may be removed by a localising argument in the spirit of [22, 35].</footnote>
In order to be able to establish our result we impose the following set of assumptions which, for clarity, are split into assumptions on the initial data, the growth terms, and the advective term, respectively.

We assume that for every $\gamma > 1$ the initial data are non-negative, integrable, and uniformly essentially bounded, i.e.,

$$ (A1-n_{\gamma}^{0}) \qquad n_{\gamma}^{0} \in BV(\mathbb{R}^{d}) \cap L^{\infty}(\mathbb{R}^{d}), \quad 0 \le n_{\gamma}^{0} \le n_{M}, \quad \text{and} \quad 0 \le p_{\gamma}^{0} \le p_{M}, $$

for some constants $n_M, p_M > 0$. Here $BV$ denotes the space of functions with bounded variation. Moreover, we assume the initial population is contained in a compact set, i.e., there exists a bounded set $K \subset \mathbb{R}^d$ such that

$$ (A2-n_{\gamma}^{0}) \qquad \text{supp}(n_{\gamma}^{0}) \subset K. $$

Let us notice that, thanks to the finite speed of propagation property of porous medium type equations, assumption (A2-$n_{\gamma}^{0}$) implies that, for any $T > 0$, there exists a bounded domain $\Omega \subset \mathbb{R}^d$ such that the supports of $n_{\gamma}(\cdot, t)$, $p_{\gamma}(\cdot, t)$ are contained in $\Omega$ for any $t \in [0, T]$, uniformly in $\gamma$, as proven in the next section, cf. Lemma 2.1.

In addition, we suppose that there exists a positive constant $C$ independent of $\gamma$ such that

$$ (A3-n_{\gamma}^{0}) \qquad \|\Delta(n_{\gamma}^{0})^{\gamma+1}\|_{L^{1}(\mathbb{R}^{d})} + \|\nabla p_{\gamma}^{0}\|_{L^{2}(\mathbb{R}^{d})} + \|\Delta p_{\gamma}^{0}\|_{-\|\Delta p_{\gamma}^{0}\|_{L^{2}(\mathbb{R}^{d})}} \le C. $$

Note, that strictly speaking, the $L^2$-bound on the pressure gradient is not required as it is a consequence of the $L^2$-control on the Laplacian of the pressure. Besides we make the biological assumption

$$ (A-G) \qquad G'(p) < -\alpha, \quad \text{and} \quad G(p_M) = 0, $$

for some $\alpha > 0$ and all $p \ge 0$, and some $p_M > 0$, to include the tendency of tissue to grow slower as the pressure increases and starts to die when the pressure exceeds the homeostatic pressure, $p_M$. Finally, we have to make the following regularity assumptions on the chemical distribution

$$ (A1-\Phi) \qquad \left\{ \begin{array}{l} \nabla(\partial_t \Phi) \in L^1((0, T); L_{\text{loc}}^{\infty}(\mathbb{R}^d)), \\ \Delta(\partial_t \Phi) \in L_{\text{loc}}^1(Q_T), \\ D^2\Phi \in L_{\text{loc}}^{\infty}(Q_T), \\ \nabla\Phi \in L_{\text{loc}}^2(Q_T) \cap L_{\text{loc}}^{\infty}(Q_T), \end{array} \right. $$

and

$$ (A2-\Phi) \qquad \nabla(\Delta\Phi) \in L_{\text{loc}}^{12/5}(Q_T). $$

Note, that the additional assumption, (A2-$\Phi$), is required solely for technical reasons to establish the control of the Laplacian of the pressure.

Under these hypotheses we are now able to state the two main theorems of this work. The first concerns the complementarity relation.

**Theorem 1.1 (Complementarity relation).** We may pass to the limit in Eq. (2) as $\gamma \to \infty$ and establish the so-called complementarity relation

$$ (6) \qquad p_{\infty}(\Delta p_{\infty} + \Delta\Phi + G(p_{\infty})) = 0, $$

in the distributional sense. Moreover, $0 \le n_{\infty} \le 1$ and $p_{\infty} \ge 0$ satisfy the equation

$$ (7a) \qquad \frac{\partial n_{\infty}}{\partial t} = \Delta p_{\infty} + n_{\infty}G(p_{\infty}) + \nabla \cdot (n_{\infty}\nabla\Phi), $$

in $D'(Q_T)$, as well as

$$ (7b) \qquad p_{\infty}(1 - n_{\infty}) = 0, $$

almost everywhere.
The complementarity relation, Eq. (6), is a crucial link that allows us to bridge the gap between the compressible model, Eq. (1), and the geometrical free boundary problem of Hele-Shaw type. Let us define the set

$$ \Omega(t) := \{x \mid p_{\infty}(x, t) > 0\}. $$

Then, the pressure satisfies

$$ \begin{cases} -\Delta p_{\infty} = \Delta\Phi + G(p_{\infty}), & \text{in } \Omega(t), \\ p_{\infty} = 0, & \text{on } \partial\Omega(t), \end{cases} $$

which coincides with the classical Hele-Shaw problem whenever $\Phi$ and $G$ are identically equal to zero.

**Theorem 1.2 (Uniqueness of the limit solution).** There exists at most one distributional solution such that for all $T > 0$ the couple $(n_{\infty}, p_{\infty}) \in L^{\infty}(Q_T) \times L^2(0, T; H^1(\Omega))$ is a solution to system (7a).

The rest of the paper is organised as follows. In Section 2 we present straigh-forward a priori estimates necessary to derive more refined bounds on the pressure. The latter are proven in Section 3. This includes both the $L^3$-version of the Aronson-Bénilan estimate as well as an $L^4$-space-time estimate on the pressure gradient. Building on the estimates derived in the previous sections, Section 4 is dedicated to the rigorous limit process in the pressure equation and to obtaining the complementarity relation. In the subsequent section, Section 5, we then proceed to proving the uniqueness of solutions to the complementarity relation.

## 2 A PRIORI ESTIMATES

We state some a priori estimates on the main quantities and their derivatives, that we need to obtain the main result of the paper.

**Lemma 2.1 (A priori estimates).** For any $T > 0$, there exists a bounded domain $\Omega \subset \mathbb{R}^d$ such that the supports of $n_{\gamma}(\cdot, t)$, $p_{\gamma}(\cdot, t)$ are contained in $\Omega$ for any $t \in [0, T]$, uniformly in $\gamma$. Moreover, the following estimates hold uniformly in $\gamma$:

(i) $n_{\gamma}, p_{\gamma} \in L^{\infty}(0, T; L^{\infty}(\Omega)),$

(ii) $\partial_t n_{\gamma}, \partial_t n_{\gamma} \in L^{\infty}(0, T; L^1(\Omega))$, for $i = 1, \dots, d,$

(iii) $\partial_i p_{\gamma}, \partial_t p_{\gamma} \in L^1((0, T) \times \Omega)$, for $i = 1, \dots, d,$

(iv) $\nabla p_{\gamma} \in L^2(0, T; L^2(\Omega)).$

*Proof.* Thanks to the comparison principle, from Eq. (1) we immediately find $n_{\gamma} \ge 0$ and, as a consequence, $p_{\gamma} \ge 0$. In order to establish uniform essential bounds, we construct a super solution. To this end we define

$$ \Pi(x, t) := C \left| R(t) - \frac{|x|^2}{2} \right|_{+} $$

where $C$ is a positive constant that satisfies

$$ (8) \qquad C \ge \frac{2}{d} (G(0) + \|\Delta\Phi\|_{\infty}), $$

and we take $R(t)$ such that

$$ (9) \qquad R'(t) \ge (2C + 1)R(t) + \frac{\|\nabla\Phi\|_{\infty}}{2}. $$

From Eq. (2) and the assumption on the growth term (A-G), we know that $p_{\gamma}$ satisfies

$$ \frac{\partial p_{\gamma}}{\partial t} - |\nabla p_{\gamma}|^2 - \nabla p_{\gamma} \cdot \nabla \Phi - \gamma p_{\gamma} (\Delta p_{\gamma} + G(0) + \|\Delta \Phi\|_{\infty}) \le 0. $$
Let us show that $\Pi(x, t)$ is a super-solution to this differential inequality. We have

$$ \frac{\partial \Pi}{\partial t} = C R'(t) \mathbf{1}_{\{R(t) \ge \frac{|x|^2}{2}\}}, $$

and

$$ \nabla \Pi = -C x \mathbf{1}_{\{R(t) \ge \frac{|x|^2}{2}\}}, $$

as well as

$$ \Delta \Pi = -C d \mathbf{1}_{\{R(t) \ge \frac{|x|^2}{2}\}} - C |x| \delta_{\{R(t) = \frac{|x|^2}{2}\}}. $$

Using Eq. (8) in conjunction with Eq. (9) we get

$$ \begin{align} & \frac{\partial \Pi}{\partial t} - |\nabla \Pi|^2 - \nabla \Pi \cdot \nabla \Phi - \gamma \Pi(\Delta \Pi + G(0) + \|\Delta \Phi\|_{\infty}) \nonumber \\ & \qquad \geq C R'(t) \mathbf{1}_{\{R(t) \ge \frac{|x|^2}{2}\}} - C^2 |x|^2 \mathbf{1}_{\{R(t) \ge \frac{|x|^2}{2}\}} + C x \cdot \nabla \Phi \mathbf{1}_{\{R(t) \ge \frac{|x|^2}{2}\}} + \gamma C \Pi \frac{d}{2} \nonumber \\ & \qquad \geq \left( R'(t) - 2C R(t) - \frac{|x|^2}{2} - \frac{\|\nabla \Phi\|_{\infty}}{2} \right) \mathbf{1}_{\{R(t) \ge \frac{|x|^2}{2}\}} \nonumber \\ & \qquad \geq 0. \end{align} $$

Taking $R(0)$ such that $K \subset B_{\sqrt{2R(0)}}$ and $C$ large enough, by the assumption on the initial data (A2-$n_\gamma^0$) we have $p_\gamma^0 \le \Pi(0)$. Then, this implies that $p_\gamma(t) \le \Pi(t)$ for all positive times by comparison. Let us show the argument for the sake of completeness.

Setting $N(\Pi) = \Pi^{1/\gamma}$, and multiplying Eq. (10) by $N'(\Pi)$ we obtain

$$ \frac{\partial N}{\partial t} - N'(\Pi) |\nabla \Pi|^2 - N'(\Pi) \nabla \Pi \cdot \nabla \Phi - \gamma N'(\Pi) \Pi \Delta \Pi \geq \gamma N'(\Pi) \Pi(G(0) + \|\Delta \Phi\|_{\infty}), $$

whence

$$ \frac{\partial N}{\partial t} - \nabla \cdot (N \nabla \Pi) - \nabla N \cdot \nabla \Phi \geq N(G(0) + \|\Delta \Phi\|_{\infty}). $$

Since, by Eq. (1), we know that $n_\gamma$ is a sub-solution to the same equation, we have $n_\gamma(t) \le N(t)$ for all $t > 0$, by the comparison principle. Therefore, we conclude that $p_\gamma(t) \le \Pi(t)$ for all positive times. We take $\Omega \subset \mathbb{R}^d$ a bounded domain such that $B_{\sqrt{2R(T)}} \subset \Omega$ and then, by the definition of $\Pi$, we infer that

$$ \mathrm{supp}(p_\gamma(t)) \subset \Omega, $$

for all $t \in [0, T]$ and any $\gamma > 1$. As consequence, both $n_\gamma$ and $p_\gamma$ are uniformly bounded in $L^\infty(\Omega_T)$, where $\Omega_T := \Omega \times (0, T)$.

Now we prove the BV-estimates on the density. Differentiating Eq. (1) with respect to the $i$-th component of the space variable, $x_i$, and multiplying by $\mathrm{sign}(\partial_{x_i} n_\gamma)$ we get

$$ \begin{align*} \frac{d}{dt} \int_{\Omega} \left| \frac{\partial n_{\gamma}}{\partial x_i} \right| dx &\le \int_{\Omega} \gamma \Delta \left( n_{\gamma} \left| \frac{\partial n_{\gamma}}{\partial x_i} \right| \right) dx + \int_{\Omega} \nabla \cdot \left( n_{\gamma} \nabla \left( \frac{\partial \Phi}{\partial x_i} \right) \right) \mathrm{sign} \left( \frac{\partial n_{\gamma}}{\partial x_i} \right) dx + G(0) \int_{\Omega} \left| \frac{\partial n_{\gamma}}{\partial x_i} \right| dx \\ &\le \sum_{j=1}^{d} \int_{\Omega} \left| \frac{\partial n_{\gamma}}{\partial x_j} \right| \left| \frac{\partial^2 \Phi}{\partial x_i \partial x_j} \right| dx + \sum_{j=1}^{d} \int_{\Omega} n_{\gamma} \left| \frac{\partial^3 \Phi}{\partial x_i \partial x_j^2} \right| dx + G(0) \int_{\Omega} \left| \frac{\partial n_{\gamma}}{\partial x_i} \right| dx, \end{align*} $$

for $i = 1, \dots, d$. We sum the inequalities over all $i = 1, \dots, d$, and obtain

$$ \frac{d}{dt} \sum_{i=1}^{d} \int_{\Omega} \left| \frac{\partial n_{\gamma}}{\partial x_i} \right| dx \le C \sum_{i=1}^{d} \int_{\Omega} \left| \frac{\partial n_{\gamma}}{\partial x_i} \right| dx + C, $$
where the constants depend on the $L^\infty$-norm of $G$ and the assumptions on the potential $\Phi$, cf. Eqs. (A-G, A1-Φ). Using Gronwall's lemma we conclude

$$ \sum_{i=1}^{d} \int_{\Omega} \left| \frac{\partial n_{\gamma}}{\partial x_{i}} \right| dx \leq C e^{Ct} \sum_{i=1}^{d} \int_{\Omega} \left| \frac{\partial n_{\gamma}^{0}}{\partial x_{i}} \right| dx \leq C(T), $$

where, in the last inequality, we have used the uniform BV-bounds on the initial data, cf. assumption (A1-$n_\gamma^0$).

Following the same line of reasoning for the time derivatives we obtain

$$ (11) \qquad \begin{aligned} \left| \frac{\partial}{\partial t} \left| \frac{\partial n_{\gamma}}{\partial t} \right| \right| \leq & \gamma \Delta \left( p_{\gamma} \left| \frac{\partial n_{\gamma}}{\partial t} \right| \right) + \nabla \cdot \left( \left| \frac{\partial n_{\gamma}}{\partial t} \right| \nabla \Phi \right) + \operatorname{sign} \left( \frac{\partial n_{\gamma}}{\partial t} \right) \nabla \cdot \left( n_{\gamma} \nabla \left( \frac{\partial \Phi}{\partial t} \right) \right) \\ & + \left| \frac{\partial n_{\gamma}}{\partial t} \right| G(p_{\gamma}) + n_{\gamma} G'(p_{\gamma}) \left| \frac{\partial p_{\gamma}}{\partial t} \right|, \end{aligned} $$

due to the fact that $\operatorname{sign}(\partial_t p_\gamma) = \operatorname{sign}(\partial_t n_\gamma)$. An integration in space yields

$$ \frac{d}{dt} \int_{\Omega} \left| \frac{\partial n_{\gamma}}{\partial t} \right| dx \leq G(0) \int_{\Omega} \left| \frac{\partial n_{\gamma}}{\partial t} \right| dx + \underbrace{\int_{\Omega} \left| \nabla \cdot \left( n_{\gamma} \nabla \left( \frac{\partial \Phi}{\partial t} \right) \right) \right|}_{\mathcal{I}} dx, $$

where we used that $G' < -\alpha$, cf. Eq. (A-G). We can estimate the term $\mathcal{I}$ as follows

$$ \begin{aligned} \mathcal{I} &= \int_{\Omega} \left| \nabla n_{\gamma} \cdot \nabla \left( \frac{\partial \Phi}{\partial t} \right) + n \Delta \left( \frac{\partial \Phi}{\partial t} \right) \right| dx \\ &\leq \int_{\Omega} \left| \nabla n_{\gamma} \cdot \nabla \left( \frac{\partial \Phi}{\partial t} \right) \right| dx + \int_{\Omega} \left| n \Delta \left( \frac{\partial \Phi}{\partial t} \right) \right| dx \\ &\leq \left\| \nabla \left( \frac{\partial \Phi}{\partial t} \right) (\cdot, t) \right\|_{L^{\infty}(\Omega)} \left\| \nabla n_{\gamma} \right\|_{L^{\infty}(0,T;L^1(\Omega))} + n_H \left\| \Delta \left( \frac{\partial \Phi}{\partial t} \right) (\cdot, t) \right\|_{L^1(\Omega)} \\ &\leq C \left\| \nabla \left( \frac{\partial \Phi}{\partial t} \right) (\cdot, t) \right\|_{L^{\infty}(\Omega)} + C \left\| \Delta \left( \frac{\partial \Phi}{\partial t} \right) (\cdot, t) \right\|_{L^1(\Omega)}, \end{aligned} $$

where we have used the BV-space regularity of $n_\gamma$ from before. Hence, we obtain

$$ \frac{d}{dt} \int_{\Omega} \left| \frac{\partial n_{\gamma}}{\partial t} \right| dx \leq G(0) \int_{\Omega} \left| \frac{\partial n_{\gamma}}{\partial t} \right| dx + C \left\| \nabla \left( \frac{\partial \Phi}{\partial t} \right) (\cdot, t) \right\|_{L^{\infty}(\Omega)} + C \left\| \Delta \left( \frac{\partial \Phi}{\partial t} \right) (\cdot, t) \right\|_{L^1(\Omega)}. $$

By assumption (A1-Φ) we know that $\|\nabla(\partial_t \Phi)(\cdot, t)\|_{L^\infty(\Omega)}$ and $\|\Delta(\partial_t \Phi)(\cdot, t)\|_{L^1(\Omega)}$ are $L^1$-integrable in time. Using Gronwall's lemma, we conclude

$$ (12) \qquad \begin{aligned} \left\| \frac{\partial n_{\gamma}}{\partial t}(t) \right\|_{L^1(\Omega)} &\leq e^{G(0)t} \left\| \left( \frac{\partial n_{\gamma}}{\partial t} \right)_0 \right\|_{L^1(\Omega)} \\ &\quad + \int_0^t C \left( \left\| \nabla \left( \frac{\partial \Phi}{\partial t} \right) (\cdot, t) \right\|_{L^\infty(\Omega)} + \left\| \Delta \left( \frac{\partial \Phi}{\partial t} \right) (s, \cdot) \right\|_{L^1(\Omega)} \right) e^{G(0)(t-s)} \, ds \\ &\leq C(T), \end{aligned} $$

for a.e. $t \in (0, T)$, i.e., $\partial_t n_\gamma \in L^\infty(0, T; L^1(\Omega))$. Let us stress that assumptions (A1-$n_\gamma^0$) and (A3-$n_\gamma^0$) imply the initial bound $\|(\partial_t n_\gamma)_0\|_{L^1(\Omega)} \leq C$.

Before establishing the BV-bounds on the pressure, let us notice that integrating Eq. (11) in space and time, we have

$$ \left\| \frac{\partial n_{\gamma}}{\partial t} (\cdot, t) \right\|_{L^1(\Omega)} + \min_{0 \leq p_{\gamma} \leq \Pi(0, T)} \left| G'(p_{\gamma}) \right| \int_0^t \int_{\Omega} n_{\gamma} \left| \frac{\partial p_{\gamma}}{\partial t} \right| dx dt \leq C(T), $$
thanks to Eq. (12). Then, it holds

$$ \left\| \frac{\partial p_{\gamma}}{\partial t} \right\|_{L^1(\Omega_T)} \le \iint_{\Omega_T \cap \{n_{\gamma} \le 1/2\}} \gamma n_{\gamma}^{\gamma-1} \left| \frac{\partial n_{\gamma}}{\partial t} \right| dx dt + 2 \iint_{\Omega_T \cap \{n_{\gamma} > 1/2\}} n_{\gamma} \left| \frac{\partial p_{\gamma}}{\partial t} \right| dx dt \le C(T). $$

The same argument can be used for the space derivatives of $p_γ$ without major changes. We can actually gain more information on the pressure gradient, by integrating Eq. (2) in space, i.e.,

$$ \int_{\Omega} \frac{\partial p_{\gamma}}{\partial t} dx = \gamma \int_{\Omega} p_{\gamma}(\Delta(p_{\gamma} + \Phi) + G(p_{\gamma})) dx + \int_{\Omega} \nabla p_{\gamma} \cdot \nabla(p_{\gamma} + \Phi) dx. $$

Integration by parts yields

$$ \int_{\Omega} \frac{\partial p_{\gamma}}{\partial t} dx = (1 - \gamma) \int_{\Omega} |\nabla p_{\gamma}|^2 dx + \gamma \int_{\Omega} p_{\gamma} G(p_{\gamma}) dx + (1 - \gamma) \int_{\Omega} \nabla p_{\gamma} \cdot \nabla \Phi dx, $$

and using Young's inequality we obtain

$$ \frac{\gamma - 1}{2} \iint_{\Omega_T} |\nabla p_{\gamma}(t)|^2 dx dt \le \|p_{\gamma}^0\|_{L^1(\Omega)} + \frac{(\gamma - 1)}{2} \iint_{\Omega_T} |\nabla \Phi|^2 dx dt + \gamma \iint_{\Omega_T} |p_{\gamma} G(p_{\gamma})| dx dt. $$

Dividing by $(\gamma - 1)$ we finally get

$$ \iint_{\Omega_T} |\nabla p_{\gamma}|^2 dx dt \le C(T), $$

which concludes the proof. □

### 3 STRONGER BOUNDS ON $p_γ$

This section is dedicated to establishing more refined estimates on the pressure, cf. Lemma 3.2 and Lemma 3.3. Upon obtaining those estimates we will then be able to proceed to proving the strong compactness of the pressure gradient, cf. Lemma 3.6, which is crucial in the overall endeavour of establishing the incompressible limit.

The first result on the pressure's regularity is the $L^4$-boundedness of its gradient. This bound was already proved in [47], although, the authors use the $L^\infty$-version of the Aronson-Bénilan estimate. Here we adapted the method used in [22], where a new method was employed, that does not require any estimate on $\Delta p_γ$. Unlike the model in [22], the convective term may not vanish at the boundary which leads to boundary terms to be considered in the subsequent analysis. In the following remark we shall see, however, that they do not pose any problems.

**Remark 3.1** (Boundary Terms and Integration by Parts). The subsequent technical lemmas (Lemma 3.2 and Lemma 3.3) are critical to establishing the regularity necessary for passing to the stiff limit. Due to several integrations by parts, boundary terms occur that need to be addressed. Since their treatment is purely technical and they are not even at the heart of the strategy we introduce the notation $O_{\partial\Omega_T}(1)$ to indicate that the traces of the respective quantities are bounded uniformly in $\gamma$. This is possible due to the elliptic regularity result presented in [33, Theorem 9.11] which states that

$$ \|u\|_{H^2(U')} \le C(\|u\|_{L^2(U)} + \|\Delta u\|_{L^2(U)}), $$

for some open $U \subset \mathbb{R}^n$ containing $U' \subset \text{compactly}$. Choosing $u = \partial_t\Phi$, for all $i = 1, \dots, d$, and using assumption (A2-Φ), it is immediate that $\nabla\Delta\Phi \in H^2(Q_T)$. With the third-order derivatives controlled in $L^2(Q_T)$ the traces of all second order derivatives appearing in the integration by parts are bounded. Let us highlight, too, that terms involving $p_γ$ and its derivatives vanish close to the boundary by the choice of $\Omega_T$. We therefore collect all boundary terms in $O_{\partial\Omega_T}(1)$ lest the notation blow up.
**Lemma 3.2** (L⁴-estimate of the pressure gradient.). Given T > 0, there exists a positive constant C, independent of γ, such that

$$\iint_{\Omega_T} p_{\gamma} \sum_{i,j=1}^{d} \left| \frac{\partial^2 p_{\gamma}}{\partial x_i \partial x_j} \right|^2 dx dt + (\gamma - 1) \iint_{\Omega_T} p_{\gamma} |\Delta p_{\gamma} + \Delta \Phi + G|^2 dx dt \leq C(T),$$

as well as

$$\iint_{\Omega_T} |\nabla p_{\gamma}|^4 dx dt \leq C(T).$$

*Proof.* We write the equation for the pressure as follows

$$(13) \qquad \frac{\partial p_{\gamma}}{\partial t} = \gamma p_{\gamma}(\Delta f_{\gamma} + G) + \nabla p_{\gamma} \cdot \nabla f_{\gamma},$$

where $f_{\gamma} := p_{\gamma} + \Phi$. We multiply Eq. (13) by $-(\Delta f_{\gamma} + G)$ and integrate in space and time to obtain

$$(14) \qquad \int_0^T \frac{d}{dt} \int_{\Omega} \frac{|\nabla p_{\gamma}|^2}{2} dx dt - \iint_{\Omega_T} \Delta \Phi \frac{\partial p_{\gamma}}{\partial t} dx dt - \iint_{\Omega_T} G \frac{\partial p_{\gamma}}{\partial t} dx dt \\ = - \underbrace{\iint_{\Omega_T} \nabla p_{\gamma} \cdot \nabla f_{\gamma} (\Delta f_{\gamma} + G)}_{\mathcal{I}} dx dt - \gamma \iint_{\Omega_T} p_{\gamma} |\Delta f_{\gamma} + G|^2 dx dt.$$

For convenience, let us define the function $\bar{G} = \bar{G}(p_{\gamma}) = \int_{0}^{p_{\gamma}} G(q) dq$. Thus, we have

$$\partial_t p_{\gamma} G(p_{\gamma}) = \partial_t \bar{G}(p_{\gamma}),$$

and thus

$$\iint_{\Omega_T} \frac{\partial p_{\gamma}}{\partial t} G(p_{\gamma}) dx dt = \int_0^T \frac{d}{dt} \int_{\Omega} \bar{G}(p_{\gamma}) dx dt.$$

Now, we need to estimate the term $\mathcal{I}$ on the right-hand side of Eq. (14). Since $p_{\gamma} = f_{\gamma} - \Phi$ we have

$$\begin{aligned} \mathcal{I} &= - \iint_{\Omega_T} \nabla p_{\gamma} \cdot \nabla f_{\gamma} (\Delta f_{\gamma} + G) dx dt \\ &= - \iint_{\Omega_T} |\nabla f_{\gamma}|^2 \Delta f_{\gamma} dx dt + \iint_{\Omega_T} \nabla \Phi \cdot \nabla f_{\gamma} \Delta f_{\gamma} dx dt - \iint_{\Omega_T} G \nabla p_{\gamma} \cdot \nabla f_{\gamma} dx dt \\ &\leq - \underbrace{\iint_{\Omega_T} |\nabla f_{\gamma}|^2 \Delta f_{\gamma} dx dt}_{\mathcal{I}_1} + \underbrace{\iint_{\Omega_T} \nabla \Phi \cdot \nabla f_{\gamma} \Delta f_{\gamma} dx dt}_{\mathcal{I}_2} + C, \end{aligned}$$

thanks to the L²-bounds of both $\nabla p_{\gamma}$ and $\nabla \Phi$. We integrate by parts twice in space the term $\mathcal{I}_1$ and obtain

$$\begin{aligned} \mathcal{I}_1 &= \iint_{\Omega_T} f_{\gamma} \Delta (|\nabla f_{\gamma}|^2) dx dt \\ &= 2 \iint_{\Omega_T} f_{\gamma} \nabla f_{\gamma} \cdot \nabla (\Delta f_{\gamma}) dx dt + 2 \iint_{\Omega_T} f_{\gamma} \sum_{i,j=1}^{d} \left| \frac{\partial^2 f}{\partial x_i \partial x_j} \right|^2 dx dt + \mathcal{O}_{\partial \Omega_T}(1) \\ &= -2 \iint_{\Omega_T} f_{\gamma} |\Delta f_{\gamma}|^2 dx dt - 2 \iint_{\Omega_T} |\nabla f_{\gamma}|^2 \Delta f_{\gamma} dx dt + 2 \iint_{\Omega_T} f_{\gamma} \sum_{i,j=1}^{d} \left| \frac{\partial^2 f}{\partial x_i \partial x_j} \right|^2 dx dt + \mathcal{O}_{\partial \Omega_T}(1). \end{aligned}$$
Let us notice that the second term on the right-hand side is equal to $-2\mathcal{I}_1$. Hence, moving it to the left-hand side of the equation and simplifying the expression we obtain

$$ 
\begin{aligned}
-\mathcal{I}_1 &= -\iint_{\Omega_T} |\nabla f_\gamma|^2 \Delta f_\gamma \,dx\,dt \\
&= \frac{2}{3} \iint_{\Omega_T} f_\gamma |\Delta f_\gamma|^2 \,dx\,dt - \frac{2}{3} \iint_{\Omega_T} f_\gamma \sum_{i,j=1}^d \left| \frac{\partial^2 f_\gamma}{\partial x_i \partial x_j} \right|^2 \,dx\,dt + \mathcal{O}_{\partial\Omega_T}(1) \\
&= \frac{2}{3} \iint_{\Omega_T} p_\gamma |\Delta f_\gamma|^2 \,dx\,dt - \frac{2}{3} \iint_{\Omega_T} p_\gamma \sum_{i,j=1}^d \left| \frac{\partial^2 f_\gamma}{\partial x_i \partial x_j} \right|^2 \,dx\,dt \\
&\quad + \frac{2}{3} \iint_{\Omega_T} \Phi |\Delta f_\gamma|^2 \,dx\,dt - \frac{2}{3} \iint_{\Omega_T} \Phi \sum_{i,j=1}^d \left| \frac{\partial^2 f_\gamma}{\partial x_i \partial x_j} \right|^2 \,dx\,dt + \mathcal{O}_{\partial\Omega_T}(1).
\end{aligned}
$$

We now compute the sum of the last two integrals of the right-hand side

$$ 
\begin{aligned}
& \frac{2}{3} \iint_{\Omega_T} \Phi |\Delta f_\gamma|^2 \,dx\,dt - \frac{2}{3} \iint_{\Omega_T} \Phi \sum_{i,j=1}^d \left| \frac{\partial^2 f_\gamma}{\partial x_i \partial x_j} \right|^2 \,dx\,dt \\
&= \frac{2}{3} \iint_{\Omega_T} \left( \sum_{i,j=1}^d \frac{\partial f_\gamma}{\partial x_j} \frac{\partial^2 \Phi}{\partial x_i \partial x_j} \frac{\partial f_\gamma}{\partial x_i} \,dx\,dt - \Delta \Phi |\nabla f_\gamma|^2 \right) \,dx\,dt \\
&\le C (\|D^2 \Phi\|_{L^\infty} \|\nabla f_\gamma\|_{L^2}^2 + \|\Delta \Phi\|_{L^\infty} \|\nabla f_\gamma\|_{L^2}^2) \\
&\le C,
\end{aligned}
$$

having used the assumptions on the velocity field, cf. (A1-Φ), and the information on the pressure gradient, cf. Lemma 2.1. Therefore, we can estimate the term $-\mathcal{I}_1$ as follows

$$ -\mathcal{I}_1 \le \frac{2}{3} \iint_{\Omega_T} p_\gamma |\Delta f_\gamma|^2 \,dx\,dt - \frac{2}{3} \iint_{\Omega_T} p_\gamma \sum_{i,j=1}^d \left| \frac{\partial^2 f_\gamma}{\partial x_i \partial x_j} \right|^2 \,dx\,dt + C. $$

Now we proceed integrating by parts and estimating the term $\mathcal{I}_2$

$$ 
\begin{aligned}
\mathcal{I}_2 &= \iint_{\Omega_T} \nabla \Phi \cdot \nabla f_\gamma \Delta f_\gamma \,dx\,dt \\
&= - \iint_{\Omega_T} \sum_{i,j=1}^d \frac{\partial f_\gamma}{\partial x_j} \frac{\partial^2 \Phi}{\partial x_i \partial x_j} \frac{\partial f_\gamma}{\partial x_i} \,dx\,dt - \iint_{\Omega_T} \sum_{i,j=1}^d \frac{\partial \Phi}{\partial x_j} \frac{\partial^2 f_\gamma}{\partial x_i \partial x_j} \frac{\partial f_\gamma}{\partial x_i} \,dx\,dt + \mathcal{O}_{\partial\Omega_T}(1) \\
&\le C \|D^2 \Phi\|_{L^\infty} \|\nabla f_\gamma\|_{L^2}^2 - \iint_{\Omega_T} \sum_{i,j=1}^d \frac{\partial \Phi}{\partial x_j} \frac{\partial^2 f_\gamma}{\partial x_i \partial x_j} \frac{\partial f_\gamma}{\partial x_i} \,dx\,dt + \mathcal{O}_{\partial\Omega_T}(1) \\
&\le C - \frac{1}{2} \iint_{\Omega_T} \nabla \Phi \cdot \nabla |\nabla f_\gamma|^2 \,dx\,dt + \mathcal{O}_{\partial\Omega_T}(1) \\
&= C + \frac{1}{2} \iint_{\Omega_T} \Delta \Phi \cdot |\nabla f_\gamma|^2 \,dx\,dt + \mathcal{O}_{\partial\Omega_T}(1) \\
&\le C + \frac{1}{2} \|\Delta \Phi\|_{L^\infty} \|\nabla f_\gamma\|_{L^2}^2 + \mathcal{O}_{\partial\Omega_T}(1) \\
&\le C.
\end{aligned}
$$
Therefore, we obtain

$$
\begin{align*}
\mathcal{I} &\leq -\mathcal{I}_1 + \mathcal{I}_2 \\
&\leq \frac{2}{3} \iint_{\Omega_T} p_{\gamma} |\Delta f_{\gamma}|^2 \, dx \, dt - \frac{2}{3} \iint_{\Omega_T} p_{\gamma} \sum_{i,j=1}^{d} \left| \frac{\partial^2 f_{\gamma}}{\partial x_i \partial x_j} \right|^2 \, dx \, dt + C \\
&\leq \frac{2}{3} \iint_{\Omega_T} p_{\gamma} |\Delta f_{\gamma} + G|^2 \, dx \, dt - \frac{2}{3} \iint_{\Omega_T} p_{\gamma} \sum_{i,j=1}^{d} \left| \frac{\partial^2 f_{\gamma}}{\partial x_i \partial x_j} \right|^2 \, dx \, dt + C,
\end{align*}
$$

where in the last inequality we used the fact that G is uniformly bounded.
Gathering all the bounds we can write Eq. (14) as

$$
\begin{align*}
& \frac{2}{3} \iint_{\Omega_T} p_{\gamma} \sum_{i,j=1}^{d} \left| \frac{\partial^2 f_{\gamma}}{\partial x_i \partial x_j} \right|^2 dx \, dt + \left( \gamma - \frac{2}{3} \right) \iint_{\Omega} p_{\gamma} |\Delta f_{\gamma} + G|^2 dx \, dt \\
& \leq \int_{0}^{T} \frac{d}{dt} \int_{\Omega} \left( \bar{G} - \frac{|\nabla p_{\gamma}|^2}{2} \right) dx \, dt + \iint_{\Omega_T} \Delta \Phi \frac{\partial p_{\gamma}}{\partial t} dx \, dt + C \\
& \leq C(T),
\end{align*}
$$

where in the last inequality we used the L¹-bound of ∂_t p_γ. Thus, we have proved the following
bound

$$
\frac{2}{3} \iint_{\Omega_T} p_{\gamma} \sum_{i,j=1}^{d} \left| \frac{\partial^2 f_{\gamma}}{\partial x_i \partial x_j} \right|^2 dx \, dt + \left( \gamma - \frac{2}{3} \right) \iint_{\Omega_T} p_{\gamma} |\Delta f_{\gamma} + G|^2 dx \, dt \leq C(T),
$$

Finally, thanks to the boundedness of $\partial_{i,j}^2 \Phi$, we have

$$
\begin{equation} \tag{15}
\begin{split}
& \iint_{\Omega_T} p_{\gamma} \sum_{i,j=1}^{d} \left| \frac{\partial^2 p_{\gamma}}{\partial x_i \partial x_j} \right|^2 dx \, dt \\
& \le 2 \iint_{\Omega_T} p_{\gamma} \sum_{i,j=1}^{d} \left| \frac{\partial^2 f_{\gamma}}{\partial x_i \partial x_j} \right|^2 dx \, dt + 2 \iint_{\Omega_T} p_{\gamma} \sum_{i,j=1}^{d} \left| \frac{\partial^2 \Phi}{\partial x_i \partial x_j} \right|^2 dx \, dt \\
& \le C(T),
\end{split}
\end{equation}
$$

and since $\gamma > 1$

$$
\begin{equation} \tag{16}
\begin{split}
\iint_{\Omega_T} p_{\gamma} |\Delta p_{\gamma}|^2 dx \, dt &\le 2 \iint_{\Omega_T} p_{\gamma} |\Delta f_{\gamma} + G|^2 dx \, dt + 2 \iint_{\Omega_T} p_{\gamma} |\Delta \Phi + G|^2 dx \, dt \\
&\le C(T),
\end{split}
\end{equation}
$$

and the first part of the lemma is proven. Now it remains to prove the L⁴-bound of the pressure
gradient. Integrating by parts we have

$$
\int_{\Omega} |\nabla p_{\gamma}|^4 dx = - \int_{\Omega} p_{\gamma} \Delta p_{\gamma} |\nabla p_{\gamma}|^2 dx - \int_{\Omega} p_{\gamma} \nabla p_{\gamma} \cdot \nabla (|\nabla p_{\gamma}|^2) dx.
$$

Applying Young's inequality to the first term, we obtain

$$
\frac{1}{2} \int_{\Omega} |\nabla p_{\gamma}|^4 dx \le \frac{1}{2} \int_{\Omega} p_{\gamma}^2 |\Delta p_{\gamma}|^2 dx - 2 \sum_{i,j=1}^{d} \int_{\Omega} p_{\gamma} \frac{\partial p_{\gamma}}{\partial x_i} \frac{\partial p_{\gamma}}{\partial x_j} \frac{\partial^2 p_{\gamma}}{\partial x_i \partial x_j} dx.
$$

Thanks to Young's inequality, the last term can be bounded from above by

$$
\left| 2 \sum_{i,j=1}^{d} \int_{\Omega} p_{\gamma} \frac{\partial p_{\gamma}}{\partial x_i} \frac{\partial p_{\gamma}}{\partial x_j} \frac{\partial^2 p_{\gamma}}{\partial x_i \partial x_j} dx \right| \le \frac{1}{4} \int_{\Omega} |\nabla p_{\gamma}|^4 dx + 4 \int_{\Omega} p_{\gamma}^2 \sum_{i,j=1}^{d} \left| \frac{\partial^2 p_{\gamma}}{\partial x_i \partial x_j} \right|^2 dx.
$$
Therefore, we obtain

$$ \frac{1}{4} \int_{\Omega} |\nabla p_{\gamma}|^4 dx \leq \frac{1}{2} \int_{\Omega} p_{\gamma}^2 |\Delta p_{\gamma}|^2 dx + 4 \int_{\Omega} p_{\gamma}^2 \sum_{i,j=1}^{d} \left| \frac{\partial^2 p_{\gamma}}{\partial x_i \partial x_j} \right|^2 dx. $$

Since $p_{\gamma} \leq \Pi(0, T)$ and thanks to Eqs. (15,16), we conclude that

$$ \iint_{\Omega_T} |\nabla p_{\gamma}|^4 dx dt \leq C(T), $$

which completes the proof.

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

Building on the $L^4$-estimate on the pressure gradient, we are now dedicated to an additional bound on the pressure which, by itself, yields $L^1$-compactness of the pressure gradient. In conjunction with the $L^4$-estimate the gradient is then shown to be strongly compact in any $L^p(\Omega_T)$, for $1 \leq p < 4$, cf. Lemma 3.6. The subsequent estimate is an $L^p$-version of the celebrated Aronson-Bénilan estimate, cf. [2, 7]. At the heart of its proof is the study of an auxiliary second-order quantity and its evolution along the flow of the pressure equation. We define $w := \Delta p_{\gamma} + G(p_{\gamma})$ and, for the reader's convenience, recall that the pressure satisfies the equation

$$ (17) \qquad \frac{\partial p_{\gamma}}{\partial t} = \gamma p_{\gamma} w + \gamma p_{\gamma} \Delta \Phi + \nabla p_{\gamma} \cdot (\nabla p_{\gamma} + \nabla \Phi). $$

**Lemma 3.3** (Aronson-Bénilan $L^3$-estimate.). For all $T > 0$ and $\gamma > \max(1, 2 - \frac{2}{d})$, there exists a positive constant $C(T)$, independent of $\gamma$, such that

$$ \iint_{\Omega_T} |w|^3 \underline{dx} dt \leq C(T). $$

*Proof.* We compute the time derivative of $w$

$$ \begin{aligned} \frac{\partial w}{\partial t} = & \gamma \Delta(p_{\gamma}w) + \gamma p_{\gamma} \Delta(\Delta\Phi) + \gamma(w-G)\Delta\Phi + 2\gamma \nabla p_{\gamma} \cdot \nabla(\Delta\Phi) + 2\nabla p_{\gamma} \cdot \nabla(w-G) \\ & + 2 \sum_{i,j=1}^{d} \left| \frac{\partial^2 p_{\gamma}}{\partial x_i \partial x_j} \right|^2 + \nabla(w-G) \cdot \nabla\Phi + \nabla p_{\gamma} \cdot \nabla(\Delta\Phi) + 2 \sum_{i,j=1}^{d} \frac{\partial^2 p_{\gamma}}{\partial x_i \partial x_j} \frac{\partial^2 \Phi}{\partial x_i \partial x_j} + G' \frac{\partial p_{\gamma}}{\partial t}. \end{aligned} $$

Young's inequality yields

$$ \left| 2 \sum_{i,j=1}^{d} \partial_{i,j}^2 p_{\gamma} \frac{\partial^2 \Phi}{\partial x_i \partial x_j} \right| \leq \sum_{i,j=1}^{d} \left| \frac{\partial^2 p_{\gamma}}{\partial x_i \partial x_j} \right|^2 + \sum_{i,j=1}^{d} \left| \frac{\partial^2 \Phi}{\partial x_i \partial x_j} \right|^2, $$

and thus, using Eq. (17), we get

$$ \begin{aligned} \frac{\partial w}{\partial t} \geq & \gamma \Delta(p_{\gamma}w) + \gamma p_{\gamma} \Delta(\Delta\Phi) + \gamma w \Delta\Phi - \gamma G \Delta\Phi + (2\gamma + 1) \nabla p_{\gamma} \cdot \nabla(\Delta\Phi) + 2 \nabla p_{\gamma} \cdot \nabla w \\ & - 2 |\nabla p|^2 G' + \sum_{i,j=1}^{d} \left| \frac{\partial^2 p_{\gamma}}{\partial x_i \partial x_j} \right|^2 - \sum_{i,j=1}^{d} \left| \frac{\partial^2 \Phi}{\partial x_i \partial x_j} \right|^2 + \nabla w \cdot \nabla \Phi - G' \nabla p \cdot \nabla \Phi \\ & + \gamma G' p_{\gamma} w + \gamma p_{\gamma} G' \Delta \Phi + G' |\nabla p_{\gamma}|^2 + G' \nabla p_{\gamma} \cdot \nabla \Phi. \end{aligned} $$

We use the fact that

$$ \sum_{i,j=1}^{d} \left| \frac{\partial^2 p_{\gamma}}{\partial x_i \partial x_j} \right|^2 \geq \frac{1}{d} |\Delta p_{\gamma}|^2 = \frac{1}{d} (w - G)^2, $$
and we obtain

$$ \frac{\partial w}{\partial t} \geq \gamma \Delta(p_{\gamma} w) + \gamma p_{\gamma} \Delta(\Delta\Phi) + \gamma w \Delta\Phi - \gamma G \Delta\Phi + (2\gamma + 1) \nabla p_{\gamma} \cdot \nabla(\Delta\Phi) + 2 \nabla p_{\gamma} \cdot \nabla w \\ - |\nabla p|^2 G' + \frac{1}{d} w^2 - \frac{2}{d} w G + \frac{1}{d} G^2 - \sum_{i,j=1}^{d} \left| \frac{\partial^2 \Phi}{\partial x_i \partial x_j} \right|^2 + \nabla w \cdot \nabla \Phi \\ + \gamma G' p_{\gamma} w + \gamma p_{\gamma} G' \Delta\Phi. $$

We multiply by $-|w|_-$, to find

$$ -\frac{\partial w}{\partial t}|_{w=-} \leq -\frac{1}{d}|w|_{-}^3 + \gamma \Delta\Phi|w|_{-}^2 - \frac{2}{d}G|w|_{-}^2 + \gamma G' p_{\gamma}|w|_{-}^2 - \frac{1}{d}G^2|w|_{-} + \gamma G \Delta\Phi|w|_{-} \\ + \sum_{i,j=1}^{d} \left| \frac{\partial^2 \Phi}{\partial x_i \partial x_j} \right|^2 |w|_{-} - \gamma p_{\gamma} G' \Delta\Phi|w|_{-} + |\nabla p_{\gamma}|^2 G' |w|_{-} \\ + \gamma \Delta(p_{\gamma}|w|_{-})|w|_{-} + 2 \nabla p_{\gamma} \cdot \nabla |w|_{-} |w|_{-} \\ - \gamma p_{\gamma} \Delta(\Delta\Phi)|w|_{-} - (2\gamma + 1) \nabla p_{\gamma} \cdot \nabla(\Delta\Phi)|w|_{-} \\ + \nabla \Phi \cdot \nabla |w|_{-} |w|_{-}. $$

Hence, using the fact that $G' < -\alpha$ and integrating in space and time, we obtain

$$ \begin{aligned} - \int_{\Omega} \frac{|w^0|_-^2}{2} dx \leq & - \frac{1}{d} \iint_{\Omega_T} |w|_-^3 dx dt + C\gamma \iint_{\Omega_T} |w|_-^2 dx dt + C\gamma \iint_{\Omega_T} |w|_- dx dt \\ & + \underbrace{\gamma \iint_{\Omega_T} \Delta(p_{\gamma}|w|_-)|w|_- + 2\nabla p_{\gamma} \cdot \nabla |w|_- |w|_-}_{\mathcal{I}_1} dx dt \end{aligned} \tag{18} $$ $$ \begin{aligned} & -\gamma \underbrace{\iint_{\Omega_T} p_{\gamma} \Delta(\Delta\Phi)|w|_- dx dt}_{\mathcal{I}_2} - (2\gamma + 1) \underbrace{\iint_{\Omega_T} \nabla p_{\gamma} \cdot \nabla(\Delta\Phi)|w|_- dx dt}_{\mathcal{I}_3} \\ & + \underbrace{\iint_{\Omega_T} \nabla\Phi \cdot \nabla |w|_- |w|_- dx dt}_{\mathcal{I}_4} \end{aligned} $$

where $C$ represents different constants depending on the $L^\infty$-norms of $G, G'$ and $\partial_{i,j}^2 \Phi$, for $i, j = 1, \dots, d$.

Now, we compute each term individually. Integration by parts yields

$$ \begin{aligned} \mathcal{I}_1 &= \gamma \iint_{\Omega_T} \Delta(p_{\gamma}|w|_-)|w|_- + 2\nabla p_{\gamma} \cdot \nabla |w|_- |w|_- dx dt \\ &= -\frac{\gamma}{2} \iint_{\Omega_T} \nabla p_{\gamma} \cdot \nabla |w|_-^2 dx dt - \gamma \iint_{\Omega_T} p_{\gamma} |\nabla |w|_-|^2 dx dt + \iint_{\Omega_T} \nabla p_{\gamma} \cdot \nabla |w|_-^2 dx dt \\ &= -\left(1 - \frac{\gamma}{2}\right) \iint_{\Omega_T} (w - G)|w|_-^2 dx dt - \gamma \iint_{\Omega_T} p_{\gamma} |\nabla |w|_-|^2 dx dt \\ &= \left(1 - \frac{\gamma}{2}\right) \iint_{\Omega_T} |w|_-^3 dx dt + \left(1 - \frac{\gamma}{2}\right) \iint_{\Omega_T} G|w|_-^2 dx dt - \gamma \iint_{\Omega_T} p_{\gamma} |\nabla |w|_-|^2 dx dt \\ &\leq \left(1 - \frac{\gamma}{2}\right) \iint_{\Omega_T} |w|_-^3 dx dt - \gamma \iint_{\Omega_T} p_{\gamma} |\nabla |w|_-|^2 dx dt + C\gamma \iint_{\Omega_T} |w|_-^2 dx dt. \end{aligned} $$
We continue by using integration by parts and Young's inequality to get

$$
\begin{align*}
\mathcal{I}_2 &= -\gamma \iint_{\Omega_T} p_\gamma \Delta(\Delta\Phi) |w|_- \,dx\,dt \\
&= \gamma \iint_{\Omega_T} p_\gamma \nabla(\Delta\Phi) \cdot \nabla|w|_- \,dx\,dt + \gamma \iint_{\Omega_T} \nabla p_\gamma \cdot \nabla(\Delta\Phi) |w|_- \,dx\,dt \\
&\le \frac{\gamma}{2} \iint_{\Omega_T} p_\gamma |\nabla|w|_-|^2 \,dx\,dt + \frac{\gamma}{2} \iint_{\Omega_T} p_\gamma |\nabla(\Delta\Phi)|^2 \,dx\,dt \\
&\quad + \gamma \left( \iint_{\Omega_T} |\nabla p_\gamma|^4 \right)^{1/4} \left( \iint_{\Omega_T} |\nabla(\Delta\Phi)| |w|_- |^4 |^3 \,dx\,dt \right)^{3/4} \\
&\le \frac{\gamma}{2} \iint_{\Omega_T} p_\gamma |\nabla|w|_-|^2 \,dx\,dt + \frac{\gamma}{2} \iint_{\Omega_T} p_\gamma |\nabla(\Delta\Phi)|^2 \,dx\,dt \\
&\quad + C\gamma \left( \iint_{\Omega_T} |\nabla(\Delta\Phi)|^{12/5} \,dx\,dt \right)^{5/12} \left( \iint_{\Omega_T} |w|_-^3 \,dx\,dt \right)^{1/3} \\
&\le \frac{\gamma}{2} \iint_{\Omega_T} p_\gamma |\nabla|w|_-|^2 \,dx\,dt + C\gamma + C\gamma \left( \iint_{\Omega_T} |w|_-^3 \,dx\,dt \right)^{1/3},
\end{align*}
$$

where we used Hölder's inequality, the $L^4$-bound of the pressure gradient of Lemma 3.2 and the assumption (A2-Φ), $\nabla(\Delta\Phi) \in L_{\text{loc}}^{12/5}(Q_T)$.

Using again Young's and Holder's inequalities we have

$$
\begin{align*}
\mathcal{I}_3 &\le (2\gamma + 1) \left( \iint_{\Omega_T} |\nabla p_\gamma|^4 \,dx\,dt \right)^{1/4} \left( \iint_{\Omega_T} |\nabla(\Delta\Phi)| |w|_- |^4 |^3 \,dx\,dt \right)^{3/4} \\
&\le C\gamma \left( \iint_{\Omega_T} |\nabla(\Delta\Phi)|^{12/5} \,dx\,dt \right)^{5/12} \left( \iint_{\Omega_T} |w|_-^3 \,dx\,dt \right)^{1/3} \\
&\le C\gamma \left( \iint_{\Omega_T} |w|_-^3 \,dx\,dt \right)^{1/3}.
\end{align*}
$$

The last term is

$$
\mathcal{I}_4 = \iint_{\Omega_T} \frac{1}{2} \nabla \Phi \cdot \nabla |w|_-^2 \,dx\,dt = -\frac{1}{2} \iint_{\Omega_T} \Delta \Phi |w|_-^2 \,dx\,dt \le C \iint_{\Omega_T} |w|_-^2 \,dx\,dt.
$$

Here we have used the fact that $\Omega$ is a compact set which contains $\text{supp}(p_\gamma)$ and large enough
such that $\Delta p_\gamma = 0$ on $\partial\Omega$, then $|w|_-=0$ on $\partial\Omega$.

Hence, gathering all the estimates and using Hölder's inequality, we can rewrite Eq. (18) as

$$
\left(\frac{\gamma}{2} - 1 + \frac{1}{d}\right) \iint_{\Omega_T} |w|_-^3 \,dx\,dt \le C\gamma \left(\iint_{\Omega_T} |w|_-^3 \,dx\,dt\right)^{1/3} + C\gamma \left(\iint_{\Omega_T} |w|_-^3 \,dx\,dt\right)^{2/3} + C\gamma,
$$

since we assumed $|w^0|_- \in L^2(\mathbb{R}^d)$. Finally, for $\gamma > \max(1, 2 - 2/d)$, we have

$$
\iint_{\Omega_T} |w|_-^3 \,dx\,dt \le C \left( \iint_{\Omega_T} |w|_-^3 \,dx\,dt \right)^{1/3} + C \left( \iint_{\Omega_T} |w|_-^3 \,dx\,dt \right)^{2/3} + C,
$$

which yields

$$
\iint_{\Omega_T} |w|_-^3 \,dx\,dt \le C(T),
$$

where $C(T)$ depends on $T$, $|\Omega|$ and previous uniform bounds, and the proof is concluded.

□
**Corollary 3.4.** It holds

$$ (19) \qquad \iint_{\Omega_T} |\Delta p_\gamma| \, dx \, dt \le C(T). $$

*Proof.* The compact support assumption yields

$$ \iint_{\Omega_T} (\Delta p_\gamma + G) \, dx \, dt \le C(T), $$

and then, thanks to Hölder's inequality, we have

$$ \begin{aligned} \iint_{\Omega_T} |\Delta p_\gamma + G| \, dx \, dt &= \iint_{\Omega_T} (\Delta p_\gamma + G) \, dx \, dt + 2 \iint_{\Omega_T} |w|^{-1} \, dx \, dt \\ &\le C(T) + C \left( \iint_{\Omega_T} |w|^{3/2} \, dx \, dt \right)^{1/3} \\ &\le C(T). \end{aligned} $$

Finally, since $G$ is bounded, we obtain

$$ \iint_{\Omega_T} |\Delta p_\gamma| \, dx \, dt \le C(T). $$

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

**Remark 3.5.** The proof of the Aronson-Bénilan estimate can be made independent of the $L^4$-bound on $\nabla p_\gamma$ imposing a stronger condition on $\Phi$, namely $\nabla(\Delta\Phi) \in L^6$ rather than $L^{12/5}$.

The bounds provided by Lemma 3.2 and Lemma 3.3 allow us to prove the strong convergence of $\nabla p_\gamma$ in $L^2(Q_T)$ thanks to compactness arguments, in particular the Fréchet-Kolmogorov theorem and the Aubin-Lions lemma.

**Lemma 3.6** (Strong convergence of the pressure gradient). For any $T > 0$ it holds

$$ \nabla p_\gamma \to \nabla p_\infty, $$

strongly in $L^2(Q_T)$.

*Proof.* Thanks to Lemma 3.2, we infer the weak convergence (up to a subsequence) of the pressure gradient

$$ (20) \qquad \nabla p_\gamma \to \nabla p_\infty, $$

weakly in $L^4(Q_T)$. From Lemma 3.3, we know that $\Delta p_\gamma$ is bounded in $L^1(Q_T)$, which is instrumental in establishing space-time compactness in any $L^r(Q_T)$, with $1 \le r < 4$. The proof of this claim is an extension of [8, Theorem 1] to a space-time setting.

To this end, let us define the continuous function $\psi$, by setting

$$ \left\{ \begin{array}{ll} \psi(s) = -\epsilon, & \text{for } s < -\epsilon, \\ \psi(s) = s, & \text{for } -\epsilon \le s \le \epsilon, \\ \psi(s) = \epsilon, & \text{for } s > \epsilon, \end{array} \right. $$

for $\epsilon > 0$. Given $\gamma, \hat{\gamma} > 1$, we compute

$$ \iint_{\Omega_T} |\nabla p_\gamma - \nabla p_{\hat{\gamma}}|^2 \psi'(p_\gamma - p_{\hat{\gamma}}) \, dx \, dt = - \iint_{\Omega_T} (\Delta p_\gamma - \Delta p_{\hat{\gamma}}) \psi(p_\gamma - p_{\hat{\gamma}}) \, dx \, dt. $$

Next we split the domain into two parts by defining the set

$$ \Omega_{T,\epsilon} := \{(x, t) \in \Omega_T \mid |p_\gamma(x, t) - p_{\hat{\gamma}}(x, t)| \le \epsilon\}. $$

Thus, since $\Delta p_\gamma$ is bounded in $L^1(Q_T)$ (uniformly with respect to $\gamma$), we have

$$ \iint_{\Omega_{T,\epsilon}} |\nabla p_\gamma - \nabla p_{\hat{\gamma}}|^2 \, dx \, dt \le C\epsilon. $$
Hence

$$ \iint_{\Omega_T} |\nabla p_{\gamma} - \nabla p_{\hat{\gamma}}| \, dx \, dt = \iint_{\Omega_{T,\epsilon}} |\nabla p_{\gamma} - \nabla p_{\hat{\gamma}}| \, dx \, dt + \iint_{\Omega_{T,\epsilon}^c} |\nabla p_{\gamma} - \nabla p_{\hat{\gamma}}| \, dx \, dt \\ \leq C\epsilon^{1/2} + 2 T^{1/2} \|\nabla p_{\gamma}\|_{L^2(Q_T)} \cdot |\Omega_{T,\epsilon}^c|^{1/2}, $$

where in the last line we used Hölder's inequality. Since $p_γ$ is compact, it is a Cauchy sequence, and there exist $\Gamma(\epsilon)$ large enough such that for $\gamma, \hat{\gamma} > \Gamma(\epsilon)$ there holds

$$ \iint_{\Omega_T} |\nabla p_{\gamma} - \nabla p_{\hat{\gamma}}| \, dx \, dt \leq C\epsilon^{1/2} + C\epsilon. $$

This implies that $\nabla p_γ$ is a Cauchy sequence in $L^1(Q_T)$. Up to a subsequence we have a.e. convergence. Thanks to Eq. (20), the pressure gradient is compact in any $L^r(Q_T)$, for $1 \leq r < 4$. $\square$

**Remark 3.7.** The tumour growth rate usually depends also on the presence of nutrients, therefore one can couple Eq. (1), with an equation on the nutrient concentration. Then, the model reads

$$ (21) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial n_{\gamma}}{\partial t} - \nabla \cdot (n_{\gamma} \nabla p_{\gamma}) - \nabla \cdot (n_{\gamma} \nabla \Phi) = n_{\gamma} G(p_{\gamma}, c_{\gamma}), \\[2ex] \displaystyle \frac{\partial c_{\gamma}}{\partial t} - \Delta c_{\gamma} = -n_{\gamma} H(c_{\gamma}), \end{array} \right. $$

where $H$ is the nutrient consumption rate. Thus, system (21) is actually an extension of the model with nutrient studied in [50].

Let us notice that the proofs of the estimates in Lemma 3.2 and Lemma 3.3 can be adapted for system (21) without any particular difficulty. In fact, the boundedness of the new terms depending on $c_γ$, $\nabla c_γ$, and $\Delta c_γ$ relies only on the $L^2$-regularity of $c_γ$ and its derivatives, which comes directly from its equation in system (21). Therefore, the strong convergence stated in Lemma 3.6 still holds for this model. We refer the reader to [50] and [22] for the complete treatment of these additional terms.

# 4 THE INCOMPRESSIBLE LIMIT

The results obtained in Section 3 allow us to finally pass to the incompressible limit in Eq. (2) and obtain the complementarity relation, Eq. (6). Let us point out that, thanks to the uniform (with respect to $\gamma$) boundness of $\nabla p_γ$ in $L^2(Q_T)$ and $\partial_t p_γ$ in $L^1(Q_T)$, the complementarity relation turns out to be equivalent to the strong convergence of $\nabla p_γ$ in $L^2(Q_T)$, given by Lemma 3.6.

**Theorem 4.1 (Complementarity relation).** We may pass to the limit in Eq. (2), as $\gamma \to \infty$, and obtain the so-called complementarity relation

$$ p_{\infty}(\Delta p_{\infty} + \Delta \Phi + G(p_{\infty})) = 0, $$

in the distributional sense. Moreover, $n_∞$ and $p_∞$ satisfy the equations

$$ (22a) \qquad \frac{\partial n_{\infty}}{\partial t} = \Delta p_{\infty} + n_{\infty} G(p_{\infty}) + \nabla \cdot (n_{\infty} \nabla \Phi), $$

in $D'(Q_T)$, as well as

$$ (22b) \qquad p_{\infty}(1 - n_{\infty}) = 0, $$

almost everywhere.

*Proof.* Thanks to the bounds in Lemma 2.1,

$$ \iint_{\Omega_T} \left| \frac{\partial p_{\gamma}}{\partial t} \right| + |\nabla p_{\gamma}| \, dx \, dt \leq C(T), $$

then, by the Fréchet-Kolmogorov Theorem, $p_γ$ is strongly compact in $L^1(Q_T)$, for all $T > 0$.
We integrate Eq. (2) against a test function $\varphi \in C_c^\infty(Q_T)$ to obtain

$$ \iint_{Q_T} \frac{\partial p_\gamma}{\partial t} \varphi \, dx \, dt = (1-\gamma) \left( \iint_{Q_T} |\nabla p_\gamma|^2 \varphi \, dx \, dt + \iint_{Q_T} \nabla p_\gamma \cdot \nabla \Phi \varphi \, dx \, dt \right) \\ - \gamma \iint_{Q_T} p_\gamma \nabla p_\gamma \cdot \nabla \varphi \, dx \, dt - \gamma \iint_{Q_T} p_\gamma \nabla \Phi \cdot \nabla \varphi \, dx \, dt \\ + \gamma \iint_{Q_T} p_\gamma G(p_\gamma) \varphi \, dx \, dt. $$

Dividing by $\gamma - 1$ and passing to the limit $\gamma \to \infty$, we obtain

$$ \lim_{\gamma \to \infty} \left[ - \iint_{Q_T} (|\nabla p_\gamma|^2 \varphi + p_\gamma \nabla p_\gamma \cdot \nabla \varphi) \, dx \, dt \right. \\ \left. - \iint_{Q_T} (\nabla p_\gamma \cdot \nabla \Phi \varphi + p_\gamma \nabla \Phi \cdot \nabla \varphi) \, dx \, dt + \iint_{Q_T} p_\gamma G(p_\gamma) \varphi \, dx \, dt \right] = 0. $$

It remains to identify the limit. By the strong convergence of $p_\gamma$ and $\nabla p_\gamma$ in $L^2(Q_T)$ we have

$$ - \iint_{Q_T} (|\nabla p_\infty|^2 \varphi + p_\infty \nabla p_\infty \cdot \nabla \varphi) \, dx \, dt - \iint_{Q_T} (\nabla p_\infty \cdot \nabla \Phi \varphi + p_\infty \nabla \Phi \cdot \nabla \varphi) \, dx \, dt \\ + \iint_{Q_T} p_\infty G(p_\infty) \varphi \, dx \, dt = 0, $$

i.e.,

$$ p_\infty(\Delta p_\infty + \Delta \Phi + G(p_\infty)) = 0, $$

in the distributional sense.

Now, we prove that Eq. (22a) and Eq. (22b) are satisfied. By Lemma 2.1, we have

$$ \iint_{\Omega_T} \left| \frac{\partial n_\gamma}{\partial t} \right| + |\nabla n_\gamma| \, dx \, dt \le C(T), $$

and then we infer the compactness of the density. Up to a subsequence, we also have almost everywhere convergence, both for $n_\gamma$ and $p_\gamma$. Passing to the limit in the relation $p_\gamma^{(1+\gamma)/\gamma} = n_\gamma p_\gamma$, we obtain

$$ p_\infty(1 - n_\infty) = 0, $$

a.e. in $Q_T$.

Now, we may pass to the limit in the distributional sense in Eq. (1) to obtain

$$ \frac{\partial n_\infty}{\partial t} = \nabla \cdot (n_\infty \nabla p_\infty) + n_\infty G(p_\infty) + \nabla \cdot (n_\infty \nabla \Phi). $$

From the following relation

$$ \frac{1+\gamma}{\gamma} n_\gamma \nabla p_\gamma = p_\gamma \nabla n_\gamma + n_\gamma \nabla p_\gamma, $$

we infer

$$ \frac{1}{\gamma} n_\gamma \nabla p_\gamma = p_\gamma \nabla n_\gamma, $$
and therefore $p_γ∇n_γ → 0$ strongly in $L^1(Q_T)$ as $γ → ∞$. Consequently, for any $φ ∈ C_c^∞(Q_T)$ we have

$$ 
\begin{align*}
\iint_{Q_T} \nabla \cdot (n_\gamma \nabla p_\gamma) \varphi \, dx \, dt &= - \iint_{Q_T} n_\gamma \nabla p_\gamma \cdot \nabla \varphi \, dx \, dt \\
&= \iint_{Q_T} n_\gamma p_\gamma \Delta \varphi \, dx \, dt + \iint_{Q_T} p_\gamma \nabla n_\gamma \cdot \nabla \varphi \, dx \, dt \\
&\rightarrow \iint_{Q_T} n_\infty p_\infty \Delta \varphi \, dx \, dt \\
&= \iint_{Q_T} p_\infty \Delta \varphi \, dx \, dt.
\end{align*}
$$

As a consequence, $n_∞$ and $p_∞$ satisfy

$$ \frac{\partial n_{\infty}}{\partial t} = \Delta p_{\infty} + n_{\infty} G(p_{\infty}) + \nabla \cdot (n_{\infty} \nabla \Phi), $$

which completes the proof. □

# 5 UNIQUENESS OF THE LIMIT PRESSURE

This section is dedicated to proving the following statement.

**Theorem 5.1 (Uniqueness of $n_∞$ and $p_∞$).** *The incompressible limit obtained in the previous section, ($n_∞, p_∞$), cf. Eq. (7a) is unique.*

*Proof.* In order to prove uniqueness, we assume that $(n_1, p_1)$ and $(n_2, p_2)$ are two solutions and let $Ω$ be a compact, simply connected Lipschitz set that contains the union of their supports. Upon subtracting the equation for $n_2$ from the equation for $n_1$ we see that difference, $n_1 - n_2$, satisfies

$$ (23) \qquad \frac{\partial(n_1 - n_2)}{\partial t} - \Delta(p_1 - p_2) - \nabla \cdot ((n_1 - n_2)\nabla\Phi) - (n_1 G(p_1) - n_2 G(p_2)) = 0. $$

For the sake of simplicity, we shall use the short-hand notation $G_i = G(p_i)$, for $i = 1, 2$, and $v = \nabla\Phi$. Multiplying Eq. (23) by a test function $\psi = \psi(x, t)$ and integrating by parts we get

$$ (24) \qquad \iint_{\Omega_T} \left[ (n_1 - n_2) \frac{\partial \psi}{\partial t} + (p_1 - p_2) \Delta \psi - (n_1 - n_2) \nabla \psi \cdot v + (n_1 G_1 - n_2 G_2) \psi \right] dx \, dt = 0. $$

The strategy is to employ Hilbert's dual method to establish uniqueness. To this end we introduce the following notation

$$ 
\left\{ 
\begin{aligned}
\mathcal{Z} &:= n_1 - n_2 + p_1 - p_2, \\
\mathcal{A} &:= \frac{n_1 - n_2}{\mathcal{Z}}, \\
\mathcal{B} &:= \frac{p_1 - p_2}{\mathcal{Z}}, \\
\mathcal{C} &:= -n_2 \frac{G_1 - G_2}{p_1 - p_2},
\end{aligned}
\right. 
$$

where we set $\mathcal{A} = \mathcal{B} = 0$, whenever $\mathcal{Z} = 0$. Using this notation we rewrite Eq. (24) which becomes

$$ (25) \qquad \iint_{\Omega_T} \mathcal{Z} \left[ A \frac{\partial \psi}{\partial t} + B \Delta \psi - A \nabla \psi \cdot \nabla v + (A G_1 - B C) \psi \right] dx \, dt = 0. $$

Note that, by definition,

$$ 0 \le \mathcal{A}, \mathcal{B} \le 1, \quad \text{as well as} \quad 0 \le \mathcal{C} \le \sup_{0 \le p \le p_M} |G'(p)|. $$
In order to apply Hilbert's duality method, we have to find a solution, $\psi$, to the dual problem

$$ (26) \qquad \mathcal{A} \frac{\partial \psi}{\partial t} + \mathcal{B} \Delta \psi - \mathcal{A} \nabla \psi \cdot v + (\mathcal{A}G_1 - \mathcal{B}\mathcal{C})\psi = \mathcal{A}\xi, $$

in $\Omega_T$, and $\psi = 0$ on $\partial\Omega \times (0, T)$. The equation is complemented by the final time condition $\psi(x, T) = 0$ for $x \in \Omega$. Here, $\xi$ is an arbitrary smooth function. If solved, substituting the solution to the dual problem, $\psi$, into Eq. (25) would yield

$$ (27) \qquad \iint_{\Omega_T} \mathcal{A}\mathcal{Z}\xi \,dx\,dt = \iint_{\Omega_T} (n_1 - n_2)\xi \,dx\,dt = 0, $$

thus proving uniqueness of the density. Subsequently, from Eq. (24), the uniqueness of the pressure follows.

However, since the coefficient of Eq. (26) are not smooth and $\mathcal{A}$ and $\mathcal{B}$ can vanish, the equation is not uniformly parabolic and we need to regularise the system first. To this end, let $\{\mathcal{A}_k\}, \{\mathcal{B}_k\}, \{\mathcal{C}_k\}, \{v_k\}, \{\mathcal{G}_{1,k}\}$ be approximating sequences of smooth and bounded functions such that

$$ (28a) \quad \|\mathcal{A} - \mathcal{A}_k\|_{L^2(\Omega_T)}, \|\mathcal{B} - \mathcal{B}_k\|_{L^2(\Omega_T)}, \|\mathcal{C} - \mathcal{C}_k\|_{L^2(\Omega_T)}, \|G_1 - \mathcal{G}_{1,k}\|_{L^2(\Omega_T)}, \|v - v_k\|_{L^2(\Omega_T)} \le \frac{1}{k}, $$

such that

$$ (28b) \qquad 1/k \le \mathcal{A}_k, \mathcal{B}_k \le 1, \quad \text{as well as} \quad 0 \le \mathcal{C}_k, |\mathcal{G}_{1,k}| \le C, $$

and

$$ (28c) \qquad \|\partial_t \mathcal{C}_k\|_{L^1(\Omega_T)}, \|\nabla \mathcal{G}_{1,k}\|_{L^2(\Omega_T)} \le C, $$

where $C > 0$ is some positive constant. Using the regularised quantities, we consider the regularised equation

$$ (29) \qquad \frac{\partial \psi_k}{\partial t} + \frac{\mathcal{B}_k}{\mathcal{A}_k} \Delta \psi_k - \nabla \psi_k \cdot v_k + \left( \mathcal{G}_{1,k} - \frac{\mathcal{B}_k \mathcal{C}_k}{\mathcal{A}_k} \right) \psi_k = \xi, $$

in $\Omega_T$, and $\psi_k = 0$, on $\partial\Omega \times (0, T)$, and $\psi_k(T, x) = 0$, in $\Omega$. Here, $\xi$ denotes an arbitrary smooth test function which is crucial for this approach, as discussed above, cf. Eq. (27). Since the coefficient $\mathcal{B}_k/\mathcal{A}_k$ is smooth and bounded from away from zero, the equation is uniformly parabolic, whence we infer the existence of a smooth solution, $\psi_k$.

Using $\psi_k$ as a test function in Eq. (25) and thanks to Eq. (29) we get

$$ \begin{aligned} 0 &= \iint_{\Omega_T} \mathcal{Z} \left( \mathcal{A} \frac{\partial \psi_k}{\partial t} + \mathcal{B} \Delta \psi_k - \mathcal{A} v \cdot \nabla \psi_k + (\mathcal{A}G_1 - \mathcal{B}\mathcal{C})\psi_k \right) \,dx\,dt \\ &= \iint_{\Omega_T} \mathcal{Z} \mathcal{A} \left( -\frac{\mathcal{B}_k}{\mathcal{A}_k} \Delta \psi_k + v_k \cdot \nabla \psi_k - \left( \mathcal{G}_{1,k} - \frac{\mathcal{B}_k \mathcal{C}_k}{\mathcal{A}_k} \right) \psi_k + \xi \right) \,dx\,dt \\ &\quad + \iint_{\Omega_T} \mathcal{Z} (\mathcal{B} \Delta \psi_k - \mathcal{A} v \cdot \nabla \psi_k + (\mathcal{A}G_1 - \mathcal{B}\mathcal{C})\psi_k) \,dx\,dt \\ &= \iint_{\Omega_T} \mathcal{Z} \mathcal{A} \xi + \iint_{\Omega_T} \mathcal{Z} \frac{\mathcal{B}_k}{\mathcal{A}_k} (\mathcal{A} - \mathcal{A}_k) (-\Delta \psi_k + \mathcal{C}_k \psi_k) \,dx\,dt \\ &\quad + \iint_{\Omega_T} \mathcal{Z} (\mathcal{B}_k - \mathcal{B}) (-\Delta \psi_k + \mathcal{C}_k \psi_k) \,dx\,dt + \iint_{\Omega_T} \mathcal{Z} \mathcal{B} (\Delta \psi_k - \mathcal{C} \psi_k) \,dx\,dt \\ &\quad + \iint_{\Omega_T} \mathcal{Z} \mathcal{B} (-\Delta \psi_k + \mathcal{C}_k \psi_k) \,dx\,dt + \iint_{\Omega_T} \mathcal{Z} \mathcal{A} \psi_k (G_1 - \mathcal{G}_{1,k}) \,dx\,dt \\ &\quad + \iint_{\Omega_T} \mathcal{Z} \mathcal{A} \nabla \psi_k \cdot (v_k - v) \,dx\,dt. \end{aligned} $$
Using the definition of $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{Z}$, we finally obtain

$$\iint_{\Omega_T} (n_1 - n_2) \xi \, dx \, dt = I_k^1 - I_k^2 + I_k^3 - I_k^4 + I_k^5,$$

where

$$I_k^1 = \iint_{\Omega_T} (n_1 - n_2 + p_1 - p_2) \frac{\mathcal{B}_k}{\mathcal{A}_k} (\mathcal{A} - \mathcal{A}_k) (\Delta\psi_k - C_k\psi_k) \, dx \, dt,$$

$$I_k^2 = \iint_{\Omega_T} (n_1 - n_2 + p_1 - p_2) (\mathcal{B} - \mathcal{B}_k) (\Delta\psi_k - C_k\psi_k) \, dx \, dt,$$

$$I_k^3 = \iint_{\Omega_T} (p_1 - p_2) (\mathcal{C} - \mathcal{C}_k) \psi_k \, dx \, dt,$$

$$I_k^4 = \iint_{\Omega_T} (n_1 - n_2) (G_1 - G_{1,k}) \psi_k \, dx \, dt,$$

$$I_k^5 = \iint_{\Omega_T} (n_1 - n_2) \nabla \psi_k \cdot (v - v_k) \, dx \, dt.$$

We aim at showing that

$$\lim_{k \to \infty} I_k^i = 0,$$

for $i = 1, \dots, 5$, in order to be able to conclude that $n_1 = n_2$. Before proving the convergence of each $I_k^i$, we need certain uniform bounds which we collect and state in the subsequent lemma.

**Lemma 5.2 (Uniform bounds).** There exist a positive constant $C > 0$, independent of $k$, such that

$$\sup_{0 \le t \le T} \|\psi_k(t)\|_{L^\infty(\Omega)} \le C, \quad \sup_{0 \le t \le T} \|\nabla \psi_k(t)\|_{L^2(\Omega)} \le C, \\ \|(\mathcal{B}_k / \mathcal{A}_k)^{1/2} (\Delta \psi_k - C_k \psi_k)\|_{L^2(\Omega_T)} \le C.$$

*Proof.* The $L^\infty$-bound comes directly from the maximum principle applied to Eq. (29), since $\xi$ is bounded and

$$\mathcal{G}_{1,k} - \frac{\mathcal{B}_k \mathcal{C}_k}{\mathcal{A}_k} \le C.$$

Now we multiply Eq. (29) by $(\Delta\psi_k - C_k\psi_k)$ and integrate in $(t, T) \times \Omega$ to obtain

$$\begin{align}
& - \int_t^T \int_\Omega \frac{\partial}{\partial t} \frac{|\nabla \psi_k|^2}{2} \, dx \, ds - \int_t^T \int_\Omega \frac{C_k}{2} \frac{\partial}{\partial t} \psi_k^2 \, dx \, ds + \int_t^T \int_\Omega \frac{\mathcal{B}_k}{\mathcal{A}_k} |\Delta \psi_k - C_k \psi_k|^2 \, dx \, ds \nonumber \\
& \qquad = \underbrace{\int_t^T \int_\Omega v \cdot \nabla \psi_k (\Delta \psi_k - C_k \psi_k) \, dx \, ds}_{\mathcal{I}_1} - \underbrace{\int_t^T \int_\Omega \mathcal{G}_{1,k} \psi_k (\Delta \psi_k - C_k \psi_k) \, dx \, ds}_{\mathcal{I}_2} \nonumber \\
& \qquad + \underbrace{\int_t^T \int_\Omega \xi (\Delta \psi_k - C_k \psi_k) \, dx \, ds}_{\mathcal{I}_3}, \tag{31}
\end{align}$$

where we shall bound each of the terms, $\mathcal{I}_i$, for $i = 1, 2, 3$, individually. First note that

$$\begin{align*}
\mathcal{I}_1 &= \int_t^T \int_\Omega v \cdot \nabla \psi_k \Delta \psi_k \, dx \, ds - \int_t^T \int_\Omega v \cdot \nabla \psi_k C_k \psi_k \, dx \, ds \\
&= \mathcal{I}_{1,1} + \mathcal{I}_{1,2}.
\end{align*}$$
Integrating by parts in the first term of $\mathcal{I}_1$ we get

$$
\begin{align*}
\mathcal{I}_{1,1} &= - \int_t^T \int_\Omega \sum_{i,j=1}^d \frac{\partial v^{(i)}}{\partial x_j} \frac{\partial \psi_k}{\partial x_i} \frac{\partial \psi_n}{\partial x_j} \mathrm{d}x \mathrm{d}s - \int_t^T \int_\Omega \sum_{i,j=1}^d v^{(i)} \frac{\partial^2 \psi_k}{\partial x_i \partial x_j} \frac{\partial \psi_k}{\partial x_j} \mathrm{d}x \mathrm{d}s \\
&= - \int_t^T \int_\Omega \sum_{i,j=1}^d \frac{\partial v^{(i)}}{\partial x_j} \frac{\partial \psi_k}{\partial x_i} \frac{\partial \psi_n}{\partial x_j} \mathrm{d}x \mathrm{d}s + \int_t^T \int_\Omega \frac{|\nabla \psi_k|^2}{2} \nabla \cdot v \mathrm{d}x \mathrm{d}s \\
&\le \left( d \|\nabla v\|_{L^\infty} + \frac{1}{2} \|\nabla \cdot v\|_{L^\infty} \right) \int_t^T \int_\Omega |\nabla \psi_k|^2 \mathrm{d}x \mathrm{d}s,
\end{align*}
$$

where $v^{(i)}$ is the $i$-th component of the vector $v$ and $\nabla v$ is the matrix with element $(\nabla v)_{i,j} = \partial_j v^{(i)}$. Similarly, we observe

$$
\begin{align*}
\mathcal{I}_{1,2} &= - \int_t^T \int_\Omega v \cdot \nabla \psi_k \mathcal{C}_k \psi_k \, \mathrm{d}x \, \mathrm{d}s \\
&\leq \frac{1}{2} \|v\|_{L^\infty(\Omega_T)} \|\mathcal{C}_k\|_{L^\infty(\Omega_T)} \|\psi_k\|_{L^2(\Omega_T)}^2 + \frac{1}{2} \|\nabla \psi_k\|_{L^2(\Omega_T)}^2 \\
&\leq C + C \|\nabla \psi_k\|_{L^2(\Omega_T)}^2,
\end{align*}
$$

with $C > 0$ independent of $k$, after applying Young's inequality. Hence

$$
\mathcal{I}_1 \leq C + C \|\nabla \psi_k\|_{L^2(\Omega_T)}^2.
$$

Next, let us address the term $\mathcal{I}_2$. We observe that

$$
\begin{align*}
\mathcal{I}_2 &= - \int_t^T \int_\Omega \mathcal{G}_{1,k} \psi_k (\Delta \psi_k - \mathcal{C}_k \psi_k) \, \mathrm{d}x \, \mathrm{d}s \\
&= \int_t^T \int_\Omega \mathcal{G}_{1,k} |\nabla \psi_k|^2 \, \mathrm{d}x \, \mathrm{d}s + \int_t^T \int_\Omega \psi_k \nabla \psi_k \cdot \nabla \mathcal{G}_{1,k} \, \mathrm{d}x \, \mathrm{d}s + \int_t^T \int_\Omega \mathcal{G}_{1,k} \mathcal{C}_k \psi_k \, \mathrm{d}x \, \mathrm{d}s.
\end{align*}
$$

We note that $\|\mathcal{G}_{1,k}\|_{L^\infty(\Omega_T)}$ whence we obtain bounds for the first and the last term, respectively.
In addition, we recall $\|\nabla \mathcal{G}_{1,k}\|_{L^2(\Omega_T)} \le C$, whence, upon using Young's inequality, we get

$$
\begin{align*}
\int_t^T \int_\Omega \psi_k \nabla \psi_k \cdot \nabla \mathcal{G}_{1,k} \, \mathrm{d}x \, \mathrm{d}s &\le \frac{1}{2} \|\psi_k\|_{L^\infty(\Omega_T)} \|\nabla \mathcal{G}_{1,k}\|_{L^2(\Omega_T)}^2 + \frac{1}{2} \|\psi_k\|_{L^\infty(\Omega_T)} \int_t^T \int_\Omega |\nabla \psi_k|^2 \, \mathrm{d}x \, \mathrm{d}s \\
&\le C + C \int_t^T \int_\Omega |\nabla \psi_k|^2 \, \mathrm{d}x \, \mathrm{d}s.
\end{align*}
$$

In combination we get

$$
\mathcal{I}_2 \leq C + C \|\nabla \psi_k\|_{L^2(\Omega_T)}^2,
$$

with $C > 0$ independent of $k$. Last, let us address the term $\mathcal{I}_3$. We readily observe

$$
\begin{align*}
\mathcal{I}_3 &= \int_t^T \int_\Omega \xi (\Delta \psi_k - \mathcal{C}_k \psi_k) \, \mathrm{d}x \, \mathrm{d}s \\
&\le C,
\end{align*}
$$
integrating by parts twice and using the $L^\infty$-bounds. Using the bounds obtained above, the right-hand side of Eq. (31) can be bounded as follows

$$
\begin{align*}
C+C\|\nabla\psi_k\|_{L^2(\Omega_T)}^2
&\geq -\int_t^T \int_\Omega \frac{\partial}{\partial t} \frac{|\nabla\psi_k|^2}{2} \,dx\,ds - \int_t^T \int_\Omega \frac{C_k}{2} \frac{\partial}{\partial t} \psi_k^2 \,dx\,ds + \int_t^T \int_\Omega \frac{B_k}{A_k} |\Delta\psi_k - C_n\psi_k|^2 \,dx\,ds \\
&\geq -\int_t^T \frac{d}{dt} \int_\Omega \frac{|\nabla\psi_k|^2}{2} \,dx\,ds + \int_t^T \int_\Omega \frac{\partial C_k}{\partial t} \frac{\psi_k^2}{2} \,dx\,ds + \int_t^T \int_\Omega \frac{B_k}{A_k} |\Delta\psi_k - C_k\psi_k|^2 \,dx\,ds \\
&\phantom{\geq} + \int_\Omega \frac{C_k(t)\psi_k^2(t)}{2} \,dx \\
&\geq \frac{1}{2} \|\nabla\psi_k(\cdot, t)\|_{L^2(\Omega)}^2 - \|\partial_t C_k\|_{L^1(\Omega_T)} \|\psi_k\|_{L^\infty(\Omega_T)}^2 + \int_t^T \int_\Omega \frac{B_k}{A_k} |\Delta\psi_k - C_k\psi_k|^2 \,dx\,ds \\
&\phantom{\geq} - \frac{1}{2} \|C_k\|_{L^\infty(\Omega_T)} \|\psi_k\|_{L^2(\Omega_T)}^2 \\
&\geq \frac{1}{2} \|\nabla\psi_k(\cdot, t)\|_{L^2(\Omega)}^2 + \int_t^T \int_\Omega \frac{B_k}{A_k} |\Delta\psi_k - C_k\psi_k|^2 \,dx\,ds - C,
\end{align*}
$$

having used the regularity assumptions on the regularised coefficients, cf. Eq. (28). Finally, since $C_k$ is positive, we get

$$ (32) \quad \frac{1}{2} \int_{\Omega} |\nabla \psi_k(t)|^2 dx + \int_{t}^{T} \int_{\Omega} \frac{B_k}{A_k} |\Delta \psi_k - C_k \psi_k|^2 dx ds \leq C + C \int_{t}^{T} \int_{\Omega} |\nabla \psi_k|^2 dx ds. $$

Introducing the notation

$$ Q(s) := \int_{\Omega} |\nabla \psi_k(s, x)|^2 dx, $$

we observe that Eq. (32) now reads

$$ Q(t) \leq C + C \int_{t}^{T} Q(s) \, ds, $$

and by Gronwall's lemma we conclude that

$$ \sup_{0 \le t \le T} Q(t) = \sup_{0 \le t \le T} \|\nabla \psi_k(t)\|_{L^2(\Omega)}^2 \le C. $$

The third bound of Eq. (30) comes a posteriori from Eq. (32), which completes proof. $\square$

Thanks to these uniform bounds, we obtain

$$
\begin{align*}
I_k^1 &= \iint_{\Omega_T} (n_1 - n_2 + p_1 - p_2) \frac{B_k}{A_k} (A - A_k) (\Delta \psi_k - C_k \psi_k) \,dx\,dt \\
&\leq C \|(B_k/A_k)^{1/2} (A - A_k)\|_{L^2(\Omega_T)} \\
&\leq C \sqrt{k} \|A - A_k\|_{L^2(\Omega_T)} \\
&\leq \frac{C}{\sqrt{k}},
\end{align*}
$$

and, similarly,

$$ I_n^2 = \iint_{\Omega_T} (n_1 - n_2 + p_1 - p_2)(B - B_k)(\Delta\psi_k - C_k\psi_k) \,dx\,dt \leq \frac{C}{\sqrt{k}} \|B - B_k\|_{L^2(\Omega_T)} \leq \frac{C}{\sqrt{k}}. $$
Finally, we have

$$I_k^3 = \iint_{\Omega_T} (p_1 - p_2)(C - C_n)\psi_k \, dx \, dt \le C\|C - C_k\|_{L^2(\Omega_T)} \le \frac{C}{k},$$

and

$$I_k^4 = \iint_{\Omega_T} (n_1 - n_2)(G_1 - G_{1,k})\psi_n \, dx \, dt \le C\|G_1 - G_{1,k}\|_{L^2(\Omega_T)} \le \frac{C}{k},$$

as well as

$$I_n^5 = \iint_{\Omega_T} (n_1 - n_2)\nabla\psi_n \cdot (v - v_k) \, dx \, dt \le C\|v - v_k\|_{L^2(\Omega_T)} \le \frac{C}{k}.$$

In summary, we have

$$\iint_{\Omega_T} (n_1 - n_2)\xi \, dx \, dt = I_k^1 - I_k^2 + I_k^3 - I_k^4 + I_k^5 \to 0,$$

as $k \to \infty$, and therefore $n_1 = n_2$. From Eq. (24) we have

$$\iint_{\Omega_T} ((p_1 - p_2)\Delta\psi + n_1(G(p_1) - G(p_2))\psi) \, dx \, dt = 0.$$

Taking a smooth approximation of $p_1 - p_2$ as test function we get

$$\iint_{\Omega_T} |\nabla(p_1 - p_2)|^2 dx \, dt = \iint_{\Omega_T} n_1(G(p_1) - G(p_2))(p_1 - p_2) dx \, dt,$$

and, by the monotonicity of $G$, cf. Eq. (A-G), we conclude that $\nabla(p_1 - p_2)$, almost everywhere. Substituting this, in conjunction with $n_1 = n_2$, into Eq. (24), we get

$$\iint_{\Omega_T} n_1(G_1 - G_2)\psi \, dx \, dt = 0,$$

whence $G(p_1) = G(p_2)$, almost everywhere on $n_1 > 0$. Since $G$ is strictly decreasing we have $p_1 = p_2$. In the case $n_1 = 0$, the uniqueness follows from the relation $p_\infty(n_\infty - 1) = 0$. $\square$

# 6 VELOCITY OF THE BOUNDARY FOR PATCHES

Let us recall that the Hele-Shaw problem is given by

$$(33) \qquad \begin{cases} -\Delta p_{\infty} = \Delta \Phi + G(p_{\infty}), & \text{in } \Omega(t), \\ V = -(\nabla p_{\infty} + \nabla \Phi) \cdot \nu, & \text{on } \partial\Omega(t), \end{cases}$$

where $\nu$ indicates the outward normal to the boundary and $\Omega(t) := \{x; p_\infty(x, t) > 0\}$. Below we give a characterisation of patch solutions, i.e., the indicator of the growing domain described by Eq. (33) satisfies the incompressible limit equation, cf. Eq. (7a). To this end, we suppose that the boundary $\partial\Omega(t)$ admits a Lipschitz parameterisation $\partial\Omega(t) = \{x(t, \alpha) | \alpha \in [0, 1], x(t, 0) = x(t, 1)\}$ that satisfies

$$(34) \qquad \frac{d}{dt}x(t, \alpha) = -(\nabla p_{\infty}(x(t, \alpha), t) + \nabla \Phi(x(t, \alpha), t)).$$

Then the characteristic function

$$(35) \qquad n_{\infty}(t) = \mathbb{1}_{\Omega(t)}.$$

satisfies the limit problem, Eq. (7a).

**Theorem 6.1** (Characterisation of the Free Boundary Velocity). Let $\Omega_0$ be a bounded and Lipschitz continuous domain. Let us consider the solution $(\Omega(t), p_\infty)$ to the free boundary problem, Eq. (33), with initial data $\Omega_0$. Then, the characteristic function in Eq. (35), satisfies Eq. (7a).
Proof. We have to show that $n_{\infty}(t) = 1_{\Omega(t)}$ satisfies

$$ \frac{\partial n_{\infty}}{\partial t} = \Delta p_{\infty} + \nabla \cdot (n_{\infty} \nabla \Phi) + n_{\infty} G(p_{\infty}), $$

in the distributional sense. Given a test function $\psi = \psi(x)$, by Reynolds' transport Theorem and Eq. (34), we have

$$ \int_{\mathbb{R}^d} \psi(x) \frac{\partial n_{\infty}}{\partial t} dx = \frac{d}{dt} \int_{\mathbb{R}^d} \psi(x) 1_{\Omega(t)} dx = \int_{\partial \Omega(t)} V \psi(x) dx = V \delta_{\partial \Omega(t)}. $$

On the other hand, it holds

$$ \Delta p_{\infty} + \nabla \cdot (n_{\infty} \nabla \Phi) + n_{\infty} G(p_{\infty}) = -(\partial_{\nu} p_{\infty} + \partial_{\nu} \Phi) \delta_{\partial \Omega(t)} = V \delta_{\partial \Omega(t)}, $$

in the sense of distributions, as can be seen by the following argument. First, by the definition of $\Omega(t)$ as the positivity set of $p_{\infty}$ and the fact that $n_{\infty} = 1_{\Omega(t)}$ we observe that the weak formulation of the left-hand side can be manipulated as follows:

$$ \int_{\mathbb{R}^d} -\nabla p_{\infty} \cdot \nabla \psi - n_{\infty} \nabla \Phi \cdot \nabla \psi + n_{\infty} G(p_{\infty}) \psi \, dx = \int_{\Omega(t)} -\nabla p_{\infty} \cdot \nabla \psi - \nabla \Phi \cdot \nabla \psi + G(p_{\infty}) \psi \, dx. $$

Integrating by parts the right-hand side, we obtain

$$ \int_{\Omega(t)} (\Delta p_{\infty} + \Delta \Phi + G(p_{\infty})) \psi \, dx - \int_{\partial \Omega(t)} \partial_{\nu} p_{\infty} \psi \, dx - \int_{\partial \Omega(t)} \partial_{\nu} \Phi \psi \, dx \\ = - \int_{\partial \Omega(t)} \partial_{\nu} p_{\infty} \psi \, dx - \int_{\partial \Omega(t)} \partial_{\nu} \Phi \psi \, dx $$

where we used $\Delta p_{\infty} + \Delta \Phi + G(p_{\infty}) = 0$, in $D'$, by Eq. (33). $\square$

## ACKNOWLEDGEMENTS

This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie (grant agreement No 754362). The authors are grateful for the discussions with Benoît Perthame at several occasions. Moreover, the authors would like to express their gratitude to Filippo Santambrogio for his interest in the work and his suggestions.

## REFERENCES

[1] D. Alexander, I. Kim, and Y. Yao. Quasi-static evolution and congested crowd transport. *Nonlinearity*, 27(4):823–858, 2014.

[2] D. G. Aronson and P. Bénilan. Régularité des solutions de l'équation des milieux poreux dans $\mathbb{R}^N$. C. R. Acad. Sci. Paris Sér. A-B, 288(2):A103–A105, 1979.

[3] F. B. Belgacem and P.-E. Jabin. Compactness for nonlinear continuity equations. *J. Funct. Anal.*, 264(1):139–168, 2013.

[4] P. Bénilan and M. G. Crandall. The continuous dependence on $\varphi$ of solutions of $u_t - \Delta\varphi(u) = 0$. *Indiana University Mathematics Journal*, 30(2):161–177, 1981.

[5] F. Berthelin, P. Degond, M. Delitala, and M. Rascle. A model for the formation and evolution of traffic jams. *Archive for Rational Mechanics and Analysis*, 187(2):185–220, 2008.

[6] F. Berthelin, P. Degond, V. Le Blanc, S. Moutari, M. Rascle, and J. Royer. A traffic-flow model with constraints for the modeling of traffic jams. *Mathematical Models and Methods in Applied Sciences*, 18(supp01):1269–1298, 2008.

[7] G. Bevilacqua, B. Perthame, and M. Schmidchen. The Aronson-Bénilan Estimate in Lebesgue Spaces. *ArXiv preprint arXiv: 2007.15267*, 2020.

[8] L. Boccardo and T. Gallouët. Non linear elliptic and parabolic equations involving measure data. *Jour. Funct. Anal.*, 87:149–169, 1989.

[9] D. Bresch and P.-E. Jabin. Global weak solutions of PDEs for compressible media: a compactness criterion to cover new physical situations. In *Shocks, singularities and oscillations in nonlinear optics and fluid mechanics*, volume 17 of Springer INdAM Ser., pages 33–54. Springer, Cham, 2017.
[10] D. Bresch and P.-E. Jabin. Global existence of weak solutions for compressible Navier-Stokes equations: thermodynamically unstable pressure and anisotropic viscous stress tensor. *Ann. of Math.* (2), 188(2):577-684, 2018.

[11] F. Bubba, B. Perthame, C. Pouchol, and M. Schmidtchen. Hele-Shaw limit for a system of two reaction-(cross-)diffusion equations for living tissues. *Arch. Rational. Mech. Anal.*, 236:735-766, 2020.

[12] H. M. Byrne and D. Drasdo. Individual-based and continuum models of growing cell populations: a comparison. *J. Math. Biol.*, 58(4-5):657-687, 2009.

[13] L. A. Caffarelli and A. Friedman. Asymptotic behavior of solutions of $u_t = \Delta u^m$ as $m \to \infty$. *Indiana University mathematics journal*, 36(4):711-728, 1987.

[14] J. A. Carrillo, K. Craig, and Y. Yao. *Aggregation-Diffusion Equations: Dynamics, Asymptotics, and Singular Limits*, pages 65–108. Springer International Publishing, Cham, 2019.

[15] J. A. Carrillo and R. S. Gvalani. Phase transitions for nonlinear nonlocal aggregation-diffusion equations. *Communications in Mathematical Physics*, 382(1):485-545, 2021.

[16] L. Chizat and S. Di Marino. A tumor growth model of Hele-Shaw type as a gradient flow. *ESAIM Control Optim. Calc. Var.*, 26:103, 2020.

[17] K. Craig, I. Kim, and Y. Yao. Congested aggregation via Newtonian interaction. *Arch. Ration. Mech. Anal.*, 227(1):1-67, 2018.

[18] K. Craig and I. Topaloglu. Aggregation-diffusion to constrained interaction: Minimizers & gradient flows in the slow diffusion limit. *Annales de l'Institut Henri Poincaré C, Analyse non linéaire*, 37(2):239-279, 2020.

[19] M. G. Crandall and P. M. Regularizing effects for $u_t = \Delta\varphi(u)$. *Trans. Amer. Math. Soc.*, 274(1):159-168, 1982.

[20] A.-L. Dalibard, G. Lopez-Ruiz, and C. Perrin. Traveling waves for the porous medium equation in the incompressible limit: asymptotic behavior and nonlinear stability. *arXiv preprint arXiv:2108.10563*, 2021.

[21] N. David, T. Dębiec, and B. Perthame. Convergence rate for the incompressible limit of nonlinear diffusion-advective equations. *arXiv preprint arXiv:2108.00787*, 2021.

[22] N. David and B. Perthame. Free boundary limit of a tumor growth model with nutrient. *Journal de Mathématiques Pures et Appliquées*, 155:62-82, 2021.

[23] G. De Philippis, A. R. Mészáros, F. Santambrogio, and B. Velichkov. Bv estimates in optimal transportation and applications. *Archive for Rational Mechanics and Analysis*, 219(2):829-860, 2016.

[24] T. Dębiec, B. Perthame, M. Schmidtchen, and N. Vauchelet. Incompressible limit for a two-species model with coupling through brinkman's law in any dimension. *Journal de Mathématiques Pures et Appliquées*, 2020.

[25] T. Dębiec and M. Schmidtchen. Incompressible limit for a two-species tumour model with coupling through brinkman's law in one dimension. *Acta Applicandae Mathematicae*, 2020.

[26] P. Degond, S. Hecht, and N. Vauchelet. Incompressible limit of a continuum model of tissue growth for two cell populations. *Networks & Heterogeneous Media*, 15(1):57-85, 2020.

[27] P. Degond, J. Hua, and L. Navoret. Numerical simulations of the Euler system with convection constraint. *Journal of Computational Physics*, 230(22):8057-8088, 2011.

[28] S. Di Marino and A. R. Mészáros. Uniqueness issues for evolution equations with density constraints. *Mathematical Models and Methods in Applied Sciences*, 26(09):1761-1783, 2016.

[29] X. Dou, J.-G. Liu, and Z. Zhou. Modeling the autophagic effect in tumor growth: a cross diffusion model and its free boundary limit. *arXiv preprint arXiv:2007.13543*, 2020.

[30] C. Elliott, M. Herrero, J. King, and J. Ockendon. The mesa problem: Diffusion patterns for $u_t = \nabla \cdot (u^m \nabla u)$ as $m \to \infty$. *IMA journal of applied mathematics*, 37(2):147-154, 1986.

[31] O. Gil and F. Quirós. Convergence of the porous media equation to Hele-Shaw. *Nonlinear Analysis: Theory, Methods & Applications*, 44(8):1111-1131, 2001.

[32] O. Gil and F. Quirós. Boundary layer formation in the transition from the porous media equation to a hele-shaw flow. In *Annales de l'IHP Analyse non linéaire*, volume 20, pages 13-36, 2003.
[33] D. Gilbarg and N. S. Trudinger. *Elliptic partial differential equations of second order*, volume 224 of *Grundlehren der Mathematischen Wissenschaften* [Fundamental Principles of Mathematical Sciences]. Springer-Verlag, Berlin, second edition, 1983.

[34] N. Guillen, I. Kim, and A. Mellet. A hele-shaw limit without monotonicity. *ArXiv preprint arXiv: 2012.02365*, 2020.

[35] P. Gwiazda, B. Perthame, and A. Świerczewska Gwiazda. A two-species hyperbolic–parabolic model of tissue growth. *Communications in Partial Differential Equations*, 44(12):1605–1618, 2019.

[36] S. Hecht and N. Vauchelet. Incompressible limit of a mechanical model for tissue growth with non-overlapping constraint. *Communications in mathematical sciences*, 15(7):1913, 2017.

[37] N. Igbida. $L^1$-theory for reaction-diffusion hele-shaw flow with linear drift. *arXiv preprint arXiv:2105.00182*, 2021.

[38] I. Kim and H. K. Lei. Degenerate diffusion with a drift potential: A viscosity solutions approach. *Discrete & Continuous Dynamical Systems - A*, 27(2):767–786, 2010.

[39] I. Kim and N. Požár. Porous medium equation to Hele-Shaw flow with general initial density. *Trans. Amer. Math. Soc.*, 370(2):873–909, 2018.

[40] I. Kim, N. Požár, and B. Woodhouse. Singular limit of the porous medium equation with a drift. *Adv. Math.*, 349:682–732, 2019.

[41] I. Kim and J. Tong. Interface dynamics in a two-phase tumor growth model. *Interfaces Free Bound.*, 23:191–304, 2021.

[42] I. Kim and Y. P. Zhang. Porous medium equation with a drift: Free boundary regularity. *Arch. Rational Mech. Anal.*, 242:1177–1228, 2021.

[43] P.-L. Lions and N. Masmoudi. On a free boundary barotropic model. In *Annales de l'Institut Henri Poincare (C) Non Linear Analysis*, volume 16, pages 373–410. Elsevier, 1999.

[44] J.-G. Liu, M. Tang, L. Wang, and Z. Zhou. Analysis and computation of some tumor growth models with nutrient: from cell density models to free boundary dynamics. *arXiv preprint arXiv:1802.00655*, 2018.

[45] J.-G. Liu, M. Tang, L. Wang, and Z. Zhou. Towards understanding the boundary propagation speeds in tumor growth models. *arXiv preprint arXiv:1910.11502*, 2019.

[46] J.-G. Liu and X. Xu. Existence and incompressible limit of a tissue growth model with autophagy. *SIAM J. Math. Anal.*, 53(5):5215–5242, 2021.

[47] A. Mellet, B. Perthame, and F. Quirós. A Hele-Shaw problem for tumor growth. *J. Funct. Anal.*, 273(10):3061–3093, 2017.

[48] O. A. Oleinik, A. S. Kalashnikov, and C. Jui-lin. The Cauchy problem and boundary problems for equations of the type of non-stationary filtration. 22(5):667–704, 1958.

[49] C. Perrin and E. Zatorska. Free/congested two-phase model from weak solutions to multi-dimensional compressible navier-stokes equations. *Communications in Partial Differential Equations*, 40(8):1558–1589, 2015.

[50] B. Perthame, F. Quirós, and J. L. Vázquez. The Hele-Shaw asymptotics for mechanical models of tumor growth. *Arch. Ration. Mech. Anal.*, 212(1):93–127, 2014.

[51] B. Perthame, F. Quirós, M. Tang, and N. Vauchelet. Derivation of a Hele-Shaw type system from a cell model with active motion. *Interfaces Free Bound.*, 16(4):489–508, 2014.

[52] B. Perthame, M. Tang, and N. Vauchelet. Traveling wave solution of the Hele-Shaw model of tumor growth with nutrient. *Math. Models Methods Appl. Sci.*, 24(13):2601–2626, 2014.

[53] B. Perthame and N. Vauchelet. Incompressible limit of a mechanical model of tumour growth with viscosity. *Philos. Trans. Roy. Soc. A*, 373(2050):20140283, 16, 2015.

[54] E. Sabinina. On the Cauchy problem for the equation of nonstationary gas filtration in several space variables. In *Doklady Akademii Nauk*, volume 136, pages 1034–1037. Russian Academy of Sciences, 1961.

[55] M. Tang, N. Vauchelet, I. Cheddadi, I. Vignon-Clementel, D. Drasdo, and B. Perthame. Composite waves for a cell population system modeling tumor growth and invasion. In *Partial Differential Equations: Theory, Control and Approximation*, pages 401–429. Springer, 2014.
[56] N. Vauchelet and E. Zatorska. Incompressible limit of the Navier-Stokes model with a growth term. *Nonlinear Anal.*, 163:34–59, 2017.

[57] J. L. Vazquez. The mesa problem for the fractional porous medium equation. *Interfaces and Free Boundaries*, 17(2):263–289, 2015.