---
abstract: |
    Both compressible and incompressible porous medium models are used in the literature to describe the mechanical properties of living tissues. These two classes of models can be related using a stiff pressure law. In the incompressible limit, the compressible model generates a free boundary problem of Hele-Shaw type where incompressibility holds in the saturated phase.

    Here we consider the case with a nutrient. Then, a badly coupled system of equations describes the cell density number and the nutrient concentration. For that reason, the derivation of the free boundary (incompressible) limit was an open problem, in particular a difficulty is to establish the so-called *complementarity relation* which allows to recover the pressure using an elliptic equation. To establish the limit, we use two new ideas. The first idea, also used recently for related problems, is to extend the usual Aronson-Bénilan estimates in $L^\infty$ to an $L^2$ setting. The second idea is to derive a sharp uniform $L^4$ estimates on the pressure gradient, independently of space dimension.
author:
- 'Noemi David[^1]'
- 'Benoît Perthame[^2]'
bibliography:
- 'biblio.bib'
title: Free boundary limit of tumor growth model with nutrient
---

.7cm

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

Introduction {#introduction .unnumbered}
============

We consider a compressible mechanical model of tumor growth, where the cell motion is driven by the pressure gradient according to Darcy’s law. The cell proliferation is governed by a bio-mechanical form of contact inhibition, that prevents cell division when the total cell density exceeds a critical threshold. The evolution of the cell population density $n\geq 0$ and the concentration of nutrients $c\geq 0$ are described by the following type of system $$\label{sys}
    \begin{cases}
    \partial_t n - \operatorname{div}(n \nabla p) = n G(p, c),\qquad \qquad x \in {\mathbb{R}}^d, \;  t\geq 0,\\[5pt]
    \partial_t c - \Delta c + n H(c) = (c_B-c)K(p),\\[5pt]
    c(x,t) \rightarrow c_B \quad \text{for }  x \rightarrow \infty.
    \end{cases}$$ The pressure within the tissue is denoted by $p$, and in the compressible setting, we use for simplicity the following law of state $$\label{pressure law}
    p= n^\gamma, \qquad  \gamma > 1.$$ The reaction term $G(p,c)$ is the cell division rate and the lowest value of pressure that prevents cell division is called *homeostatic pressure*, and we denote it by $p_H$. The concentration $c_B>0$ is the level of nutrients at the source, namely the network of blood vessels. Here, we consider the vascular phase of tumor growth, after *angiogenesis* has occurred, therefore the vasculature is present both outside and inside the tumor. The term $K\geq 0$ is the rate of nutrient release, which decreases with respect to the pressure. As clinical observation have shown, the mechanical stress generated by the cells shrinks the vessels inside the tumor and effects the blood flow and, by consequence, the nutrients delivery, see [@MMACCL] for further details. Finally, the term $H\geq 0$ is an increasing function of $c$ and represents the consumption rate of the nutrient by the tumor cells.

The specific form of the reaction term in the equation on $c$ is not fully relevant for our analysis, and we only need the possibility to derive some generic a priori estimates, mostly in $L^2$. Our study covers, for example, the terms in [@PQV] where the authors take $H=H(p,c)$, $K=0$ and those in [@PTV] where $K=\mathds{1}_{\{n=0\}}$, since the authors are considering the avascular phase of tumor growth. For our study, only some general conditions are needed, which are detailed in the next sections.

#### Motivation and previous works. {#motivation-and-previous-works. .unnumbered}

Models of tumor growth, including , possibly with more biological relevance, have been widely used recently. Several surveys are available, as [@RCM]. Numerical schemes for the model at hand, with AP property (asymptotic preserving), have been proposed in [@LTWZ].

Mechanical models of tumor growth are focused on the effect of the internal pressure which governs the dynamics of the cell population density. This kind of description was initiated in [@Green] by Greenspan and further developed by Byrne and Chaplain, [@BC96], Friedman, [@Fri], and Lowengrub et al., [@LFJCMWC], among the others. The leading assumption is that the birth of a cell generates a mechanical stress on the surrounding cells which start to move under a gradient of pressure. By consequence, the motion of the cells is usually described by Darcy’s law $$\label{Darcy}
    \vec{v}=-\nabla p,$$ which relates the velocity to the pressure gradient. This type of models have been extensively used to describe the early stage of tumor growth, the so-called *avascular phase*, see for example [@BCGR; @byrneChaplain96; @SCh]. Models of tumor growth that include the effect of viscosity, [@PV; @DeSc; @RJPJ], or more than one species of tissue cells, [@CFSS; @LLP], are also well-developed. For a comprehensive review on this topic we refer the reader to [@Fri; @LFJCMWC; @PT; @RSCB].

The equation for the density in the system is based on the continuous mechanical model presented in [@BD], in which the dynamics of tumor growth are governed by competition for space and contact inhibition. The *homeostatic pressure* is determined by the maximum level of stress that the cells can tolerate, see [@BD] for further details on the individual-based model that leads to the continuous one.

As explained above, this type of models are usually referred to as *compressible*, since they relate the density and the pressure through a compressible constitutive law, in a fluid mechanical view. A second class of models commonly used to describe tumor growth are free boundary problems, [@frihie]. They are also called geometric or *incompressible* models and describe the tumor as a moving domain where the density is constant. Free boundary problems arise also from the theory of mixture applied to tumor growth, [@BKMP2; @BP].

Building a link between these two classes of models has attracted the attention of many researchers in recent years. This result has first been achieved in [@PQV] for a purely mechanical model, passing to the so-called *incompressible limit*, as the pressure becomes stiff. Later, it has been studied for a lot of models, which included viscosity [@PV; @DeSc], different laws of state [@DHV] and more than one species of cells [@BPPS]. In each case the limit model turns out to be a free boundary model of Hele-Shaw type.

Our goal is to study the limit $\gamma \to \infty$ in the law of state , and prove that the limit solution satisfies a free boundary problem. It has been proved in [@PQV] that (the norms are specified in the next section and we now use the notation $n_\gamma, \,p_\gamma, \, {c_\gamma}$ in place of $n, \, p, \, c$ to indicate the dependency upon $\gamma$) $$n_\gamma \to n_\infty, \qquad p_\gamma \to p_\infty, \qquad c_\gamma \to c_\infty,$$ and the limits satisfy the system $$\label{syslim}
    \begin{cases}
    \partial_t n_\infty - \operatorname{div}(n_\infty \nabla p_\infty) = n_\infty G(p_\infty, c_\infty),\qquad \qquad x \in {\mathbb{R}}^d, \;  t\geq 0,\\[5pt]
    \partial_t c_\infty - \Delta c_\infty + n_\infty H(c_\infty) = (c_B-c_\infty)K(p_\infty),\\[5pt]
    c_\infty(x,t) \rightarrow c_B \quad\text{ for } x \rightarrow \infty,
    \end{cases}$$ with a *graph* relation between $p_\infty$ and $n_\infty$ given by $$\label{pressurelilmit}
    0\leq n_\infty \leq 1, \qquad  p_\infty (n_\infty-1)=0  .$$ A remarkable result is the uniqueness of the weak solutions of this system.

However, it was left open in [@PQV] to establish the so-called *complementarity condition*, which reads (in the sense of distributions) $$\label{crd}
    p_\infty \big(\Delta p_\infty + G(p_\infty, c_\infty) \big)=0 \quad\text{ in }\quad {\mathcal{D}}'({\mathbb{R}}^d\times (0,\infty)),$$ which follows formally from the equation on $n$ written for the pressure, namely $$\label{pressure}
    \partial_t {p_\gamma}= \gamma {p_\gamma}\big(\Delta {p_\gamma}+ G({p_\gamma},{c_\gamma}) \big) +|\nabla {p_\gamma}|^2.$$ The complementarity condition is fundamental because it relates the weak solutions defined by the equations and to the geometric form of the Hele-Shaw problem, where the set ${\cal O}(t):=\{x;\, p(x,t) >0\}$ evolves with the speed determined by the normal component of $\nabla p_\infty$. The limit pressure is a solution to the elliptic equation with Dirichlet boundary conditions $$-\Delta p_\infty = G(p_\infty , c_\infty ) \qquad \text{in }\; {\cal O}(t)=\{x;\, p_\infty(x,t) >0\}.$$ The Hele-Shaw problem is a widely studied free boundary model. Although we are only interested in the weak formulation, the regularity of the boundary is also a challenging issue, see [@SalsaCaffa; @HS; @MePeQu].

#### Difficulties and strategies. {#difficulties-and-strategies. .unnumbered}

To handle this problem, we make use of two new estimates which hold because the cell population density satisfies a Porous Medium Equation, which reads $$\label{pme}
    \partial_t n_\gamma -\frac{\gamma}{\gamma +1}\Delta n_\gamma^{\gamma+1} = n_\gamma G(p_\gamma,c_\gamma).$$ $\bullet$ The first estimate results from the famous Aronson-Bénilan (AB in short) inequalities for the porous media, [@AB; @CrPi1982], which have been extended in various contexts (see [@VaVi] for another example). It was used in the purely mechanical case, [@PQV], and it gives the lower bound $\Delta p_\gamma(t) +G( p_\gamma(t)) \geq - C/\gamma t$, with $C$ positive constant. Here, unlike in the case without nutrients, it cannot hold. In fact, as shown in [@PTV], where a semi-explicit travelling wave solution was found, there exists a region where $p$ is constantly equal to zero and $G$ is negative.

Therefore, we show a weaker, but still sufficient, condition $$\int_0^T \int_{{\mathbb{R}}^d} |\min (0, \Delta p_\gamma )|^3 \leq C(T).$$ This is proved by working in $L^2$, rather than with a sub-solution, as it has been recently initiated in [@GPS; @BPPS]. This method has the advantage to be compatible with the $L^2$ estimates on $c_\gamma$ and its derivatives. We recall here that $\Delta p_\infty$ is a bounded measure due to the free boundary of the set ${\cal O}(t)$ where the pressure is positive.\
$\bullet$ The second new estimate is an $L^4$ bound on $\nabla p_\gamma$, independent of the dimension $d$. In the simple case, where $G$ depends only on $p$, it results from the kinetic energy relation combined to the AB inequality in $L^\infty$, which is wrong here. We have a new and more general way to derive it, independently of the AB inequality.

#### Plan of the paper. {#plan-of-the-paper. .unnumbered}

The paper is organized as follows. The next section is devoted to explain the notations and assumptions and to state the main result of the paper, namely that the complementarity condition holds. The rest of the paper is dedicated to prove this result. We begin in Section \[Sec:est1\] presenting standard bounds which are useful for deriving the main new estimates that are stated and proved in Section \[Sec:est2\]. Finally, in Section \[Sec:4\] we give the proof of the complementarity relation.

Notations, assumptions and main result {#sec:NAMR}
======================================

#### Notations. {#notations. .unnumbered}

We denote $Q={\mathbb{R}}^d \times (0,\infty)$, and for $T>0$ we set $Q_T={\mathbb{R}}^d\times (0,T)$. We frequently use the abbreviation form $n(t):=n(x,t),\, p(t):=p(x,t),\, c(t):=c(x,t)$. We denote $$\operatorname{sign}_+{\{w\}}=\mathds{1}_{\{w>0\}}\quad \text{ and } \quad \operatorname{sign}_-{\{w\}}=-\mathds{1}_{\{w<0\}}.$$ We also define the positive and negative part of $w$ as follows $$|w|_+ :=
    \begin{cases}
    w, &\text{ for } w>0,\\
    0, &\text{ for } w\leq 0,
    \end{cases}
    \quad  \text{ and }\quad
    |w|_- :=
    \begin{cases}
    -w, &\text{ for } w<0,\\
    0, &\text{ for } w\geq 0.
    \end{cases}$$

#### Assumptions. {#assumptions. .unnumbered}

Considering the growth/reaction terms, the functions $G$, $H$ and $K$ are assumed to be smooth and we make the following assumption. There exist positive constants $\beta$, $p_H$, $p_B$ (reference pressure of blood vessels) such that $$\begin{aligned}
\label{ReTe}
    \partial_p G &<-\beta, & \partial_c G&\geq 0, & G(p,c_B)\leq 0,& \text{ for } p\geq p_H,\\
    K'(p) &\leq 0,    & 0 \leq K(p) &\leq 1,  & \quad K(p)=0,& \text{ for } p\geq p_B,\\
    H'(c) &\geq 0,   & 0 \leq H(c) &,  & H(0)=0 ,& 
    \end{aligned}$$ Furthermore, for a given pressure $p$, $G(p,c)<0$ for $c$ small enough. This assumption indicates that the tumor cells die by *necrosis* when the concentration of nutrients is below a survival threshold.

Some standard choices for the reaction terms are $$\begin{aligned}
    G(p,c)&=g(p)(c+c_1)-c_2, & H(c)&=c, & K(p)&=\left|1-\frac{p}{p_B}\right|_+,
    \end{aligned}$$ where $c_{1,2}$ are positive constants and $g$ is a decreasing function of $p$, see [@CBCB; @MMACCL; @PQV].

#### Initial data. {#initial-data. .unnumbered}

The system is completed with initial data. We assume that for some $n^0, c^0$, the initial data $n^0_\gamma, {c_\gamma}^0$ satisfy $$\begin{aligned}
    \label{inf n}  &0\leq n^0_\gamma \leq n_H:=p_H^{1/\gamma}, &\|n^0_\gamma-n^0\|_{L^1({\mathbb{R}}^d)}\xrightarrow[\gamma    \rightarrow \infty]{}0, \qquad \qquad &n^0 \in L_+^1({\mathbb{R}}^d), \\
    \label{inf c}  &0\leq {c_\gamma}^0\leq c_B,  &\|c^0_\gamma-c^0\|_{L^1({\mathbb{R}}^d)}\xrightarrow[\gamma  \rightarrow \infty]{}0,\qquad \qquad & c^0-c_B \in L_+^1({\mathbb{R}}^d).
    \end{aligned}$$ We also assume that there is a positive constant $C$ such that $$\begin{aligned}
    \label{2 deltap}         \|\nabla p^0_\gamma\|_{L^2({\mathbb{R}}^d)} + \|\Delta p^0_\gamma\|_{L^2({\mathbb{R}}^d)} & \leq C,\\
    \label{timezero} \|(\partial_t {n_\gamma})^0\|_{L^1({\mathbb{R}}^d)}+\|(\partial_t {c_\gamma})^0\|_{L^1({\mathbb{R}}^d)}&\leq C,\\
    \label{2 gradc}        \|\nabla {c_\gamma}^0\|_{L^2({\mathbb{R}}^d)} & \leq C.
    \end{aligned}$$ Set these conditions on the initial data, we give the definition of weak solution of the system as follows.

Given $T>0$, a weak solution of the system is a triple $({n_\gamma},{p_\gamma},{c_\gamma})$ such that, $${n_\gamma}, {p_\gamma}, {c_\gamma}\in L^{\infty}((0,T),L^p({\mathbb{R}}^d)) \quad \forall p\geq 1, \qquad  \nabla {c_\gamma}, \; \nabla {p_\gamma}\in  L^{2}({\mathbb{R}}^d \times (0,T)),$$ and for all $\varphi\in C^1_{\text{comp}}({\mathbb{R}}^d \times [0,T))$, $$\begin{aligned}
        &\int_0^T \int_{{\mathbb{R}}^d}\left(-{n_\gamma}\partial_t \varphi + {n_\gamma}\nabla{p_\gamma}\nabla\varphi-{n_\gamma}G({p_\gamma},{c_\gamma})\varphi\right)=\int_{{\mathbb{R}}^d} {n_\gamma}^0 \varphi(0),\\
        &\int_0^T \int_{{\mathbb{R}}^d}\left(-c_\gamma \partial_t \varphi + \nabla c_\gamma \nabla\varphi +{n_\gamma}H({c_\gamma})\varphi - (c_B-c) K(p)\varphi\right)=\int_{{\mathbb{R}}^d} c_\gamma^0 \varphi(0).
        \end{aligned}$$

From [@V] we know that a weak solution exists for all $T>0$.

#### Compact support. {#compact-support. .unnumbered}

Because our arguments rely on technical calculations, we first simplify the setting assuming that there exists a smooth bounded open domain $\Omega_0\subset {\mathbb{R}}^d$, independent of $\gamma$, such that for all $\gamma > 1$ $$\text{supp}({n_\gamma}^0)\subset \Omega_0.$$ Unlike the solutions of the heat equation, the PME’s solutions have a finite speed of propagation, see [@V]. This means that, for all $T>0$, there exists a smooth bounded open domain $\Omega_T$ independent of $\gamma$ such that $$\
    \text{supp}({n_\gamma}(t)) \subset \Omega_T, \qquad \forall t  \in [0,T],$$ see Appendix \[appA\] for the proof. From now on, we consider a solution $({n_\gamma},{p_\gamma})$ with compact support for all $\gamma > 1$. In the Appendix \[app\], we show how to extend the result to more general solutions.

#### Main result. {#main-result. .unnumbered}

We now state the main result of the paper, namely the weak formulation of the complementarity relation.

\[CR1\] With all the previous assumptions, the limit pressure $p_\infty$ satisfies the relation , that means, for all test functions $\zeta \in {\mathcal{D}}(Q)$ , we have $$\iint_{Q}  \left( -|\nabla p_\infty|^2\zeta - p_\infty\nabla p_\infty\nabla\zeta + p_\infty G(p_\infty,c_\infty)\zeta\right)=0.$$ Furthermore the following estimates hold uniformly in $\gamma$ $$\int_0^T \int_{\Omega_T} | \Delta {p_\gamma}+ G({p_\gamma},{c_\gamma})|_-^3 \leq C(T), \qquad  \qquad  \int_0^T\int_{\Omega_T} |\nabla {p_\gamma}|^4 \leq C(T).$$

Preliminary Estimates {#Sec:est1}
=====================

Let $({n_\gamma}, {p_\gamma}, {c_\gamma})$ be a weak solution to the system . We recall some standard preliminary bounds on ${n_\gamma}, {p_\gamma}, {c_\gamma}$ and their derivatives, gathered in the following Proposition.

\[prop1\] Given $({n_\gamma}, {p_\gamma}, {c_\gamma})$ a weak solution of the system for $\gamma > 1$, and $T>0$, there exists a constant $C(T)$, independent of $\gamma$, such that for all $0\leq t\leq T$ $$\begin{aligned}
        &0\leq {n_\gamma}\leq n_H, & &0\leq {p_\gamma}\leq p_H, & &0\leq {c_\gamma}\leq c_B,\\
        \label{L1bds}   
        &\|{n_\gamma}(t)\|_{L^1({\mathbb{R}}^d)}\leq C(T), & &\|{p_\gamma}(t) \|_{L^1({\mathbb{R}}^d)}\leq C(T), & &\|{c_\gamma}(t)-c_B\|_{L^1({\mathbb{R}}^d)} \leq C(T), \\
        \label{2 c}  &\|\nabla{c_\gamma}(t)\|_{L^2({\mathbb{R}}^d)} \leq C(T), &  &\|\Delta{c_\gamma}\|_{L^2(Q_T)} \leq C(T), &  &\|\partial_t {c_\gamma}\|_{L^2(Q_T)} \leq C(T),\\
        \label{timeder}
        &\|\partial_t {n_\gamma}\|_{L^1(Q_T)} \leq C(T),&  &\|\partial_t {p_\gamma}\|_{L^1(Q_T)} \leq C(T),&  &\|\partial_t {c_\gamma}\|_{L^1(Q_T)} \leq C(T), \\
        \label{c4p2}   &\|\nabla{c_\gamma}\|_{L^4(Q_T)} \leq C(T), & &\|\nabla{p_\gamma}\|_{L^2(Q_T)} \leq C(T).
        \end{aligned}$$

For the sake of completeness, we now recall the derivation of these bounds.

#### $L^\infty$ bounds for ${n_\gamma}, {p_\gamma}, {c_\gamma}$. {#linfty-bounds-for-n_gamma-p_gamma-c_gamma. .unnumbered}

The $L^\infty$ bounds are just consequences of our assumptions on $G$ using the comparison principle.

#### $L^1$ bounds on ${n_\gamma}, {p_\gamma}, {c_\gamma}$. {#l1-bounds-on-n_gamma-p_gamma-c_gamma. .unnumbered}

These are also standard estimates, noting that $$\| p(t) \|_{L^1({\mathbb{R}}^d)} = \| n(t) p(t)^\frac{\gamma- 1}{\gamma}  \|_{L^1({\mathbb{R}}^d)} \leq p_H^\frac{\gamma- 1}{\gamma} \| n(t) \|_{L^1({\mathbb{R}}^d)}.$$

#### $L^2$ bounds for the derivatives of ${c_\gamma}$. {#l2-bounds-for-the-derivatives-of-c_gamma. .unnumbered}

We now prove the $L^2$ bounds for $\nabla {c_\gamma}, \Delta {c_\gamma}$ and $\partial_t {c_\gamma}$. We multiply the equation for ${c_\gamma}$ by $-\Delta {c_\gamma}$ and we integrate in space and time $$- \int_0^t\int_{{\mathbb{R}}^d} \partial_t {c_\gamma}\Delta {c_\gamma}+\int_0^t\int_{{\mathbb{R}}^d} |\Delta {c_\gamma}|^2 = \int_0^t\int_{{\mathbb{R}}^d}({n_\gamma}H({c_\gamma})-(c_B-c)K({p_\gamma}))\Delta {c_\gamma}.$$ Integrating by parts and using Young’s inequality we obtain $$\int_0^t \int_{{\mathbb{R}}^d} \partial_t(\nabla {c_\gamma}) \nabla {c_\gamma}+\int_0^t\int_{{\mathbb{R}}^d} |\Delta {c_\gamma}|^2 \leq \int_0^t\int_{{\mathbb{R}}^d}\frac{|{n_\gamma}H({c_\gamma})-(c_B-c)K({p_\gamma})|^2}{2}+ \int_0^t\int_{{\mathbb{R}}^d}\frac{|\Delta {c_\gamma}|^2}{2}.$$ Hence, we have $$\frac{1}{2}\int_{{\mathbb{R}}^d} |\nabla {c_\gamma}(t)|^2  +\frac{1}{2}\int_0^t\int_{{\mathbb{R}}^d} |\Delta {c_\gamma}|^2 \leq C \int_0^t \left(\|{n_\gamma}(s)\|^2_{L^1({\mathbb{R}}^d)}+\|{c_\gamma}(s)-c_B\|^2_{L^1({\mathbb{R}}^d)}\right)ds+ \frac{1}{2}\|\nabla {c_\gamma}^0\|^{2}_{L^2({{\mathbb{R}}^d})},$$ where $C$ is a positive constant depending on $n_H, c_B$ and the $L^\infty$ norms of $H$ and $K$.

Finally, using the $L^1$ bounds , we obtain $$\begin{aligned}
    \int_{{\mathbb{R}}^d} |\nabla {c_\gamma}(t)|^2  +\int_0^t\int_{{\mathbb{R}}^d} |\Delta {c_\gamma}|^2 \leq C(T) + \|\nabla {c_\gamma}^0\|^{2}_{L^2({{\mathbb{R}}^d})} ,
    \end{aligned}$$ for $0<t\leq T$, and thanks to we conclude the proof of the first and second estimates in .

At last, considering the equation for ${c_\gamma}$ $$\partial_t {c_\gamma}= \Delta {c_\gamma}-{n_\gamma}H({c_\gamma}) + (c_B -{c_\gamma})K({p_\gamma}),$$ and using the previous bounds on ${n_\gamma}, {c_\gamma}$ and $\Delta {c_\gamma}$ we conclude that $\partial_t {c_\gamma}\in L^2(Q_T)$.

#### $L^1$ bounds for the time derivatives of ${n_\gamma}$ and ${p_\gamma}$. {#l1-bounds-for-the-time-derivatives-of-n_gamma-and-p_gamma. .unnumbered}

We differentiate the equation for ${n_\gamma}$ and we multiply it by $\operatorname{sign}{\{\partial_t {n_\gamma}\}}$ $$\label{dndt}
    \partial_t |\partial_t {n_\gamma}|-\gamma\Delta 
    ({n_\gamma}^{\gamma}|\partial_t {n_\gamma}|)\leq |\partial_t {n_\gamma}|G 
    + {n_\gamma}\partial_p G |\partial_t {p_\gamma}|+
    {n_\gamma}\partial_c G \partial_t {c_\gamma}\operatorname{sign}\{\partial_t {n_\gamma}\}.$$ We integrate in space using the monotonicity of $G$ $$\frac{d}{dt}\|\partial_t {n_\gamma}(t)\|_{L^1({\mathbb{R}}^d)}\leq \|G\|_{L^\infty(Q_T)}\|\partial_t {n_\gamma}(t)\|_{L^1({\mathbb{R}}^d)}+\|\partial_c G\|_{L^\infty(Q_T)}\|{n_\gamma}(t)\|_{L^2({\mathbb{R}}^d)}\|\partial_t{c_\gamma}(t)\|_{L^2({\mathbb{R}}^d)}.$$ Thanks to and , Gronwall’s lemma gives $$\|\partial_t {n_\gamma}(t)\|_{L^1({\mathbb{R}}^d)}\leq C(T)  \|(\partial_t {n_\gamma})^0\|_{L^1({\mathbb{R}}^d)}\leq C(T),$$ where in the last inequality we used .

By integrating in $Q_t:={\mathbb{R}}^d\times (0,t)$ the equation , we obtain $$\|\partial_t {n_\gamma}(t)\|_{L^1({\mathbb{R}}^d)}+\min|\partial_p G|\iint_{Q_t}{n_\gamma}|\partial_t {p_\gamma}| \leq C(T),$$ thanks to and the $L^1$ bounds proved above. Then, for the time derivative of the pressure it holds $$\|\partial_t {p_\gamma}\|_{L^1(Q_T)}\leq \iint_{Q_T \cap \{{n_\gamma}\leq 1/2\}}\gamma {n_\gamma}^{\gamma-1}|\partial_t {n_\gamma}| + 2 \iint_{Q_T \cap \{{n_\gamma}\geq 1/2\}} {n_\gamma}|\partial_t {p_\gamma}|\leq C(T).$$ We differentiate the equation for ${c_\gamma}$ and multiply it by $\operatorname{sign}{\{\partial_t {c_\gamma}\}}$ $$\partial_t |\partial_t {c_\gamma}|-\Delta 
    (|\partial_t {c_\gamma}|)\leq -\partial_t {n_\gamma}H \operatorname{sign}\{\partial_t{c_\gamma}\} -{n_\gamma}H' |\partial_t{c_\gamma}|  -|\partial_t{c_\gamma}| K+(c_B-c)K' \partial_t {p_\gamma}\operatorname{sign}\{\partial_t{c_\gamma}\} .$$ Integrating in space we obtain $$\begin{aligned}
    \frac{d}{dt}\|\partial_t {c_\gamma}(t)\|_{L^1({\mathbb{R}}^d)}\leq &\,\|H\|_{L^\infty(Q_T)}\|\partial_t {n_\gamma}(t)\|_{L^1({\mathbb{R}}^d)}+ n_H \|H'\|_{L^\infty(Q_T)}\|\partial_t {c_\gamma}(t)\|_{L^1({\mathbb{R}}^d)}\\&+c_B\|K'\|_{L^\infty(Q_T)}\|\partial_t p(t)\|_{L^1({\mathbb{R}}^d)},
    \end{aligned}$$ and thanks to the previous bounds and Gronwall’s lemma we have $$\|\partial_t {c_\gamma}(t)\|_{L^1({\mathbb{R}}^d)}\leq C(T) \|(\partial_t {c_\gamma})^0\|_{L^1({\mathbb{R}}^d)}\leq C(T),$$ and this concludes the proof of .

#### $L^4$ bound for the gradient of ${c_\gamma}$. {#l4-bound-for-the-gradient-of-c_gamma. .unnumbered}

Now, we prove that the gradient of ${c_\gamma}$ is bounded in $L^4$. Integration by parts gives $$\begin{aligned}
    \int_{{\mathbb{R}}^d} |\nabla {c_\gamma}|^4 = - \int_{{\mathbb{R}}^d} {c_\gamma}\Delta {c_\gamma}|\nabla {c_\gamma}|^2 - \int_{{\mathbb{R}}^d} {c_\gamma}\nabla {c_\gamma}\cdot \nabla(|\nabla {c_\gamma}|^2).
    \end{aligned}$$ We use Young’s inequality on the first term of the RHS and we get $$\begin{aligned}
    \frac{1}{2} \int_{{\mathbb{R}}^d} |\nabla {c_\gamma}|^4 \leq \frac{1}{2} \int_{{\mathbb{R}}^d} {c_\gamma}^2 |\Delta {c_\gamma}|^2 - \int_{{\mathbb{R}}^d} {c_\gamma}\nabla {c_\gamma}\cdot \nabla(|\nabla {c_\gamma}|^2).
    \end{aligned}$$ We write the last term as $$\begin{aligned}
    - \int_{{\mathbb{R}}^d} {c_\gamma}\nabla {c_\gamma}\cdot \nabla(|\nabla {c_\gamma}|^2) &= -2 \sum_{i,j} \int_{{\mathbb{R}}^d}{c_\gamma}\, \partial_i {c_\gamma}\, \partial_j {c_\gamma}\, \partial^2_{i,j} {c_\gamma}\\
    &\leq \frac{1}{4} \int_{{\mathbb{R}}^d} |\nabla {c_\gamma}|^4 + 4c_B^2\int_{{\mathbb{R}}^d}  \sum_{i,j}(\partial^2_{i,j}{c_\gamma})^2
    \\
    &= \frac{1}{4} \int_{{\mathbb{R}}^d} |\nabla {c_\gamma}|^4 + 4c_B^2\int_{{\mathbb{R}}^d}   |\Delta {c_\gamma}|^2.  
    \end{aligned}$$ Thus, we have $$\frac{1}{4} \int_{{\mathbb{R}}^d} |\nabla {c_\gamma}|^4 \leq \left(\frac{1}{2}+4\right)  c_B^2\int_{{\mathbb{R}}^d}   |\Delta {c_\gamma}|^2 .$$ and the $L^4$ estimate is proved.

#### $L^2$ bound for the pressure gradient. {#l2-bound-for-the-pressure-gradient. .unnumbered}

Since the pressure satisfies the equation , integrating it in space we get $$\begin{aligned}
    \frac{d}{dt}\int_{{\mathbb{R}}^d} {p_\gamma}(t)= -\gamma\int_{{\mathbb{R}}^d} |\nabla {p_\gamma}(t)|^2 + \gamma \int_{{\mathbb{R}}^d} {p_\gamma}(t)G({p_\gamma}(t),{c_\gamma}(t)) +\int_{{\mathbb{R}}^d} |\nabla {p_\gamma}(t)|^2.
    \end{aligned}$$ Then, we integrate in time $$\begin{aligned}
    (\gamma -1)\int_0^T\int_{{\mathbb{R}}^d} |\nabla {p_\gamma}|^2&=  \|{p_\gamma}(0)\|_{L^1({\mathbb{R}}^d)}-\| {p_\gamma}(T)\|_{L^1({\mathbb{R}}^d)} +\gamma \int_0^T\int_{{\mathbb{R}}^d} {p_\gamma}G({p_\gamma},{c_\gamma}),\\
    (\gamma -1)\int_0^T\int_{{\mathbb{R}}^d} |\nabla {p_\gamma}|^2&\leq C_0 + \gamma C(T),
    \end{aligned}$$ and this gives, since $\gamma >1$, $$\int_0^T\int_{{\mathbb{R}}^d} |\nabla {p_\gamma}|^2\leq  C(T).$$

Stronger a priori estimates on $p_\gamma$ {#Sec:est2}
=========================================

To establish the complementarity condition is equivalent to prove the strong compactness of $|\nabla p_\gamma |^2$. One step towards this goal is to prove compactness in space using the classical AB estimate, [@AB; @CrPi1982]. Here, major difficulties arise. As explained in the Introduction, since the reaction term can change sign the usual Aronson-Bénilan lower bound cannot hold true, see [@PQV; @PTV]. Moreover, we cannot apply the comparison principle because of the bad coupling in the system . Since the $L^\infty$ bound from below in the AB estimate is missing, we prove an $L^3$ version, adapting the method presented in [@GPS]. Then, we show that the gradient of the pressure is bounded in $L^4(Q_T)$, which gives the compactness needed to pass to the limit.

Our first goal is to prove the AB estimate on the functional $$\label{w}
    w:= \Delta {p_\gamma}+ G({p_\gamma},{c_\gamma}),$$ which is a variation of the Laplacian in order to take into account the source term, at the same order of $\Delta p_\gamma$, in equation .

\[Thm1\] With the assumptions of Section \[sec:NAMR\] and with $\gamma >\max(1,2-\frac 4 d)$, for all $T>0$ there is a constant $C(T)$ depending on $T$ and the previous bounds and independent of $\gamma$ such that $$\label{AB}
        \int_0^T \int_{\Omega_T} |w|_-^3 \leq C(T), \quad 
        \int_0^T \int_{\Omega_T}  |\nabla {p_\gamma}|^2|w|_-\leq C(T),\quad
        \int_0^T \int_{{\mathbb{R}}^d} |\Delta {p_\gamma}| \leq C(T).$$

Let us point out that because the free boundary is where $p_\infty$ vanishes, it is important that $w$ itself is controlled and not merely $p w$ as in the next estimate.

\[L4\] With the same assumptions as before, given $T>0$, it holds $$\label{t4}
        (\gamma -1)    \int_0^T\int_{\Omega_T} {p_\gamma}|\Delta {p_\gamma}+G|^2 + \int_0^T\int_{\Omega_T} {p_\gamma}\sum_{i,j} (\partial^2_{i,j}{p_\gamma})^2 \leq C(T),$$ $$\label{grad p 4}
        \int_0^T\int_{\Omega_T} |\nabla {p_\gamma}|^4 \leq C(T),$$ where $C$ depends on $T$ and previous bounds and is independent of $\gamma$.

We recall that in the model independent of $c$, [@PQV], the AB estimate is much stronger and gives $\Delta {p_\gamma}(t) + G({p_\gamma}(t)) \geq -\frac{1}{\gamma t}$, and the major difficulty is the control of $\Delta {p_\gamma}$ which is provided by Theorem \[AB\]. As proved in [@MePeQu], the $L^4$ estimate follows from the total energy control when $G=G(p)$, but this uses the strong form of the AB estimate. Therefore, we have invented another proof, which is reminiscent of the energy control, but uses a different treatment of the ’dissipation’ terms.

#### Proof of Theorem \[Thm1\]. {#proof-of-theoremthm1. .unnumbered}

For the sake of simplicity we forget the index $\gamma$ in this proof. We compute the time derivative of $w$ and obtain $$\partial_t w =\Delta(|\nabla p|^2)+\gamma \Delta(p w)+ \partial_p G (|\nabla p|^2+\gamma p w)+\partial_c G\, \partial_t c.$$ The first term is $$\Delta (|\nabla p|^2)= 2\sum_{i,j} (\partial^2_{i,j}p)^2  + 2\nabla p \cdot\nabla(\Delta p) \geq \frac{2}{d}(\Delta p)^2 + 2\nabla p \cdot\nabla(\Delta p).$$ By definition of $w$ we have $$\begin{aligned}
    2\nabla p \cdot \nabla(\Delta p)=2\nabla p \cdot \nabla (w-G)=2\nabla p\cdot\nabla w -2\partial_p G |\nabla p|^2 -2 \partial_c G\nabla p \cdot \nabla c.
    \end{aligned}$$ Hence, the time derivative satisfies $$\begin{aligned}
\label{w_}
    \partial_t w\geq &\frac{2}{d}(w-G)^2+2\nabla p\cdot\nabla w -\partial_p G |\nabla p|^2 -2 \partial_c G\nabla p \cdot \nabla c\\\nonumber
    &+\gamma \Delta(p w)
    +\gamma p w \,\partial_p G +\partial_c G \,\partial_t c.
    \end{aligned}$$ Multiplying by $-|w|_-$, we obtain $$\begin{aligned}
    -\partial_t w\, |w|_- \leq &-\frac{2}{d}|w|_-^3-\frac{4}{d} G|w|^2_--\frac{2}{d}G^2|w|_- +  \nabla p\cdot \nabla |w|_-^2 +\partial_p G |\nabla p|^2|w|_- \\
    &+2 \partial_c G\nabla p \cdot \nabla c |w|_-
    +\gamma \Delta(p |w|_-)|w|_-+\gamma p\, \partial_p G  |w|_-^2 - \partial_c G\, \partial_t c |w|_-.
    \end{aligned}$$ Hence, using the fact that $\partial_p G< -\beta$ from , we integrate in space to obtain $$\begin{aligned}
    \frac{d}{{dt}}\int_{\Omega_T}\frac{|w|_-^2}{2}\leq &-\frac{2}{d}\int_{\Omega_T}|w|_-^3 -\frac{2}{d}\int_{\Omega_T}   G^2 |w|_- -\beta \int_{\Omega_T} |\nabla p|^2|w|_- \\&-\frac{4}{d} \int_{\Omega_T}   G |w|_-^2 
    +\underbrace{\int_{\Omega_T}\left[ \nabla p \cdot \nabla |w|_-^2 + \gamma \Delta(p |w|_-)|w|_- \right] }_{{\mathcal{A}}} 
    \\&\underbrace{-\int_{\Omega_T} \partial_c G\,\partial_t c |w|_- }_{{\mathcal{B}}} +\underbrace{2\int_{\Omega_T} \partial_c G\nabla p \cdot \nabla c |w|_-}_{{\mathcal{C}}},
    \end{aligned}$$ where $C$ is a positive constant depending on $\|G\|_\infty$ and $d$. Now we proceed integrating by parts each term. $$\begin{aligned}
        \nonumber  {\mathcal{A}}&=-\int_{\Omega_T}\left[ \Delta p |w|_-^2  + \gamma  \nabla p \nabla |w|_- |w|_- + \gamma p |\nabla|w|_-|^2\right]
    \\
    \nonumber
    &= \int_{\Omega_T} |w|_-^3  +\int_{\Omega_T} G |w|_-^2+ \frac{\gamma}{2}\int_{\Omega_T} \Delta p |w|_-^2-\gamma \int_{\Omega_T}p |\nabla|w|_-|^2\\
    \label{rel1}
    &=\left(1-\frac{\gamma}{2}\right) \int_{\Omega_T}|w|_-^3 +\left(1-\frac{\gamma}{2}\right) \int_{\Omega_T}G|w|_-^2-\gamma\int_{\Omega_T}p|\nabla |w|_-|^2.
    \end{aligned}$$ Next, using and the Cauchy-Schwarz inequality, we obtain $${\mathcal{B}}\leq\, C \int_{\Omega_T}|w|_-^2+C.$$ Thanks to Young’s inequality and , we compute $$\begin{aligned}
    {\mathcal{C}}&\leq \frac \beta 2 \int_{\Omega_T}|\nabla p|^2  |w|_- + C \int_{\Omega_T}| \nabla c |^4 + C \int_{\Omega_T}| w|_-^2 \\
    &\leq \frac \beta 2 \int_{\Omega_T}|\nabla p|^2  |w|_- + C \int_{\Omega_T}| w|_-^2 +C.
    \end{aligned}$$ We may now come back to the control of $\frac{d}{dt}\int_{\Omega_T}\frac{|w|_-^2}{2}$. Gathering all the previous bounds, we get the following estimate $$\begin{aligned}
    \frac{d}{{dt}}\int_{\Omega_T}\frac{|w|_-^2}{2}&\leq  -\left(\frac{2}{d}-1+\frac{\gamma}{2}\right)\int_{\Omega_T}|w|_-^3 - \frac{\beta}{2} \int_{\Omega_T}|\nabla p|^2|w|_- +C(\gamma+1)\int_{\Omega_T} |w|_-^2 +C.
    \end{aligned}$$ Hence integrating in time we have $$\begin{aligned}
\label{time}
    \frac{\beta}{2} \int_0^T\int_{\Omega_T}|\nabla p|^2|w|_- +\left(\frac{2}{d}-1+\frac{\gamma}{2}\right) \int_0^T\int_{\Omega_T}|w|_-^3
    &\leq   C\, (\gamma +1) \int_0^T\int_{\Omega_T} |w|_-^2  +C(T)+ \int_{\Omega_T} \frac{|w^0|_-^2}{2}\\\nonumber
    &\leq  C\, (\gamma +1) \left(\int_0^T\int_{\Omega_T}|w|_-^3\right)^\frac{2}{3} +C(T),
    \end{aligned}$$ where we used the assumption and C represents different constants depending on $T$, $\left|\Omega(T) \right|$ and previous bounds. This is the place where we strongly use the compact support assumption.

At last, with our assumption that $\gamma$ is large enough, we obtain $$\begin{aligned}
    \int_0^T \int_{\Omega_T}|w|_-^3&\leq  C\left(\int_0^T\int_{\Omega_T}|w|_-^3\right)^\frac{2}{3}+C(T),
    \end{aligned}$$ and hence we have proved our main result, that is the first estimate of , $$\int_0^T\int_{\Omega_T}|w|_-^3\leq  C(T).$$ By consequence, we find $$\int_0^T\int_{\Omega_T}|\nabla p|^2|w|_- \leq C(T).$$ To prove the last estimate of , we argue as follows. Since $$\int_0^T \int_{\Omega_T} (\Delta p + G)\leq C(T),$$ we can also control the positive part of $w$ $$\int_0^T \int_{\Omega_T} |w|_+\leq C(T) +  \int_0^T \int_{\Omega_T} |w|_- \leq C(T)+ \left( \int_0^T \int_{\Omega_T} |w|_-^3\right)^{\frac 13}.$$ Thus it holds $$\int_0^T \int_{\Omega_T} |\Delta p + G| \leq C(T).$$ Hence, we finally obtain the $L^1$ estimate for the Laplacian of the pressure $$\int_0^T \int_{\Omega_T} |\Delta p| \leq C(T),$$ that concludes the proof of Theorem \[Thm1\].

#### Proof of Theorem \[L4\]. {#proof-of-theoreml4. .unnumbered}

We consider the equation for the pressure , we multiply it by $-(\Delta {p_\gamma}+ G({p_\gamma},{c_\gamma}))$ and integrate in space. We find successively $$-\int_{\Omega_T} \partial_t {p_\gamma}\Delta {p_\gamma}- \int_{\Omega_T} \partial_t {p_\gamma}\, G = -\gamma \int_{\Omega_T} {p_\gamma}|\Delta {p_\gamma}+ G|^2 - \int_{\Omega_T} |\nabla {p_\gamma}|^2 \Delta {p_\gamma}- \int_{\Omega_T} |\nabla {p_\gamma}|^2 G ,$$ $$\frac{d}{{dt}} \int_{\Omega_T} \frac{|\nabla {p_\gamma}|^2}{2} -\int_{\Omega_T} \partial_t {p_\gamma}G+\gamma \int_{\Omega_T} {p_\gamma}|\Delta {p_\gamma}+ G |^2 +\int_{\Omega_T} |\nabla {p_\gamma}|^2 \Delta {p_\gamma}\leq \|G\|_{L^\infty} \|\nabla {p_\gamma}(t)\|^2_{L^2}.$$ We integrate by parts the last term of the LHS and obtain $$\begin{aligned}
    \int_{\Omega_T} |\nabla {p_\gamma}|^2 \Delta {p_\gamma}&=\int_{\Omega_T} {p_\gamma}\Delta(|\nabla {p_\gamma}|^2)\\
    & = 2 \int_{\Omega_T} {p_\gamma}\nabla {p_\gamma}\cdot \nabla (\Delta {p_\gamma}) +2 \int_{\Omega_T}{p_\gamma}\sum_{i,j} (\partial_{i,j}^2{p_\gamma})^2\\
    &= -2 \int_{\Omega_T}{p_\gamma}|\Delta {p_\gamma}|^2 - 2\int_{\Omega_T} |\nabla {p_\gamma}|^2\Delta {p_\gamma}+2 \int_{\Omega_T}{p_\gamma}\sum_{i,j} (\partial^2_{i,j}{p_\gamma})^2.
    \end{aligned}$$ Hence, we conclude that $$\int_{\Omega_T} |\nabla {p_\gamma}|^2 \Delta {p_\gamma}= - \frac{2}{3} \int_{\Omega_T}{p_\gamma}|\Delta {p_\gamma}|^2 +\frac{2}{3} \int_{\Omega_T}{p_\gamma}\sum_{i,j} (\partial^2_{i,j}{p_\gamma})^2.$$ Thus, we have $$\label{diseq}
    \frac{d}{{dt}} \int_{\Omega_T} \frac{|\nabla {p_\gamma}|^2}{2} \underbrace{-\int_{\Omega_T} \partial_t {p_\gamma}\, G}_{I_1}+\underbrace{\gamma \int_{\Omega_T} {p_\gamma}|\Delta {p_\gamma}+ G|^2 -\frac{2}{3}\int_{\Omega_T} {p_\gamma}|\Delta {p_\gamma}|^2}_{I_2}  +\frac{2}{3} \int_{\Omega_T}{p_\gamma}\sum_{i,j} (\partial^2_{i,j}{p_\gamma})^2 \leq C(T).$$ We can define the function $\overline{G}=\overline{G}({p_\gamma},{c_\gamma})=\int_0^{{p_\gamma}} G(q,{c_\gamma}) dq$ and then $$\partial_t {p_\gamma}\,G({p_\gamma},{c_\gamma})=\partial_t \overline{G}({p_\gamma},{c_\gamma})-\partial_t {c_\gamma}\, \partial_c \overline{G}({p_\gamma},{c_\gamma}).$$ Using this relation the term $I_1$ can be written as $$I_1=-\int_{\Omega_T}\partial_t \overline{G}+\int_{\Omega_T}\partial_c\overline{G}\, \partial_t {c_\gamma}\geq -\int_{\Omega_T}\partial_t \overline{G}-  C$$ thanks to the $L^2$ bound on $\partial_t {c_\gamma}$ in  and because $|\partial_c\overline{G} | \leq C {p_\gamma}$. We can estimate the term $I_2$ from below as follows $$I_2\geq\left(\gamma - 1 \right)\int_{\Omega_T}{p_\gamma}|\Delta {p_\gamma}+G|^2- C \int_{\Omega_T}{p_\gamma}|G|^2.$$ Therefore $$\label{I1I2}
    I_1+I_2\geq \left(\gamma - 1  \right)\int_{\Omega_T} {p_\gamma}|\Delta {p_\gamma}+G|^2  -\int_{\Omega_T}\partial_t \overline{G} -C(T).$$ Combining and , we obtain $$\frac{d}{dt}\int_{\Omega_T} \left[ \frac{|\nabla {p_\gamma}|^2}{2} - \overline{G} \right] +(\gamma -1) \int_{\Omega_T}{p_\gamma}|\Delta {p_\gamma}+G|^2 + \frac{2}{3}\int_{\Omega_T} {p_\gamma}\sum_{i,j} (\partial^2_{i,j}{p_\gamma})^2 \leq C(T).$$ Finally, integrating in time, we obtain the estimate  and this proves the first step of Theorem \[L4\].

Furthermore, this bound also implies $$\label{pDelta}
    (\gamma -1) \int_0^T\int_{\Omega_T} {p_\gamma}|\Delta {p_\gamma}|^2 \leq C(T) .$$ Now, we compute the $L^4$ norm of the gradient of ${p_\gamma}$, as we did for the gradient of ${c_\gamma}$. $$\begin{aligned}
    \int_{\Omega_T} |\nabla {p_\gamma}|^4 = - \int_{\Omega_T} {p_\gamma}\Delta {p_\gamma}|\nabla {p_\gamma}|^2 - \int_{\Omega_T} {p_\gamma}\nabla {p_\gamma}\cdot \nabla(|\nabla {p_\gamma}|^2).
    \end{aligned}$$ Applying Young’s inequality to the first term, we obtain $$\begin{aligned}
    \frac{1}{2} \int_{\Omega_T} |\nabla {p_\gamma}|^4 \leq \frac{1}{2} \int_{\Omega_T} {p_\gamma}^2 |\Delta {p_\gamma}|^2 - 2 \sum_{i,j} \int_{\Omega_T}{p_\gamma}\, \partial_i {p_\gamma}\, \partial_j {p_\gamma}\, \partial^2_{i,j} {p_\gamma}.
    \end{aligned}$$ The last term can be upper bounded by $$\begin{aligned}
     \frac{1}{4} \int_{\Omega_T} |\nabla {p_\gamma}|^4 + 4\int_{\Omega_T} {p_\gamma}^2 \sum_{i,j}(\partial^2_{i,j}{p_\gamma})^2.  
    \end{aligned}$$ Therefore, we obtain $$\frac{1}{4} \int_{\Omega_T} |\nabla {p_\gamma}|^4 \leq \frac{1}{2} \int_{\Omega_T} {p_\gamma}^2 |\Delta {p_\gamma}|^2  + 4\int_{\Omega_T} {p_\gamma}^2 \sum_{i,j}(\partial^2_{i,j}{p_\gamma})^2.$$ Since $p\leq p_H$, by  and we conclude $$\begin{aligned}
    \int_0^T\int_{\Omega_T} |\nabla {p_\gamma}|^4 \leq& C(T),
    \end{aligned}$$ and this completes the proof of Theorem \[L4\].

Complementarity relation {#Sec:4}
========================

Thanks to the bounds provided by Theorem \[Thm1\] and Theorem \[L4\], we may obtain the desired compactness on the pressure gradient. This allows us to pass to the incompressible limit and prove the complementarity relation as we state it now.

\[CR\] With the assumptions of Theorem \[Thm1\], the complementarity condition  holds. More precisely, for all test functions $\zeta \in {\mathcal{D}}(Q)$, the limit pressure $p_\infty$ satisfies $$\label{compl}
        \iint_{Q}  \left( -|\nabla p_\infty|^2\zeta - p_\infty\nabla p_\infty\cdot\nabla\zeta + p_\infty G(p_\infty,c_\infty)\zeta\right)=0.$$

This result is related to the geometric form of the Hele-Shaw free boundary problem (while  is the weak form). It tells us that the limit solution satisfies $$\begin{cases}
    -\Delta p_\infty =  G(p_\infty, c_\infty) &\text{ in } {\cal O}(t) :=\{x;\, p_{\infty}(x,t) >0\} ,\\
    p_\infty=0 &\text{ on } \partial {\cal O}(t),
    \end{cases}$$ where, for every $t>0$, the set ${\cal O}(t)$ represents the region occupied by the tumor. Moreover, in the limit, the pressure and the cell population density satisfy the relation $$\label{relpn}
    p_\infty(1-n_\infty)=0.$$ In fact, we may expect that the set ${\cal O}(t)$ coincides a.e. with the set where $n_\infty=1$, hence the classification of *incompressible* model. See [@MePeQu] for the proof in the case without nutrient. It is not obvious to extend the result in the case at hand.

#### Proof of Theorem \[CR\]. {#proof-of-theoremcr. .unnumbered}

Thanks to the bounds in , and , ${p_\gamma}$ and ${c_\gamma}$ are locally compact and thus, after the extraction of sub-sequences, $$\begin{aligned}
    {p_\gamma}\rightarrow p_\infty \text{ strongly in } L^1(Q_T),\qquad 
    {c_\gamma}\rightarrow c_\infty \text{ strongly in } L^1(Q_T),
    \end{aligned}$$ when $\gamma \rightarrow \infty$, for all $T>0$. From Theorem \[L4\], we also recover the weak convergence of the gradient of the pressure, up to a sub-sequence, $$\nabla {p_\gamma}\rightharpoonup \nabla p_\infty  \text{ weakly in } L^4(Q_T).$$ From Theorem \[Thm1\], we know that $\Delta {p_\gamma}$ is bounded in $L^1$. Then, we have local compactness in space for the pressure gradient. To gain compactness in time we use the Aubin-Lions lemma. From the equation for the pressure , we have $$\partial_t (\nabla{p_\gamma})=\nabla[\gamma {p_\gamma}(\Delta  {p_\gamma}+G) + |\nabla {p_\gamma}|^2],$$ where the RHS is a sum of space derivatives of functions bounded in $L^1$. In fact, since by and $\partial_t {p_\gamma}$ and $|\nabla {p_\gamma}|^2$ are in $L^1$, from the term $\gamma {p_\gamma}(\Delta{p_\gamma}+G)$ is also bounded in $L^1$. Thus, we can extract a sub-sequence such that $$\nabla {p_\gamma}\rightarrow \nabla p_\infty \text{ strongly in } L^q(Q_T), \text{ for } 1 \leq q < \frac{d}{d-1}.$$ After extraction of a sub-sequence we obtain convergence almost everywhere for $\nabla {p_\gamma}$. Then, using the $L^4$ bound of Theorem \[L4\], we have $$\nabla {p_\gamma}\rightarrow \nabla p_\infty \text{ strongly in } L^q(Q_T), \text{ for } 1 \leq q < 4,$$ hence, in particular, also for $q=2$.

Let $\zeta \in {\mathcal{D}}(Q)$ be a test function. We consider the equation for ${p_\gamma}$ $$\partial_t {p_\gamma}=\gamma {p_\gamma}(\Delta {p_\gamma}+ G({p_\gamma}, {c_\gamma})) + |\nabla {p_\gamma}|^2,$$ we multiply it by $\zeta$ and we integrate in $Q$ $$\label{relation}
    - \frac{1}{\gamma} \iint_{Q} \left( {p_\gamma}\partial_t\zeta+|\nabla {p_\gamma}|^2\zeta \right)=\iint_{Q}  \left( -|\nabla {p_\gamma}|^2\zeta - {p_\gamma}\nabla {p_\gamma}\cdot\nabla\zeta + {p_\gamma}G({p_\gamma},{c_\gamma})\zeta\right)$$ Hence, passing to the limit for $\gamma \rightarrow \infty$ we obtain the complementarity relation $$\iint_{Q}  \left( -|\nabla p_\infty|^2\zeta - p_\infty\nabla p_\infty\cdot\nabla\zeta + p_\infty G(p_\infty,c_\infty)\zeta\right)=0.$$ This is equivalent to $$\iint_{Q}  p_\infty \left( \Delta p_\infty + G(p_\infty,c_\infty)\right) \zeta=0,$$ which means $$p_\infty \left( \Delta p_\infty + G(p_\infty,c_\infty)\right) =0, \;\text{ in } {\mathcal{D}}'(Q),$$ and the proof of Theorem \[CR\] is complete.

Acknowledgments {#acknowledgments .unnumbered}
===============

![image](euro.png){width="0.7cm"} 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).\
B.P. has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 740623).

Compact support property {#appA}
========================

We now give the proof of the finite speed of propagation property of the solutions of system . Our goal is to show that if the initial data satisfy $$\text{supp}(n_\gamma^0)\subset \Omega_0, \quad \forall \gamma >1,$$ with $\Omega_0$ independent of $\gamma$, then the solutions ${n_\gamma}(t)$, ${p_\gamma}(t)$ are compactly supported, uniformly in $\gamma$ and $t\in [0,T]$, for all $T>0$. This means that there exists a bounded open domain $\Omega_T$ independent of $\gamma$ such that $$\text{supp}(n_\gamma(t))\subset \Omega_T, \quad \forall\gamma>1,\, \forall t\in [0,T].$$ For every $\gamma>1$, the pressure ${p_\gamma}$ is a sub-solution to the equation $$\partial_t {p_\gamma}\leq |\nabla {p_\gamma}|^2+ \gamma {p_\gamma}(\Delta {p_\gamma}+ G(0,c_B)),$$ therefore, finding a super-solution with compact support we can control the supports of ${p_\gamma}$ and ${n_\gamma}$.

We consider the function $$\Pi (x,t) = G(0,c_B)\left|S(t)-\frac{|x|^2}{2}\right|_+,$$ where we choose the function $S$ such that it satisfies $$S'(t)\geq 2 G(0,c_B) S(t).$$ We compute the derivatives of $\Pi$ and we find $$\begin{aligned}
        &\partial_t \Pi (x,t)=G(0,c_B)S'(t)\mathds{1}_{\{S(t)\geq \frac{|x|^2}{2}\}   }   ,\\
        &\nabla \Pi(x,t)=-G(0,c_B) x \mathds{1}_{\{S(t)\geq \frac{|x|^2}{2}\}   }  ,\quad
        \Delta \Pi (x,t)\leq-d\, G(0,c_B)\mathds{1}_{\{S(t)\geq \frac{|x|^2}{2}\} }   .
        \end{aligned}$$ Therefore $\Pi$ satisfies $$\begin{aligned}
        \partial_t \Pi -|\nabla \Pi|^2-\gamma \Pi(\Delta \Pi +G(0,c_B))&\geq
        (  G(0,c_B)S'(t)  -G(0,c_B)^2 x^2)\mathds{1}_{\{S(t)\geq \frac{|x|^2}{2}\}}+\gamma \Pi G(0,c_B)(d+1)\\
        &\geq ( 2 G(0,c_B)^2 S(t)  -G(0,c_B)^2 x^2)\mathds{1}_{\{S(t)\geq \frac{|x|^2}{2}\} }\\
        &\geq 0.
        \end{aligned}$$ Hence, we have proved that for all $T>0$ $$\text{supp}({p_\gamma}(t))\subset \text{supp}(\Pi(t))\subset B_{T},\;\forall \gamma >1, \forall t\in[0,T],$$ where $B_{T}$ is the open ball with radius $\sqrt{2 S(T)}$.

Removing the compact support assumption {#app}
=======================================

The proof of the main result of the paper is built on the compact support assumption stated in Section \[sec:NAMR\]. Our goal is to generalize the result removing this condition. Let us note that it is sufficient to extend the Theorem \[Thm1\], since it is the only one for which we used the compact support assumption. Moreover, let us notice that Proposition \[prop1\] holds true in this framework. We define the functional $w$ as in and we state the following result.

Let $\Phi$ be a test function in $C_{\text{comp}}^2({\mathbb{R}}^d)$. With the assumptions from to , and with $\gamma > \max(1,2-\frac 4 d)$, for all $T>0$ there exists a constant $C(T)$ depending on the previous bounds and independent of $\gamma$ such that $$\int_0^T\int_{{\mathbb{R}}^d}|w|_-^3\Phi\leq C(T),\quad \int_0^T\int_{{\mathbb{R}}^d}|\nabla p|^2 |w|_-\Phi \leq C(T), \quad  \int_0^T \int_{{\mathbb{R}}^d} |\Delta {p_\gamma}| \Phi \leq C(T).$$

Computing the time derivative of the negative part of $w$, we have $$\begin{aligned}
            -\partial_t \left(\frac{|w|_-^2}{2}\right) \leq    
            &-\frac{4}{d}|w|_-^3-\frac{2}{d} G|w|^2_--\frac{2}{d}G^2|w|_-+\nabla |w|_-^2 \cdot \nabla p+\partial_p G |\nabla p|^2|w|_-\\ &+2 \partial_c G\nabla p \cdot \nabla c |w|_-
            +\gamma \Delta(p |w|_-)|w|_- - \partial_c G\, \partial_t c |w|_-.
            \end{aligned}$$ as in the proof of Theorem \[Thm1\]. We multiply the inequality by $\Phi$ and integrate in space $$\begin{aligned}
\label{phi} 
            \frac{d}{{dt}}\int_{\Omega_T}\frac{|w|_-^2}{2}\Phi\leq &-\frac{2}{d}\int_{\Omega_T}|w|_-^3\Phi -\frac{2}{d}\int_{\Omega_T}   G^2 |w|_-\Phi-\beta\int_{\Omega_T} |\nabla p|^2|w|_- \Phi
            \\\nonumber & -\frac{4}{d} \int_{\Omega_T}   G |w|_-^2 \Phi +\underbrace{\int_{\Omega_T}\left[\nabla p\cdot \nabla \left(|w|_-^2\right) \Phi + \gamma \Delta(p |w|_-)|w|_-\Phi \right] }_{{\mathcal{A}}} \\
            \nonumber&\underbrace{-\int_{\Omega_T} \partial_c G\,\partial_t c |w|_- \Phi}_{{\mathcal{B}}} +\underbrace{2\int_{\Omega_T} \partial_c G\nabla p \cdot \nabla c |w|_-\Phi}_{{\mathcal{C}}}.
            \end{aligned}$$ Now we proceed computing each term. $$\begin{aligned}
            {\mathcal{A}}=&\int_{{\mathbb{R}}^d}\nabla p \cdot\nabla \left(|w|_-^2\right) \Phi -\gamma \int_{{\mathbb{R}}^d} \nabla (p |w|_-)\cdot\nabla |w|_- \Phi -\gamma \int_{{\mathbb{R}}^d}|w|_- \nabla (p |w|_-)\cdot\nabla \Phi \\
            = &-\int_{{\mathbb{R}}^d}\Delta p |w|_-^2 \Phi - \int_{{\mathbb{R}}^d}|w|_-^2 \nabla p \nabla\Phi -\gamma \int_{{\mathbb{R}}^d}|w|_-\nabla p\nabla|w|_-\Phi \\
            &-\gamma \int_{{\mathbb{R}}^d}p |\nabla |w|_-|^2\Phi+\gamma \int_{{\mathbb{R}}^d}p|w|_-^2\Delta\Phi+\gamma \int_{{\mathbb{R}}^d} p \nabla \left(\frac{|w|_-^2}{2}\right)\cdot\nabla \Phi\\
            = &-\int_{{\mathbb{R}}^d}\Delta p |w|_-^2 \Phi - \int_{{\mathbb{R}}^d}|w|_-^2 \nabla p \cdot\nabla\Phi +\frac\gamma 2 \int_{{\mathbb{R}}^d}\Delta p|w|_-^2\Phi +\frac\gamma 2 \int_{{\mathbb{R}}^d}|w|_-^2\nabla p\cdot\nabla\Phi \\
            &-\gamma \int_{{\mathbb{R}}^d}p |\nabla |w|_-|^2\Phi +\frac{\gamma}{2} \int_{{\mathbb{R}}^d} p |w|_-^2 \Delta \Phi-\frac \gamma 2 \int_{{\mathbb{R}}^d}|w|_-^2 \nabla p \cdot\nabla \Phi\\
            =&\left(1-\frac\gamma 2\right)\int_{{\mathbb{R}}^d} |w|_-^3\Phi  +\left(1- \frac \gamma 2\right)\int_{{\mathbb{R}}^d} G |w|_-^2\Phi -\gamma \int_{{\mathbb{R}}^d} p |\nabla|w|_-|^2\Phi+{\mathcal{A}}_1,
            \end{aligned}$$ with $${\mathcal{A}}_1=\frac \gamma 2 \int_{{\mathbb{R}}^d}p|w|_-^2\Delta\Phi -\int_{{\mathbb{R}}^d}|w|_-^2\nabla p\cdot\nabla\Phi.$$ By the Cauchy-Schwarz inequality we have $$\begin{aligned}
                        {\mathcal{B}}&\leq\,  \int_{{\mathbb{R}}^d}|w|_-^2\Phi+C\int_{{\mathbb{R}}^d}|\partial_t c|^2\Phi\leq \int_{{\mathbb{R}}^d}|w|_-^2\Phi+C.
            \end{aligned}$$ Using Young’s inequality and , we find $$\begin{aligned}
            {\mathcal{C}}&\leq \frac{\beta}{2}\int_{{\mathbb{R}}^d} | \nabla p|^2|w|_- \Phi+ C\int_{{\mathbb{R}}^d} |\nabla c|^2|w|_-\Phi\\
            &\leq  \frac{\beta}{2}\int_{{\mathbb{R}}^d} |\nabla p|^2|w|_- \Phi+ C\int_{{\mathbb{R}}^d} |\nabla c|^4\Phi+ C\int_{{\mathbb{R}}^d}|w|_-^2\Phi \\
            &\leq  \frac{\beta}{2}\int_{{\mathbb{R}}^d} |\nabla p|^2|w|_- \Phi+ C\int_{{\mathbb{R}}^d}|w|_-^2\Phi +C.
            \end{aligned}$$ It remains to treat the term containing the derivatives of $\Phi$ $${\mathcal{A}}_1= -\int_{{\mathbb{R}}^d}|w|_-^2\nabla p\cdot \nabla\Phi+\frac \gamma 2 \int_{{\mathbb{R}}^d}p|w|_-^2\Delta\Phi.$$ We choose a positive function $\Phi$ with exponential decay, such that $|\nabla\Phi|\leq C\Phi$ and $|\Delta\Phi|\leq C\Phi$. Now, we integrate by parts and use Young’s inequality $$\begin{aligned}
            {\mathcal{A}}_1=&\,2\int_{{\mathbb{R}}^d}p|w|_-\nabla|w|_- \cdot\nabla\Phi +\left(1+\frac{\gamma}{2}\right)\int_{{\mathbb{R}}^d}p|w|_-^2  \Delta \Phi \\
            \leq & \,\frac{1}{2} \int_{{\mathbb{R}}^d}p|\nabla|w|_-|^2\Phi +  C (\gamma+1) \int_{{\mathbb{R}}^d}|w|_-^2\Phi.
            \end{aligned}$$ Finally, the inequality can be written as follows $$\frac{d}{{dt}}\int_{{\mathbb{R}}^d}|w|_-^2\Phi+ \left(\frac 2 d+\frac \gamma 2 -1\right)\int_{{\mathbb{R}}^d}|w|_-^3\Phi+\frac \beta 2\int_{{\mathbb{R}}^d}|\nabla p|^2 |w|_-\Phi\leq C(\gamma+1)\int_{{\mathbb{R}}^d} |w|_-^2\Phi  +C,$$ then, for $ \gamma>2-\frac 4 d$, integrating in time we have $$\int_0^T\int_{{\mathbb{R}}^d}|w|_-^3\Phi\leq \left(\int_0^T\int_{{\mathbb{R}}^d}|w|_-^3\Phi\right)^\frac{2}{3}+C(T),$$ and then we have proved $$\int_0^T\int_{{\mathbb{R}}^d}|w|_-^3\Phi\leq  C(T) .$$ By consequence $$\int_0^T\int_{{\mathbb{R}}^d}|w|_-^2\Phi\leq  C(T),\quad\int_0^T \int_{{\mathbb{R}}^d}|w|_-\Phi\leq  C(T),$$ and $$\int_0^T\int_{{\mathbb{R}}^d}|\nabla p|^2 |w|_-\Phi \leq C(T).$$ Since $\Phi$ is a smooth function with compact support $$\int_0^T \int_{{\mathbb{R}}^d} (\Delta p + G)\Phi\leq C,$$ and then also $$\int_{{\mathbb{R}}^d}\Phi|\Delta p + G|_+=\int_{{\mathbb{R}}^d}\Phi(\Delta p + G)+\int_{{\mathbb{R}}^d}\Phi|\Delta p + G|_-\leq C(T).$$ Therefore we recover the local $L^1$ estimate for the Laplacian of the pressure $$\int_0^T \int_{{\mathbb{R}}^d} |\Delta p|\Phi \leq C.$$

Optimality of the bound $\nabla p \in L^4$
==========================================

In Theorem \[L4\], we have established the uniform bound $\nabla p \in L^4$, see . Here, we aim to show that the exponent $4$ cannot be increased. We use the so-called *focusing solution* from [@AG] that we adapt to the limit $\gamma \to \infty$, i.e., the Hele-Shaw problem. We recall that, for the porous medium equation, the focusing solution consists in a spherical hole filling which generates a stronger singularity than the Barenblatt solution, see [@V].

Consider $\alpha>0$ such that $\nabla p \in L^\alpha(Q_T)$, where $p$ is a solution of the Hele-Shaw problem with Dirichlet boundary conditions in a spherical shell $\{R(t) < | x| < R_1\}$, for a given $R_1>0$ and $R(0)$ small enough. Then, to simplify the problem, we fix the external radius $R_1$ and let $p$ satisfy $$\label{p'} 
\begin{cases}
-\Delta p = 1,\quad  &\text{ for  }  \quad R(t) < | x| < R_1,
\\[2pt]
p(x)=0, \quad  &\text{ for } \quad |x|= R(t) \ \text{ or } |x|= R_1 ,
\\[2pt]
R'(t)=-\nabla p \cdot \vec{n},  &\text{ for } \quad |x|= R(t).
\end{cases}$$ Here, $\vec{n}$ denotes the inner normal to the ball $B_{R(t)}(0)$. As in [@AG], $R(t) $ diminishes and vanishes in finite time, generating a singularity $|\nabla p| \to \infty$. The power $4$ turns out to be the highest possible integrability in time at this singular time. We treat the case of dimension 2. In higher dimension the radial solutions are more regular and the worse singularity would be obatined for a cylinder with a 2 dimensional basis.

#### Case $d=2$. {#case-d2. .unnumbered}

With spherical symmetry, we set $p:=p(r)$, $r=|x|$ and equation reads $$\frac 1 r (rp')' =1.$$ Integrating once, we get, for some $a(t)$ $$p'=-\frac{r}{2}+\frac{a(t)}{r},$$ and the second integration yields $$p=-\frac{r^2}{4}+a(t) \ln r +b(t).$$ Imposing $p(R_1)=p(R(t))=0,$ we find $$\begin{aligned}
b(t)&=\frac{R_1^2}{4}-a(t) \ln R_1,\\
\frac{R(t)^2}{4}- a(t) \ln R(t) &= \frac{R_1^2}{4}-a(t) \ln R_1.\end{aligned}$$ Hence for $R(t)\approx 0$, we have $$\label{a}
a(t) \approx - \frac{R_1^2}{4 \ln R(t)}, \qquad \qquad R'(t)\approx\frac{1}{ R(t) \; \ln R(t) } .$$ Therefore there is $T>0$ when $R(T^-)=0$ and as $t \approx T$, we compute $$\int_0^T \int_{B_{R_1}(0)} |\nabla p(x)|^\alpha dx dt = \int_0^T \int_{R(t)}^{R_1} |p'(r)|^\alpha r   dr dt\approx\int_0^T \int_{R(t)}^{R_1} \frac{|a(t)|^\alpha}{r^{\alpha -1}}   dr dt.$$ The singularity at $T$ is thus driven by $$\int_0^T\frac{|a(t)|^\alpha}{R(t)^{\alpha -2}} dt\approx  \int_0^T\frac{1}{| \ln R(t) |^{\alpha}\, R(t)^{\alpha-2}} dt \approx 
\int_{0} ^{R(0)} \frac{1}{|\ln R|^{\alpha-1} R^{\alpha-3}} dR$$ by the change of variable $R=R(t)$ and using equation . We recall that we have chosen $R(0)$ small enough.

This integral is finite for $1\leq \alpha \leq 4$ and infinite for $\alpha >4$.

[^1]: Sorbonne Universit[é]{}, Inria, CNRS, Université de Paris, Laboratoire Jacques-Louis Lions UMR7598, F-75005 Paris. Email : noemi.david@ljll.math.upmc.fr

[^2]: Sorbonne Universit[é]{}, CNRS, Université de Paris, Inria, Laboratoire Jacques-Louis Lions UMR7598, F-75005 Paris. Email : Benoit.Perthame@sorbonne-universite.fr
