# Free boundary regularity for tumor growth with nutrients and diffusion

Carson Collins, Matt Jacobs and Inwon Kim

###### Abstract.

In this paper, we study a tumor growth model where the growth is driven by
nutrient availability and the tumor expands according to Darcy’s law with a
mechanical pressure resulting from the incompressibility of the cells. Our
focus is on the free boundary regularity of the tumor patch that holds beyond
topological changes. A crucial element in our analysis is establishing the
regularity of the hitting time $T(x)$, namely the first time the tumor patch
reaches a given point. We achieve this by introducing a novel Hamilton-Jacobi-
Bellman (HJB) interpretation of the pressure, which is of independent
interest. The HJB structure is obtained by viewing the model as a limit of the
Porous Media Equation (PME) and building upon a new variant of the AB
estimate. Using the HJB structure, we establish a new Hopf-Lax type formula
for the pressure variable. Combined with barrier arguments, the formula allows
us to show that $T$ is $C^{\alpha}$ with $\alpha=\alpha(d)$, which translates
into a mild nondegeneracy of the tumor patch evolution. Building on this and
obstacle problem theory, we show that the tumor patch boundary is regular in
${\mathbb{R}}^{d}\times(0,\infty)$ except on a set of Hausdorff dimension at
most $d-\alpha$. On the set of regular points, we further show that the tumor
patch is locally $C^{1,\alpha}$ in space-time. This conclusively establishes
that instabilities in the boundary evolution do not amplify arbitrarily high
frequencies.

##  1\. Introduction

In this paper, we consider the following tumor growth model:

(1.1) $\partial_{t}\rho-\nabla\cdot(\rho\nabla p)=n\rho,\quad
p(1-\rho)=0,\quad\rho\leq 1,$

where $\rho$ denotes the density of tumor cells, $p$ denotes the pressure, and
$n$ is a nutrient variable that evolves according to the diffusion equation

(1.2) $\partial_{t}n-\Delta n=-n\rho.$

The form of the pressure-density relation reflects the incompressibility of
the tumor cells, namely the pressure variable $p$ acts as the Lagrange
multiplier for the constraint $\rho\leq 1$. In short, the system (1.1-1.2)
describes a cell growth system where the growth rate is mediated by nutrient
availability and the tumor region expands according to Darcy’s law with a
mechanical pressure driven by the incompressibility of the cells. Models of
the form (1.1-1.2) have been extensively studied by both the mathematical and
biological communities with various different assumptions on the growth term
and density pressure coupling [BKMP03, PT08, RBE+10, MRCS14, PQV14], to name
just a few. Nonetheless, many mathematical questions remain outstanding, in
particular, those regarding the long-time behavior of the tumor boundary
region.

Our focus on the specific source term $n\rho$ is due to the fact that the
model (1.1-1.2) generates particularly interesting behavior of the tumor patch
despite the apparent simplicity of the coupling between the tumor and
nutrient. It is well-known in the biology literature (through numerical and
physical experiments) that the tumor patch generated by this model exhibits a
fingering instability (c.f. the discussion in [Kit97], [MVS02], [MRCS14],
[GDSK22], [JKT23]). In particular, it has been unclear whether this fingering
phenomenon occurs at some discrete scale or whether it leads to an immediate
or eventual loss of regularity in the tumor boundary. Investigating this
behavior will be the main goal of this paper.

Although the tumor system nearly corresponds to that of the classical Hele-
Shaw flow, a mathematically rigorous study of the boundary behavior has
remained elusive, due to the difficulties presented by the source term $\rho
n$. In the classical setting, which we will call the injection problem, the
Hele-Shaw flow is given with no source (namely $n=0$) and with a fixed
boundary from which the flow is injected at a given rate. For the injection
problem, the global structure of the boundary $\partial\\{\rho=1\\}$ is well
understood by now, mainly through comparison principle type arguments [CJK07,
CJK09, DGN21] or via connections to the obstacle problem [BCMP73, Mon03, FS19,
FROS20].

For our problem, the comparison approach is immediately ruled out, as the full
system (1.1-1.2) does not have comparison (though note that the individual
equations when considered separately do have comparison principles). As such,
we shall proceed via the obstacle problem analysis. However, there is a highly
nontrivial roadblock that must be overcome. Indeed, the source term $\rho n$
necessarily depends on the space-time geometry of the free boundary, while for
the injection case, the source is concentrated at a fixed boundary that is
safely away from the free boundary. This makes the analysis of the tumor
system considerably more difficult, as the influence of the source term cannot
be ignored when blowing up the problem at free boundary points (the
fundamental technique for the obstacle problem approach). In particular, to
use the obstacle problem toolbox, one must first establish the regularity of
the hitting time  $T(x)$, which records the first time that the tumor patch
reaches the point $x$ (ignoring the regularity issues, one can formulate
$T(x)$ as $\inf\\{t>0:\rho(t,x)=1\\}$, see equation (1.9) for a more careful
definition). This is essentially equivalent to establishing a quantitative
non-degeneracy property for the tumor expansion speed, a highly nontrivial
task.

To establish the regularity of $T(x)$ we first derive a novel Hopf-Lax type
estimate for the pressure (c.f. Theorem 1.1). To the best of our knowledge,
such Hopf-Lax type formulas have not previously appeared in the Hele-Shaw
literature, perhaps in part due to the difficulty of controlling the time
derivative of $p$. We get around this by viewing equation (1.1) as the
incompressible limit of the Porous Media Equation (PME). Given some parameter
$\gamma\in(1,\infty)$, the PME analogue of (1.1) is the equation

(1.3) $\partial_{t}\rho_{\gamma}-\nabla\cdot(\rho_{\gamma}\nabla
p_{\gamma})=\rho_{\gamma}n_{\gamma},\quad p_{\gamma}=\rho_{\gamma}^{\gamma},$

where $n_{\gamma}$ will solve (1.2) with $\rho$ replaced by $\rho_{\gamma}$,
and (1.1) can be recovered by sending $\gamma\to\infty$ (see for instance
[PQV14, DP21, Jac21]). Since the pressure-density coupling
$p_{\gamma}=\rho_{\gamma}^{\gamma}$ is explicit for PME, one can rewrite (1.3)
solely in terms of the pressure, namely,

(1.4) $\partial_{t}p_{\gamma}-|\nabla p_{\gamma}|^{2}-\gamma p_{\gamma}(\Delta
p_{\gamma}+n_{\gamma})=0.$

Interestingly, we ignore the parabolic structure of this equation and instead
focus on the Hamilton-Jacobi-Bellman (HJB) structure of the first two terms.
We then build upon the recent improved versions of the Aronson-Benílan
estimate introduced in [Jac23] to show that the positive part of
$u_{\gamma}:=-\gamma(\Delta p_{\gamma}+n_{\gamma})$ is uniformly bounded with
respect to $\gamma$ in a BMO type space, implying that our limiting $p$ must
be a supersolution to the HJB equation

(1.5) $\partial_{t}p-|\nabla p|^{2}+pu_{+}\geq 0$

where $u:=\lim_{\gamma\to\infty}u_{\gamma}.$ From here we finally obtain the
Hopf-Lax formula by adapting the techniques of [CG15] for HJB equations with
unbounded coefficients. It is highly intriguing to speculate whether it is
possible to obtain (1.5) or Hopf-Lax estimates directly from (1.1), however,
we will not consider this line of inquiry further in this work.

Once we have the Hopf-Lax formula, we combine this with a powerful barrier-
type argument to prove that for any point $x\notin\textup{spt}(\rho_{0})$ and
any sufficiently small radius $r>0$ there exists an explicit time
$t_{r}(x)<T(x)$ such that the tumor patch does not occupy any point in
$B_{r}(x)$. From here, it will follow that the hitting time is Hölder
continuous with an exponent that depends on the dimension only. With the
Hölder continuity of $T$ in hand, we can turn to the obstacle problem
formulation to address the regularity of the free boundary. Here, the novelty
in our analysis lies in establishing the global space-time regularity of the
free boundary, with data that is far less regular than the typical injection
problems that have previously been considered.

Ultimately, through the obstacle problem analysis, we are able to show that
the free boundary is regular except at topological singularities, which are
unavoidable for general initial data. This conclusively demonstrates that the
observed instabilities for the system (1.1-1.2) do not amplify arbitrarily
high frequencies and must occur at some fixed scale. In particular, we show
that the tumor patch boundary is regular in ${\mathbb{R}}^{d}\times(0,\infty)$
except on a relatively closed set of Hausdorff dimension at most $d-\alpha$
for some $\alpha\in(0,1)$ depending only on the dimension. On the set of
regular points, we further show that the tumor patch is $C^{1,\alpha}$ in
space, locally uniformly in time. It then follows that the associated pressure
gradient at regular boundary points is well-defined and uniformly positive in
space-time. Moreover, the direction of the pressure gradient on the set of
regular points is continuous in space-time.

In the remainder of the introduction, we give a more complete explanation of
the obstacle formulation of our problem and the connection to the hitting
time. We then summarize our main results and give a roadmap for the rest of
the paper.

###  1.1. The obstacle problem and the hitting time

To better understand the aforementioned difficulties and the importance of the
hitting time, let us describe some properties of the tumor patch and formally
introduce the obstacle problem associated to (1.1-1.2). Since our main
interest is the regularity properties of the tumor patch, throughout the
paper, we will assume that

(1.6) $\rho(x,0)\hbox{ is a characteristic function and }n(x,0)\hbox{ is
uniformly positive.}$

Under these assumptions, $\rho$ will remain a characteristic function for all
times and $t\mapsto\rho(x,t)$ will be nondecreasing for a.e.
$x\in\mathbb{R}^{d}$.

Transitioning to the obstacle problem formulation, if we integrate the
pressure variable in time,

(1.7) $w(x,t):=\int_{0}^{t}p(x,s)ds,$

the new variable $w$, the so-called Baoicchi transform, will satisfy an
obstacle problem [BCMP73]. Since the density is nondecreasing in time, the
relation $(1-\rho)p=0$ implies that $(1-\rho)w=0$. Using the patch property
for the density, this coupling can be upgraded to the even stronger relation
that the sets $\\{w>0\\}$ and $\\{\rho=1\\}$ coincide spacetime almost
everywhere (c.f. Lemma 2.7). This key relation can then be combined with the
time integral of (1.1) to see that $w$ solves the elliptic obstacle problem

(1.8) $\Delta w=(1-\rho_{0}-\eta)\chi_{\\{w>0\\}},$

where $\eta(x,t):=\int_{0}^{t}\rho(x,s)n(x,s)\,ds$ (c.f. Lemma 2.9).

The main challenge in analyzing (1.8) is the presence of the term $\eta$,
which is absent in the obstacle formulation of the classical injection case
(due to local regularity results, $\rho_{0}$ does not affect the free boundary
regularity at positive times away from the support of $\rho_{0}$). Since
$\rho$ is a characteristic function, it is not clear whether $\eta$ has any
nice regularity. This is crucial, as obstacle problem regularity theory breaks
down without Dini continuity of the coefficients (see [Bla00]). Hence, one
must hope that the time integral induces some smoothing effect. At the very
least, this can only happen if the tumor boundary is strictly expanding.
Indeed, if any part of the free boundary stagnates in time, then $\eta$ will
become discontinuous across that portion of the boundary. Note that such
stagnation would correspond to a jump in the values of the hitting time
function $T$ introduced earlier. Hence, the smoothness of $\eta$ and $T$ are
highly intertwined. In fact, it will turn out that we can express $\eta$
solely in terms of the hitting time $T$ and $n$.

To see the connection between $\eta$ and $T$, we need to first give a proper
definition of the hitting time. Recall that the hitting time $T(x)$ records
the first time that the tumor patch arrives at a point $x$. We will formally
define it using $w$, the most regular variable at our disposal. Given a point
$x\in\mathbb{R}^{d}$ we set

(1.9) $T(x):=\inf\\{t>0:w(x,t)>0\\}.$

Since the positivity set of $w$ coincides almost everywhere with the tumor
patch, we have $\rho(x,t)=\operatorname{sgn}_{+}(t-T(x))$ almost everywhere.
Hence, $\eta$ can be rewritten in terms of $T$ and $n$ as

(1.10) $\eta(x,t)=\operatorname{sgn}_{+}(t-T(x))\int_{T(x)}^{t}n(x,s)\,ds.$

From the above formula, we now see that the spatial regularity of $\eta$ is
more or less equivalent to the regularity of $T$ and $n$.

Note that generically $T$ is at best Lipschitz continuous, as it is easy to
cook up a scenario where two different parts of the tumor patch collide with
different velocities. In addition, topological changes of the tumor boundary
can cause the pressure to suddenly jump with highly nonlocal effects. For
instance, the merger of two portions of the boundary can cause far away parts
of the boundary to instantaneously start moving faster. Since $n$ is much
better than Lipschitz continuous, it is $T$ that will determine the regularity
of $\eta$. While we are inclined to believe that the Lipschitz continuity of
$T$ is true, our methods are only able to show that $T$ is Hölder continuous
with a dimensionally dependent exponent. Nonetheless, the Hölder continuity is
sufficient for us to deduce free boundary regularity using the obstacle
problem approach. However, let us note that we are forced to work in a much
lower regularity regime than what is typically considered in the obstacle
problem literature, requiring us to develop new arguments.

###  1.2. Main results

We are now ready to present the main results of our paper. All of our results
will use the following mild assumptions on the initial data.

  1. (A1)

$\rho(\cdot,0)\in L^{1}(\mathbb{R}^{d})\cap\textup{BV}(\mathbb{R}^{d})$ and
$\rho(x,0)\in\\{0,1\\}$ for almost every $x\in\mathbb{R}^{d}$.

  2. (A2)

$n(\cdot,0)\in W^{1,\infty}(\mathbb{R}^{d})$ and there exists $c>0$ such that
$n(x,0)\geq c$ for all $x\in\mathbb{R}^{d}$.

The main results of the first half of the paper are the HJB structure and
Hopf-Lax formula for the pressure, along with the Hölder continuity of the
hitting time.

######  Theorem 1.1.

The following holds for the unique weak solution $p$ to the system
(1.1)-(1.2).

  * (a)

[Cor 3.4] $p$ solves, in the sense of weak solutions,

$\partial_{t}p-|\nabla p|^{2}+pu_{+}\geq 0,$

where for any $\tau>0$ there exists $b=b(\tau,d)>0$ such that
$bu_{+}e^{bu_{+}}\in L^{1}([0,\tau];{\mathbb{R}}^{d})$.

  * (b)

[Prop. 3.7] Given points $(x_{1},t_{1})$, $(x_{0},t_{0})$ with $t_{0}<t_{1}$
and any decreasing function $\lambda\in L^{1}([0,t_{1}-t_{0}])$, there exist
constants $C=C(t_{1},d)$ and $b=b(t_{1},d)$ such that

$p(x_{0},t_{0})\leq
e^{\Lambda(t_{1}-t_{0})}\Big{(}p(x_{1},t_{1})+\frac{|x_{1}-x_{0}|^{2}}{4\int_{0}^{t_{1}-t_{0}}e^{\Lambda(s)}\,ds}+C(t_{1}-t_{0})^{7/10}e^{-\lambda(t_{1}-t_{0})}\Big{)}$

where
$\Lambda(t):=\frac{5}{4b}\int_{0}^{t}\lambda(s)\,ds+\frac{t}{b}\log(1+\frac{C}{t})$.

  * (c)

[Theorem 4.2] $T$ is locally Hölder continuous on the set
$\\{x\in\mathbb{R}^{d}:0<T(x)<\infty\\}$ with an exponent that depends only on
the dimension.

Let us note that Theorem 1.1 parts (a) and (b) represent a significant
improvement to our understanding of the Hele-Shaw equation. In particular, any
control on the time derivative of the pressure has been previously missing in
the literature. Furthermore, the delicate control that we obtain from the
Hopf-Lax formula in part (b) is completely new and unexpected.

As we mentioned earlier, we establish the HJB structure by first going through
the PME (3.2). For the classic PME without a source term, bounds on the
negative part of $\Delta p_{\gamma}$ are known through the celebrated Aronson-
Benilan estimate [AB79]. In the presence of a source term, AB-type bounds on
quantities taking a similar form to $u_{\gamma}=-\gamma(\Delta p_{\gamma}+n)$
have been studied in the literature [PQV14, GPŚG19, DP21, BPS22, Jac23],
however except for [Jac23], these bounds do not scale well with respect to
$\gamma$. We adapt the arguments from [Jac23] to show that $[u_{\gamma}]_{+}$
can be bounded uniformly with respect to $\gamma$ in a BMO-type space. Once we
have the uniform control on $u_{\gamma,+}$ we can pass to the limit in (3.4)
to obtain the result (a). A direct derivation of (a) from the Hele-Shaw flow
or the meaning of the singular limit $u=\lim_{\gamma\to\infty}u_{\gamma}$ in
terms of the Hele-Shaw flow remains open.

To obtain (b), we cannot take the usual approach to proving Hopf-Lax type
formulas (i.e. differentiating $p$ along paths) due to the potential
unboundedness of $u_{+}$. To overcome this, we adapt the approach developed in
[CG15], which handles unbounded coefficients by averaging over paths indexed
by the unit ball. Our calculation is somewhat different however, as we can
exploit the specific structure of $pu_{+}$ to decompose $pu_{+}\leq\lambda
p+p(u-\lambda)_{+}$ for some scalar $\lambda\geq 0$. By choosing $\lambda$
appropriately we can force $p(u-\lambda)_{+}$ to be small while using a
Gronwall argument to handle $\lambda p$. This allows us to obtain a much more
favorable error term in our Hopf-Lax formula compared to [CG15].

Although Theorem 1.1 (a) and (b) are stated for our particular system
(1.1-1.2), equivalent results can be proved for more general tumor growth
models where the growth term $\rho n$ is replaced by $\rho G$ for some general
growth rate $G$. In particular, our arguments only need $G\geq c(\tau)>0$
along with some control on $[\partial_{t}G]_{-}$.

As a consequence of the Hopf-Lax formula, we obtain the $C^{\alpha}$
regularity of $T$ where $\alpha=\alpha(d)$. We do this by combining the
formula with a novel barrier type argument. Given a point $(x,T(x))$ on the
free boundary and some time $t_{0}\in(0,T(x))$, we use the Hopf-Lax formula
and the values of $p$ at time $T(x)$ to construct an explicit supersolution
$\psi$ that dominates $p$ on $(t_{0},T(x))\times\mathbb{R}^{d}$. The key is
that the Hopf-Lax formula allows us to choose the values in such a way that
$\psi$ is zero in a neighborhood of $x$ up until the hitting time $T(x)$.
Since we are able to explicitly calculate and invert
$t(r)=\inf\\{t\in(t_{0},T(x)):\sup_{y\in B_{r}(x)}\psi(y,t)>0\\}$ we obtain an
upper bound on $\sup_{y\in B_{r}(x)}T(x)-T(y)$, which implies the Hölder
continuity of $T$.

Some remarks on the previous literature for hitting times are in order.
Quantitative regularity of the hitting time for PME has been obtained in
[CF80] for the classical PME and in [KZ21] for the PME with a source term and
drift. Nevertheless, both of these results obtain estimates that blow up as
$\gamma$ tends to infinity, due to the lack of a uniform AB estimate on
$u_{\gamma}$. As a result, their approaches are not suitable for our problem.
Let us also note that these papers used a rather different approach that did
not involve the Hopf-Lax approach that we use here. Estimates on the hitting
time for a simpler version of (1.1) where $n$ is replaced by a decreasing
function of $p$, were obtained in [MPQ17] for dimensions $d\leq 3$. Their
proof strongly relies on the specific structure of their growth term, which
allows them to relate the Hölder continuity of $T$ to that of the pressure
through a clever trick. Again, this approach is not applicable to our problem.
Although we also focus on a specific source term, our method is much more
general and can be applied to other instances of the Hele-Shaw or Porous Media
equation.

The remaining analysis in the paper is devoted to the study of the obstacle
problem (1.8), based on the $C^{\alpha}$ regularity of $T$. We build on the
low-regularity obstacle problem analysis of Blank [Bla00] to establish the
space-time regularity of the tumor patch. A crucial fact we use is that the
solution of the obstacle problem with $C^{\alpha}$ data has a unique blow-up
limit at each point, allowing us to decompose the boundary into a regular part
and a singular part (the regular points have blow up limits that look like
half-planes). A direct application of this dichotomy yields that the boundary
has locally finite $H^{d-1}$ measure for each time, as mentioned for instance
in [MPQ17]. However, this standard description lacks the geometric information
of the free boundary over time. Indeed, the main novelty of our obstacle
problem analysis is that we are able to stitch together information from each
time $t$ to obtain regularity of the full space time boundary
$\Gamma:=\\{(x,T(x):x\in\mathbb{R}^{d}\\}\subset\mathbb{R}^{d}\times(0,\infty)$.
In particular, we show that $\Gamma$ is regular in space-time outside of a set
of at most Hausdorff dimension $d-\alpha$, and its outward normal is Hölder
continuous in space-time. While space-time analysis of the singular set has
been carried out before for the injection problem ([Mon03], [FS19], [FROS20]),
these results have utilized smoothness (at least $C^{4}$) of the fixed
boundary data in an essential way. A more general time-varying source term was
considered in [SS17], but only for a short range of time that ensures that no
topological singularity occurs during the evolution.

Our results are summarized in the following Theorem.

######  Theorem 1.2.

Let $\Gamma$ denote the space-time boundary set of the tumor region i.e.
$\Gamma=\\{(x,T(x)):x\in\mathbb{R}^{d}\\}$

  * (a)

[Prop. 5.4, Prop. 5.15] The set $\\{0<T(x)<\infty\\}\subset{\mathbb{R}}^{d}$
decomposes as $R\cup\Sigma$, where the set $R$ of regular points is open in
${\mathbb{R}}^{d}$ and the set $\Sigma$ of singular points is locally
contained in a $C^{1}$ manifold of dimension $d-1$.

  * (b)

[Prop. 5.12] At any $x\in R$, the free boundary near $(x,T(x))\in\Gamma$ can
be locally represented as a graph $\\{x_{n}=f(x^{\prime},t)\\}$ where $f$ is
$C^{1,1}$ in $x^{\prime}$ and Lipschitz in time.

  * (c)

[Prop. 5.9, Cor.5.10 ] $p(\cdot,T(x))$ has linear growth at $x\in R$, with
locally uniform growth rates. In particular $T$ is Lipschitz in $R$.

  * (d)

[Prop. 5.6] The map $\nu:R\to\mathcal{S}^{d}$, where $\nu(x)$ denotes the
spatial outward normal of $\Gamma$ at $(x,T(x))$, is Hölder continuous. In
particular $\nabla p(x,T(x))$ is well-defined for $x\in R$ and has continuous
direction.

Note that, while Theorem 1.2 (d) yields the continuity of the direction of
$\nabla p$ on $\Gamma$, we cannot expect the same for $|\nabla p|$: this can
be easily seen from examples where a topological change occurs far away from
the given free boundary point. In terms of the quadratic blow up limit of $w$,
Theorem 1.2 (a) and (d) yield its continuity at free boundary points, along
$\Sigma$ and along $R$. As a consequence, $D^{2}w(x,T(x))$ exists on $\Sigma$
and exists in a one-sided sense on $R$, and is continuous on each set (though
not necessarily on their union).

Our last theorem discusses the Hausdorff measure of the free boundary in
space-time coordinate. Let us introduce the notation

(1.11) $\Omega_{t}:=\\{w(\cdot,t)>0\\},\quad\Gamma_{t}:=\partial\Omega_{t}.$

######  Theorem 1.3 (Corollary 5.16).

  * (a)

The free boundary $\partial\\{w>0\\}$ has Hausdorff dimension d in
$(x,t)$-coordinates.

  * (b)

$Graph(R)=\\{(x,t):x\in\Gamma_{t},x\in R\\}$ is relatively open with locally
finite $H^{d}$ measure.

  * (c)

$Graph(\Sigma)=\\{(x,t):x\in\Gamma_{t},x\in\Sigma\\}$ has locally finite
$H^{d-\alpha}$ measure.

Let us mention that we expect $T$ to be Lipschitz for all points, not just in
$R$. For instance, in the classical setting with constant Dirichlet fixed
boundary data, the Lipschitz continuity of $T$ was shown by [Mon03] by a
simple comparison principle. The remaining challenge in our setting lies in
the analysis of the singular points. This is an intriguing question as the
blow-up profile of the tumor patch at these points suggests that the evolution
at these points should be non-degenerate in general. In fact, one might even
expect the gradient of $T$ to vanish at these points. Nonetheless, accurately
capturing the hitting time behavior near singular points appears to be out of
reach for the moment. It would also be interesting to improve upon our
estimate of the singular set, using a generic notion of initial data. While it
seems plausible, new ideas seem to be necessary to obtain such a result.

The rest of the paper is organized as follows. In Section 2, we review the
basic properties of the system (1.1-1.2) and the connection to the obstacle
problem. In Section 3, we establish the HJB structure of the pressure along
with the Hopf-Lax formula. In Section 4, we construct a barrier supersolution
using the Hopf-Lax formula which allows us to establish the Hölder regularity
of the hitting time map $T(x)$. Section 5 builds on the regularity of $T$ and
the existing obstacle theory to investigate the global regularity of the free
boundary $\Gamma:=\partial\\{w(x,t)>0\\}$.

##  2\. Basic properties of the system

Here we recall the notion and show basic properties of solutions for
(1.1)-(1.2). We first introduce the weak notion of our solutions, parallel to
those introduced in [JKT23] for a similar model.

######  Definition 2.1.

A triple $(\rho,p,n)$ is a weak solution to (1.1)-(1.2) for initial data
$\rho_{0}\in L^{1}({\mathbb{R}}^{d})\cap BV({\mathbb{R}}^{d})$ and $n_{0}\in
L^{\infty}({\mathbb{R}}^{d})\cap BV({\mathbb{R}}^{d})$ if for any $\tau>0$,

  1. (i)

$p(1-\rho)=0$ in $\mathcal{D}^{\prime}({\mathbb{R}}^{d}\times[0,\tau])$

  2. (ii)

For any $\psi\in H^{1}({\mathbb{R}}^{d}\times[0,\tau])$ vanishing at time
$\tau$,

(2.1) $\int_{0}^{\tau}\int_{{\mathbb{R}}^{d}}\nabla\psi\cdot\nabla
p-\rho\partial_{t}\psi\,dx\,dt=\int_{{\mathbb{R}}^{d}}\psi(x,0)\rho_{0}(x)\,dx+\int_{0}^{\tau}\int_{{\mathbb{R}}^{d}}\psi
n\rho\,dx\,dt.$

  3. (iii)

$\partial_{t}n-\Delta n=-\rho n\hbox{ in
}\mathcal{D}^{\prime}({\mathbb{R}}^{d}\times[0,\tau]),\quad n(x,0)=n_{0}(x)$

  4. (iv)

We have $\rho\in C([0,\tau];L^{1}({\mathbb{R}}^{d}))\cap
L^{\infty}([0,\tau];BV({\mathbb{R}}^{d}))$, $p\in
L^{2}([0,\tau];H^{1}({\mathbb{R}}^{d}))$, and $n\in L^{\infty}(Q_{\tau})\cap
L^{\infty}([0,\tau];BV({\mathbb{R}}^{d}))$.

Here, $BV({\mathbb{R}}^{d})$ is the space of (not necessarily integrable)
functions with finite total variation.

We also record a few useful properties of a weak solution.

######  Lemma 2.2.

For any $\tau>0$, $\varepsilon\in(0,1)$, we have

  1. (i)

$\rho\in C^{0,1}_{t}L^{1}_{x}([0,\tau];{\mathbb{R}}^{d})$.

  2. (ii)

The support of $\rho$ in ${\mathbb{R}}^{d}\times[0,\tau]$ is compact, and
$0\leq\rho\leq 1$.

  3. (iii)

For a.e. $x\in{\mathbb{R}}^{d}$, $\rho(x,\cdot)$ is increasing in time.

  4. (iv)

$n\in
L^{\infty}([0,\tau]\times\mathbb{R}^{d}),\;\partial_{t}n\in\textup{BMO}([0,\tau]\times\mathbb{R}^{d}),\;D^{2}n\in\textup{BMO}([0,\tau]\times\mathbb{R}^{d}).$

  5. (v)

We have

(2.2) $p(\Delta p+n)=0\hbox{ in
}\mathcal{D}^{\prime}({\mathbb{R}}^{d}\times[0,\tau])$

Also, $p\in L^{\infty}({\mathbb{R}}^{d}\times[0,\tau])$.

###### Proof.

Statements (i) and (ii) are proved in [JKT23], Theorem 2.2 and Proposition
3.6.

(iii) follows from Lemma 3.11 of [JKT23], which provides comparison for the
equation

$\partial_{t}\rho^{i}-\nabla\cdot(\rho^{i}\nabla p^{i})=f^{i}$

Namely, if $\rho^{0}(x,0)\leq\rho^{1}(x,0)$, $f^{0}\leq f^{1}$, and
$p^{i}(1-\rho^{i})=0$, then $\rho^{0}(x,t)\leq\rho^{1}(x,t)$. For our system,
since $n\rho\geq 0$, comparison to the system with $p^{0}=f^{0}=0$ implies
that $\rho(x,t_{0})\leq\rho(x,t_{1})$ for any $t_{0}\leq t_{1}$ and a.e. $x$.
The measure zero set where this fails depends on $t_{0},t_{1}$, so we conclude
by applying this with a countable basis of intervals.

Item (iv) follows from parabolic estimates for the heat equation with
$L^{\infty}$ coefficients (see e.g. [OS22]).

The distributional equation $p(\Delta p+n)=0$ is also proved in [JKT23],
Theorem 2.2. The $L^{\infty}$ bound for $p$ follows from the compact support
of $\rho$, and the boundedness of $n$; one can take a sufficiently large
paraboloid supersolution to $\Delta p=-n$ to get an upper bound. ∎

From the definition of the weak solution, we can derive that $w$ satisfies an
elliptic equation at each time.

######  Lemma 2.3.

We have $w(1-\rho)=0$ a.e. in space-time. For each $t>0$, $w$ solves

(2.3) $\Delta w(x,t)=\rho(x,t)-\rho(x,0)-\int_{0}^{t}n(x,s)\rho(x,s)\,ds\hbox{
in }\mathcal{D}^{\prime}({\mathbb{R}}^{d})$

In particular, for any $\tau>0$ and $\varepsilon\in(0,1)$, we have $w\in
C^{0,1}({\mathbb{R}}^{d}\times[0,\tau])\times
L^{\infty}_{t}C^{1,1-\varepsilon}([0,\tau],{\mathbb{R}}^{d})$.

###### Proof.

For the first statement, we simply note that the monotonicity of $\rho$ in
time from Lemma 2.2 implies

$0\leq
w(x,t)(1-\rho(x,t))=\int_{0}^{t}p(x,s)(1-\rho(x,t))\,ds\leq\int_{0}^{t}p(x,s)(1-\rho(x,s))\,ds\equiv
0$

For the second statement, we first note that both $\rho$ and the function
$\int_{0}^{t}n\rho\,ds$ are continuous in time into any $L^{p}$ with
$p<\infty$; this follows from the weak solution definition, since we have
$\rho\in C_{t}L^{1}$ with values in $[0,1]$ and compact support for bounded
time intervals, while $n$ is spacetime continuous from Lemma 2.2. Then to
derive a distributional equation for $w$, we consider Definition 2.1(ii) with
$\psi$ of the form $\psi(x,t)=\varphi(x)\chi(t)$, where $\varphi\in
C^{\infty}_{c}({\mathbb{R}}^{d})$ and $\chi\in C^{\infty}_{c}([0,\tau))$ with
$\chi\equiv 1$ near 0. Using that $w=\int p\,dt$, we obtain

$\int_{{\mathbb{R}}^{d}}\chi\nabla\varphi\cdot\nabla
w(\cdot,\tau)\,dx-\int_{{\mathbb{R}}^{d}}\int_{0}^{\tau}\partial_{t}\chi\varphi\rho\,dx\,dt=\int_{{\mathbb{R}}^{d}}\varphi(x)\rho_{0}(x)\,dx+\int_{{\mathbb{R}}^{d}}\varphi\int_{0}^{\tau}\chi
n\rho\,dt\,dx$

If we take for $\chi$ a sequence of cutoffs valued in $[0,1]$ and converging
pointwise to the indicator of $[0,T)$, then we can apply the aforementioned
time continuity of $\rho$ and $\int_{0}^{t}n\rho\,ds$ to obtain the limiting
equation

$\int_{{\mathbb{R}}^{d}}\nabla\varphi\cdot\nabla
w(\cdot,\tau)+\varphi\rho(\cdot,\tau)\,dx=\int_{{\mathbb{R}}^{d}}\varphi(\rho_{0}+\eta(\cdot,\tau))\,dx$

from which we conclude (2.3).

Since $w(x,t)=\int_{0}^{t}p(x,s)\,ds$, the upper bound for $p$ implies that
$w$ is Lipschitz in time uniformly in space. This will improve to Lipschitz in
spacetime once we have $w\in L^{\infty}_{t}C^{1,1-}_{x}$.

Since $n$ is bounded on ${\mathbb{R}}^{d}\times[0,\tau]$, $\eta$ is bounded on
${\mathbb{R}}^{d}\times[0,\tau]$. Then since $\Delta
w=(1-\eta)\chi_{\\{w>0\\}}$ is uniformly bounded in $L^{\infty}$, and up to
time $\tau$, $w(\cdot,t)$ is compactly supported in
$\overline{\Omega_{\tau}}$, it follows that for any $p\in(1,\infty)$, $\Delta
w(\cdot,t)\in L^{p}({\mathbb{R}}^{d})$. Calderon-Zygmund estimates then give
$w(\cdot,t)\in W^{2,p}({\mathbb{R}}^{d})$, and thus $w(\cdot,t)\in
C^{1,1-\varepsilon}({\mathbb{R}}^{d})$ for any $\varepsilon>0$, uniformly in
$t\in[0,\tau]$. ∎

From now on we make the assumptions (A1) and (A2) on our initial data
$\rho_{0},n_{0}$.

######  Lemma 2.4.

Let $\bar{n}(t):=\inf_{x\in\mathbb{R}^{d}}n(t,x)$ For any $t>0$,
$\bar{n}(t)\geq e^{-t}\bar{n}(0).$

###### Proof.

Suppose that $\tilde{n}$ satisfies the equation
$\partial_{t}\tilde{n}-\alpha\Delta\tilde{n}=-\tilde{n}$ with constant initial
data $\tilde{n}(0,x)=\bar{n}(0).$ The comparison principle for the heat
equation implies that $\tilde{n}\leq n$ almost everywhere. If we define
$\tilde{N}=e^{t}\tilde{n}$, then $\tilde{N}$ satisfies
$\partial_{t}\tilde{N}-\alpha\Delta\tilde{N}=0$ with initial data
$\tilde{N}(0,x)=\bar{n}(0)$. Hence, $\tilde{N}(t,x)=\tilde{N}(0,x)=\bar{n}(0)$
and thus it follows that $\tilde{n}(t,x)=e^{-t}\bar{n}(0)$, which implies the
result. ∎

The following characterization of the pressure variable replaces the formal
description of $p$ solving the elliptic problem $-\Delta p=n$ in
$\\{\rho=1\\}$ with zero Dirichlet data, to avoid ambiguity rising from
potentially irregular boundary of $\\{\rho=1\\}$. The argument is similar to
ones that have previously appeared in the literature [PQV14, MRCS10, GKM22,
Jac21], here we include a proof since our setting is slightly different.

######  Lemma 2.5.

For almost every time $t$, the pressure is a solution to the variational
problem

$p(\cdot,t)=\operatorname*{argmin}_{\varphi(\cdot)(1-\rho(\cdot,t))=0,\;\varphi\geq
0}\,\int_{\mathbb{R}^{d}}\frac{1}{2}|\nabla\varphi(x)|^{2}-\varphi(x)n(x,t)\,dx.$

###### Proof.

Given $\epsilon>0$ define
$p_{\epsilon}(x,t):=\frac{1}{\epsilon}\int_{t-\epsilon}^{t}p(x,s)\,ds$ where
we set $p(x,s)=p(x,0)$ if $s<0$ and
$p^{\epsilon}(x,t):=\frac{1}{\epsilon}\int_{t}^{t+\epsilon}p(x,s)\,ds$. Fix a
time $t_{0}$ such that $p_{\epsilon}(t_{0},\cdot),p^{\epsilon}(t_{0},\cdot)$
converge to $p(t_{0},\cdot)$ in $H^{1}(\mathbb{R}^{d}).$ Choose some
nonnegative function $\varphi\in H^{1}(\mathbb{R}^{d})$ such that
$\varphi(x)(1-\rho(x,t_{0}))=0$ for almost every $x\in\mathbb{R}^{d}$ (note
that space integrals of $\rho(\cdot,t_{0})$ against functions in
$H^{1}(\mathbb{R}^{d})$ are well defined at any time $t_{0}$ since
$\partial_{t}\rho\in L^{2}([0,T];H^{-1}(\mathbb{R}^{d}))$, which itself is a
consequence of the continuity equation and $p\in
L^{2}([0,T];H^{1}(\mathbb{R}^{d}))$ ).

Integrating equation (1.1) from time $t_{0}-\epsilon$ to $t_{0}$, dividing by
$\epsilon$, and integrating against $\varphi$ we see that

$\int_{\mathbb{R}^{d}}\varphi(x)\frac{\rho(x,t_{0})-\rho(x,t_{0}-\epsilon)}{\epsilon}+\nabla\varphi(x)\cdot\nabla
p_{\epsilon}(x,t_{0})\,dx=\int_{\mathbb{R}^{d}}\varphi(x)\frac{1}{\epsilon}\int_{t_{0}-\epsilon}^{t_{0}}\rho(x,s)n(x,s)\,ds\,dx$

The condition $\varphi(x)(1-\rho(x,t_{0}))=0$ implies that
$\varphi(x)\frac{\rho(x,t_{0})-\rho(x,t_{0}-\epsilon)}{\epsilon}=\varphi(x)\frac{1-\rho(x,t_{0}-\epsilon)}{\epsilon}$.
Combined with the constraint $\rho\leq 1$, we can conclude that

$\int_{\mathbb{R}^{d}}\nabla\varphi(x)\cdot\nabla
p_{\epsilon}(x,t_{0})\,dx\leq\int_{\mathbb{R}^{d}}\varphi(x)\frac{1}{\epsilon}\int_{t_{0}-\epsilon}^{t_{0}}\rho(x,s)n(x,s)\,ds\,dx.$

Applying the same logic to the time integral over the interval
$[t_{0},t_{0}+\epsilon]$, we find that

$\int_{\mathbb{R}^{d}}\nabla\varphi(x)\cdot\nabla
p^{\epsilon}(x,t_{0})\,dx\geq\int_{\mathbb{R}^{d}}\varphi(x)\frac{1}{\epsilon}\int^{t_{0}+\epsilon}_{t_{0}}\rho(x,s)n(x,s)\,ds\,dx.$

Sending $\epsilon\to 0$ we can conclude that

(2.4) $\int_{\mathbb{R}^{d}}\nabla\varphi(x)\cdot\nabla
p(x,t_{0})=\int_{\mathbb{R}^{d}}\varphi(x)\rho(x,t_{0})n(t_{0},x)=\int_{\mathbb{R}^{d}}\varphi(x)n(x,t_{0})$

where the final equality follows from the fact that
$\varphi(x)(1-\rho(x,t_{0}))=0$ almost everywhere. The above equation is the
Euler-Lagrange equation for the variational problem, thus, combined with the
strong convexity of the variational problem, we see that $p$ solves the
variational problem at every time $t_{0}$ where
$p_{\epsilon}(\cdot,t_{0}),p^{\epsilon}(\cdot,t_{0})$ converge to
$p(\cdot,t_{0})$ in $H^{1}(\mathbb{R}^{d}).$ Since this must hold for almost
every $t_{0}\in[0,T]$ we are done. ∎

A straightforward consequence of the previous Lemma is the following Lemma
which gives a crude comparison between the pressure values at different times.
We will obtain a much sharper comparison property in Section 3 when we
establish the Hopf-Lax type formula for the pressure.

######  Lemma 2.6.

Fix some time $\tau>0$. Given almost any times $s,t\in[0,\tau]$ such that
$s<t$, there exists a constant $C(\tau)$ such that

$p(s,x)\leq C(\tau)p(t,x).$

In particular, this implies

$\\{x\in\mathbb{R}^{d}:w(s,x)>0\\}\subset\\{x\in\mathbb{R}^{d}:p(s,x)>0\\}\subset\\{x\in\mathbb{R}^{d}:w(t,x)>0\\}.$

###### Proof.

By Lemma 2.4, there exists a constant $C(\tau)>0$ such that
$n(x,s)<C(\tau)n(x,t)$ for all $x\in\mathbb{R}^{d}$. Since $\rho$ is
increasing with respect to time and $\rho\leq 1$, we know that
$\varphi(x)(1-\rho(x,t))=0$ for any nonnegative function $\varphi(x)$ such
that $\varphi(x)(1-\rho(x,s))=0$. Let us choose
$\varphi(x)=(p(s,x)-C(\tau)p(t,x))_{+}$. It then follows from Lemma 2.5 that

$\int_{\mathbb{R}^{d}}\nabla(p(x,s)-C(\tau)p(x,t))_{+}\cdot\nabla
p(x,t)\,dx=\int_{\mathbb{R}^{d}}(p(x,s)-C(\tau)p(x,t))_{+}n(x,t)\,dx,$

and

$\int_{\mathbb{R}^{d}}\nabla(p(x,s)-C(\tau)p(x,t))_{+}\cdot\nabla
p(x,s)\,dx=\int_{\mathbb{R}^{d}}(p(x,s)-C(\tau)p(x,t))_{+}n(x,s)\,dx.$

Hence,

$\int_{\mathbb{R}^{d}}\nabla(p(x,s)-C(\tau)p(x,t))_{+}\cdot\nabla(p(x,s)-C(\tau)p(x,t))\,dx=\int_{\mathbb{R}^{d}}(p(x,s)-C(\tau)p(x,t))_{+}(n(x,s)-C(\tau)n(x,t))\,dx.$

The left-hand side of the above equation is nonnegative while the right-hand
side of the equation is nonpositive. This is only possible if
$(p(x,s)-C(\tau)p(x,t))_{+}=0$ almost everywhere. ∎

######  Lemma 2.7.

Up to a set of measure zero, for any $t>0$ we have
$\\{x\in\mathbb{R}^{d}:\rho(x,t)=1\\}=\\{x\in\mathbb{R}^{d}:w(x,t)>0\\}$.

###### Proof.

This is nearly Lemma 4.6 of [JKT23], except that in the diffusion case we lack
an explicit formula for the nutrient. Nevertheless, we proceed along the same
lines.

From Lemma 2.3, we have $w(1-\rho)=0$, and thus

$\\{x\in\mathbb{R}^{d}:w(x,t)>0\\}\subset\\{x\in\mathbb{R}^{d}:\rho(x,t)=1\\}$

Thus, we must show that the set $A_{t}:=\\{x:\rho(x,t)=1,w(x,t)=0\\}$ has
measure zero.

For this, we observe that $\Delta w$ vanishes a.e. where $w$ vanishes, and
thus (2.3) implies that

$\rho(x,0)+\int_{0}^{t}n(x,s)\rho(x,s)\,ds=\rho(x,t)=1\hbox{ a.e. on }A_{t}$

From the pressure equation (2.2), any interior point of $\\{\rho(x,0)=1\\}$
has positive pressure at every positive time, and our assumptions on the
initial data provide that the boundary of this set has zero measure.

Thus, we need only consider the case where $\int_{0}^{t}n\rho\,ds=1$. Since
the nutrient is uniformly positive due to Lemma 2.4, this occurs for at most
one time for a given $x$. On the other hand, since the nutrient is uniformly
bounded, this function is continuous in time, and $x$ must be in $A_{t}$ for
an open set of times before $\int_{0}^{t}n\rho\,ds=1$ is satisfied. It follows
directly that $A_{t}$ (and, in fact, $\bigcup_{t}A_{t}$) is null.

∎

######  Lemma 2.8.

$T(x):=\inf\\{t\geq 0:w(x,t)>0\\}=\inf\\{t\geq 0:\rho(x,t)=1\\}\hbox{ for a.e.
}x\in{\mathbb{R}}^{d}$

###### Proof.

For conciseness, write $\widetilde{T}(x)=\inf\\{t\geq 0:\rho(x,t)=1\\}$.

Suppose that for some $x$, we have $\widetilde{T}(x)<T(x)$. Then we have
$w(x,t)=0$ for all $t<T(x)$. Then, since Lemma 2.2 gives that $\rho$ is
monotone in time except for a null set of $x$ which we ignore, we have
$x\in\\{y:w(y,t)=0,\rho(y,t)=1\\}$ for all $t\in(\widetilde{T}(x),T(x))$. Then
any such $x$ is contained in
$\bigcup_{t\in\mathbb{Q}\cap(0,\infty)}\\{y:w(y,t)=0,\rho(y,t)=1\\}$, which is
null by Lemma 2.7.

The other direction is similar. Suppose instead that for some $x$, we have
$T(x)<\widetilde{T}(x)$. Then $w(x,t)>0$ for all $t>T(x)$, while monotonicity
implies that for a.e. $x$ we have $\rho(x,t)=0$ for all $t<\widetilde{T}(x)$.
Then any such $x$ is contained in
$\bigcup_{t\in\mathbb{Q}\cap(0,\infty)}\\{y:w(y,t)>0,\rho(y,t)=0\\}$, which is
null by Lemma 2.7. ∎

######  Lemma 2.9.

$w(\cdot,t)$ solves the obstacle problem (1.8).

###### Proof.

This follows from applying Lemma 2.8 to Lemma 2.3. We have

$\rho(x,t)=\chi_{\\{T<t\\}}(x)=\chi_{\\{w(\cdot,t)>0\\}}(x)$

for a.e. $x$ and all $s\neq T(x)$. Thus, to obtain (1.8) we modify (2.3) by
replacing the occurrence of $\rho$ in $\int_{0}^{t}n\rho\,ds$ with
$\chi_{\\{T<t\\}}$ and the other occurrence of $\rho$ with
$\chi_{\\{w(\cdot,t)>0\\}}$. ∎

Recalling the notation of (1.11), we additionally define

(2.5) $\Omega_{\infty}:=\\{0\leq
T(x)<\infty\\}=\bigcup_{t>0}\Omega_{t},\quad\mathcal{O}:=\\{0<T(x)<\infty\\}=\Omega_{\infty}\setminus\overline{\Omega_{0}}.$

We now prove that the hitting time $T$ is continuous. While this justifies the
characterization of the level sets of $T$ as the free boundary $\Omega_{t}$,
it also is an important first step that initiates the regularity analysis of
$T$ in section 4. The main idea will be to show that a discontinuity must
result in a point $x_{0}$ and times $t_{0}<t_{1}$ such that $x_{0}$ is in
$\partial\Omega_{t}$ for $t_{0}\leq t\leq t_{1}$. Then
$w(\cdot,t_{1})-w(\cdot,t_{0})$ is a positive superharmonic function on
$\Omega_{t_{0}}$, so one would like to apply the Hopf lemma to draw a
contradiction between $w(x_{0},t_{0})=w(x_{0},t_{1})=0$ and $\nabla
w(x_{0},t_{0})=\nabla w(x_{0},t_{1})=0$. Unfortunately, $\Omega_{t_{0}}$ does
not a priori have the regularity needed to apply the Hopf lemma, so we must
first use obstacle problem techniques to shift to a setting where we do have
such regularity. The key tool in doing so will be the quadratic blowup of $w$
at free boundary points:

######  Lemma 2.10 ([Bla00] Corollary 2.5).

Let $u$ be a nonnegative solution to $\Delta u=f\chi_{\\{u>0\\}}$, for some
$f$ which is strictly positive and bounded near the free boundary
$\partial\\{u>0\\}$. Then if $x_{0}$ is a free boundary point, then the
quadratic blowup sequence $r^{-2}w(r(x-x_{0})+x_{0},t)$ is compact in
$C^{1,\alpha}(B_{1}(x_{0}))$ as $r\to 0^{+}$. Moreover, if $f$ is continuous
at $x_{0}$, then the subsequential limits solve $\Delta
v=f(x_{0})\chi_{\\{v>0\\}}$.

The subsequential limit enjoys better geometry, due to the following property
of global solutions to the constant-source obstacle problem:

######  Lemma 2.11 ([Caf98] Corollary 7).

A nonnegative solution to $\Delta u=\chi_{\\{u>0\\}}$ on ${\mathbb{R}}^{d}$ is
convex.

######  Proposition 2.12.

  * (a)

$T$ is continuous.

  * (b)

$x\in\partial\Omega_{t}$ if and only if $x\in\mathcal{O}$ and $t=T(x)$, for
all $x\in{\mathbb{R}}^{d}$ and $t>0$.

###### Proof.

First, we verify that the $\overline{\Omega_{t}}$ are continuous from above,
in the sense that for any $t$ we have:

(2.6)
$\overline{\Omega_{t}}=\bigcap_{\varepsilon>0}\overline{\Omega_{t+\varepsilon}}$

The forward inclusion is trivial by the monotonicity of $w$. For the reverse
inclusion, we suppose for contradiction that there exists
$x\in\bigcap_{\varepsilon>0}\overline{\Omega_{t+\varepsilon}}\setminus\overline{\Omega_{t}}$.
Let $r$ be sufficiently small that $B_{r}(x)\subset\\{w(\cdot,t)=0\\}$. From
(2.3), $\Delta w(\cdot,t+\varepsilon)\geq 1-\varepsilon\|n_{0}\|_{\infty}$ on
$B_{r}(x)\cap\\{w(\cdot,t+\varepsilon)>0\\}$, and by assumption we have
$x\in\overline{\Omega_{t+\varepsilon}}$ for all $\varepsilon>0$. It follows by
quadratic nondegeneracy for the obstacle problem (Lemma 6.1) that if
$\varepsilon<\frac{\|n_{0}\|_{\infty}}{2}$, then

$\sup_{B_{r}(x)}w(\cdot,t+\varepsilon)\geq Cr^{2}$

uniformly in $\varepsilon$. Since $w(\cdot,t)\equiv 0$ on $B_{r}(x)$, we get a
contradiction with the Lipschitz continuity of $w$ in time by shrinking
$\varepsilon$.

Now, we introduce

(2.7) $T_{0}(x):=\inf\\{t>0:x\in\overline{\Omega_{t}}\\}$

It is immediate that $T_{0}(x)\leq T(x)=\inf\\{t>0:x\in\Omega_{t}\\}$. We
claim that $x\in\partial\Omega_{t}$ if and only if $t\in[T_{0}(x),T(x)]$. It
is clear that for $t<T_{0}(x)$, $x\notin\overline{\Omega_{t}}$, and for
$t>T(x)$, $x\in\Omega_{t}$. Since $x\in\Omega_{T(x)+\varepsilon}$ for every
$\varepsilon>0$, (2.6) gives $x\in\overline{\Omega_{T(x)}}$. On the other
hand, by continuity of $w$ and minimality of $T$, we have
$x\notin\Omega_{T(x)}$, so $x\in\partial\Omega_{T(x)}$. Monotonicity implies
that $x$ is a boundary point for all $t\leq T(x)$ for which
$x\in\overline{\Omega_{t}}$. We have $x\in\overline{\Omega_{T_{0}(x)}}$ by
using (2.6) with the definition of $T_{0}$, so we get the claim.

For purely topological reasons related to how each is defined, $T_{0}$ is
lower semicontinuous and $T$ is upper semicontinuous. To check lower
semicontinuity of $T_{0}$, let $(x_{n})$ be a sequence converging to $x$ with
$\liminf T_{0}(x_{n}):=t$. Then for any $\varepsilon>0$, the $x_{n}$ are
eventually in $\overline{\Omega_{t+\varepsilon}}$, and thus
$x\in\overline{\Omega_{t+\varepsilon}}$. It follows that $T_{0}(x)\leq t$. To
check upper semicontinuity of $T$, we note that
$x\in\Omega_{T(x)+\varepsilon}$ for all $\varepsilon>0$. Any sequence $x_{n}$
converging to $x$ eventually has $T(x_{n})\leq T(x)+\varepsilon$ since
$\Omega_{T(x)+\varepsilon}$ is open, so we conclude that $\limsup T(x_{n})\leq
T(x)$.

Therefore, both parts of the proposition will follow if we can show that
$T_{0}\equiv T$. For this, we will need the following useful property:

(2.8) $\lim_{\Omega_{T_{0}(x)}\ni x_{n}\to x}T(x_{n})=T_{0}(x)$

To see this, we note that if $(x_{n})$ is such a sequence, then $T(x_{n})\leq
T_{0}(x)$ for each $n$. On the other hand, by minimality of $T_{0}$,
$x\notin\overline{\Omega_{T_{0}(x)-\varepsilon}}$, and so we eventually have
$T(x_{n})\geq T_{0}(x)-\varepsilon$ for any $\varepsilon>0$.

Finally, we proceed to the proof that $T_{0}=T$. Suppose that
$x_{0}\in\partial\Omega_{t}$ with $T_{0}(x_{0})<T(x_{0})$. We will use the
obstacle problem theory to compare blowups of $w$ at $(x_{0},T_{0}(x_{0}))$
and at $(x_{0},t)$ with $t>T_{0}(x_{0})$ to derive a contradiction. First let
us ensure that the blow-up profiles are well-defined. Due to (1.8) it follows
that

$0\leq\eta(x,t)\leq(t-T(x))_{+}\|n_{0}\|_{\infty}$

We also get a continuity estimate. Assuming $T(y)\leq T(x)$, we have either
$T(y)\leq t\leq T(x)$, giving

$|\eta(x,t)-\eta(y,t)|=\eta(x,t)\leq\|n_{0}\|_{\infty}|T(x)-T(y)|$

or else $t\leq T(y)\leq T(x)$, giving

$\displaystyle|\eta(x,t)-\eta(y,t)|$
$\displaystyle=\left|\int_{T(x)}^{t}n(x,s)-n(y,s)\,ds-\int_{T(y)}^{T(x)}n(y,s)\,ds\right|$
$\displaystyle\leq|t-T(x)|\sup_{T(x)\leq s\leq
t}|n(x,s)-n(y,s)|+\|n_{0}\|_{\infty}|T(x)-T(y)|$

Thus, using the nutrient regularity from Lemma 2.2 and the result of (2.8), we
conclude that $\eta(\cdot,T_{0}(x_{0}))$ restricted to $\Omega_{T_{0}(x_{0})}$
is continuous at $x_{0}$, and in a sufficiently small neighborhood of $x_{0}$,
we can ensure that it is less than $\frac{1}{2}$. Then Lemma 2.10 gives that
the family of rescalings $x\mapsto r^{-2}w(r(x-x_{0})+x_{0},T_{0}(x_{0}))$ are
compact as $r\to 0$ in $C^{1,\alpha}_{loc}$, and their subsequential limits
are nonzero global solutions of

(2.9) $\Delta u=\chi_{\\{u>0\\}}.$

Now, choose $\tau\in(T_{0}(x),T(x))$, sufficiently small such that we still
have $\eta(x,\tau)<\frac{1}{2}$ in some neighborhood of $x_{0}$. By taking a
further subsequence, the discussion above yields a sequence $r_{n}\to 0$ such
that

$r_{n}^{-2}w(r_{n}(x-x_{0})+x_{0},T(x_{0}))\to u\hbox{ and
}r_{n}^{-2}w(r_{n}(x-x_{0})+x_{0},\tau)\to v,$

for some $u,v$ in $C^{1,\alpha}_{loc}({\mathbb{R}}^{d})$. Unlike with $u$,
$\eta(\cdot,\tau)$ restricted to $\Omega_{\tau}$ is not known to be continuous
at $x_{0}$, since we do not yet know that $T$ is continuous. In particular, we
do not know that $v$ solves a constant Laplacian obstacle problem. However, we
do get $u(x_{0})=v(x_{0})=0$ and $\nabla u(x_{0})=\nabla v(x_{0})=0$ from the
convergence, and we also have that $v\geq u$ since $w(\cdot,t)\geq
w(\cdot,T_{0}(x_{0}))$.

We will apply the Hopf lemma to $v-u$ in the domain $U:=\\{u>0\\}$. First
observe that from the definition of $u$, we have
$r_{n}(x-x_{0})+x_{0}\in\Omega_{T(x_{0})}$ if $x\in U$ and if $n$ is
sufficiently large depending on $x$. We have checked above that
$\eta(\cdot,\tau)$ restricted to $\Omega_{T_{0}(x_{0})}$ is continuous at
$x_{0}$. Thus, for any $x\in U$, we have
$\eta(r_{n}(x-x_{0})+x_{0},\tau)\to\eta(x_{0},\tau)$, and so we conclude that

$\Delta v=1-\eta(x_{0},\tau)\hbox{ in }U.$

Comparing this equation to (2.9), it follows that $v-u$ satisfies

$\Delta(v-u)=-\eta(x_{0},\tau)\leq-(\tau-
T_{0}(x_{0}))\inf_{t\in[T_{0}(x_{0}),\tau]}n(x_{0},t)<0$

with the last inequality following from the nutrient lower bound in Lemma 2.4
and the assumption that the initial nutrient is bounded away from 0. This
implies that $v-u$ is strictly superharmonic inside $U$, and so our previous
observation that $v-u\geq 0$ by the monotonicity in time of $w$ improves to
$v-u>0$ inside $U$. Lastly let us observe that, from Lemma 2.11, the
complement of $U$ is convex and so $U$ satisfies the interior ball condition
at $x_{0}$. Putting together the above information, the Hopf lemma applied at
$x_{0}$ implies that $\nabla v(x_{0})-\nabla u(x_{0})\neq 0$, which is a
contradiction. It follows that $T_{0}=T$, so we finish. ∎

######  Remark 2.13.

An important consequence of Proposition 2.12 is that the spacetime interface
is exactly the graph of $T$ on $\mathcal{O}$. In other words,

(2.10)
$\\{(x,t):t\in(0.\infty),x\in\partial\Omega_{t}\\}=\mathrm{Graph}_{T}(\mathcal{O}):=\\{(x,T(x)):x\in\mathcal{O}\\}$

This also means that the interface is a $d$-dimensional topological manifold,
and the regularity of its parametrization in $d$ spatial variables is exactly
that of $T$. We will use the notation $\mathrm{Graph}_{T}$ with subsets of
$\mathcal{O}$, which may be understood in this light as projections of the
spacetime interface into ${\mathbb{R}}^{d}$.

Finally, in light of the regularity of $w$ and $T$, we note a natural way to
standardize $\rho$ on measure zero sets. A corresponding standardization of
the pressure will need to wait until the next Section, due to the need to
preserve certain delicate structures.

######  Lemma 2.14.

$\partial\Omega_{t}$ has zero measure in ${\mathbb{R}}^{d}$ for all $t>0$. The
weak solution $(\rho,p,n)$ can be taken such that $\rho$ is upper
semicontinuous in space and time, with
$\\{\rho(\cdot,t)=1\\}=\overline{\Omega_{t}}$ for each $t$. In particular, the
support of $p$ is then contained in $\\{\rho=1\\}$.

###### Proof.

By Proposition 2.12, for any $t>0$ and $\varepsilon\in(0,t)$, we have
$\partial\Omega_{t}\subset\Omega_{t+\varepsilon}\setminus\Omega_{t-\varepsilon}$.
By Lemma 2.7, up to measure zero sets we can replace the right-hand-side with
$\\{\rho(\cdot,t+\varepsilon)=1\\}\setminus\\{\rho(\cdot,t-\varepsilon)=1\\}$,
and by time continuity of $\rho$ in $L^{1}$, the measure of this set goes to 0
with $\varepsilon$. Thus, $\partial\Omega_{t}$ has zero measure.
Then we claim that $(\rho,p,n)$ with $\rho$ redefined as
$\chi_{\overline{\\{w>0\\}}}$ and $p$ redefined to vanish outside
$\overline{\\{w>0\\}}$ remains a weak solution as defined in Definition 2.1.
Indeed, since this changes $\rho,p$ by measure zero sets for each time,
equation (2.1) is unaffected. On the other hand, we have $p(1-\rho)=0$ by
construction. To check that $\rho$ is spacetime upper semicontinuous, we only
need to show that the set $\\{\rho=1\\}:=\\{(x,t):t\geq
0,x\in\overline{\Omega_{t}}\\}$ is closed. Let $(x_{n},t_{n})$ be a sequence
in this set converging to $(x,t)$. Since the $t_{n}$ converge, they are
bounded, and so the sequence is contained in an $\overline{\Omega_{\tau}}$ for
$\tau$ sufficiently large. Then $T(x)<\infty$, so we have $T(x_{n})\to T(x)$
by continuity, and since $T(x_{n})\leq t_{n}$ for each $n$ by Proposition
2.12, we have $T(x)\leq t$. It follows that $x\in\overline{\Omega_{t}}$, so we
conclude. ∎

##  3\. AB estimates and the Hopf-Lax bound

In this section, we will show that there exists a nonnegative function $u_{+}$
such that the pressure is a super solution to the following Hamilton-Jacobi
equation

(3.1) $\partial_{t}p-|\nabla p|^{2}\geq pu_{+}.$

We will then use (3.1) to obtain a Hopf-Lax type formula for the pressure. In
particular, given a fixed time $t_{0}$, this will allow us to give lower
bounds for the pressure at times $t>t_{0}$ and upper bounds for the pressure
at times $t<t_{0}$ in terms of $p(t_{0},\cdot)$. This will give us a very
precise way of constructing pressure super solutions that lead to powerful
barrier-type arguments and eventually Hölder regularity of the hitting times
(c.f. Section 4).

Let us emphasize that to the best of our knowledge, the Hopf-Lax type bounds
we obtain have not previously appeared in the literature for Hele-Shaw type
equations and they require some highly nontrivial efforts to obtain. First, to
establish (3.1), we go through the Porous Media Equation (PME) and use the
fact that our solution $(\rho,p)$ can be obtained as the incompressible limit
of solutions $(\rho_{\gamma},p_{\gamma},n_{\gamma})$ of the PME-nutrient
system

(3.2) $\partial_{t}\rho_{\gamma}-\nabla\cdot(\rho_{\gamma}\nabla
p_{\gamma})=\rho_{\gamma}n_{\gamma},\quad p_{\gamma}=\rho_{\gamma}^{\gamma},$
(3.3) $\partial_{t}n_{\gamma}-\Delta n_{\gamma}=-\rho_{\gamma}n_{\gamma}$

as the scalar parameter $\gamma$ is sent to infinity, and where the nutrient
variable from our original system is held fixed. The advantage of the PME
system is that it is possible to use the relation
$p_{\gamma}=\rho_{\gamma}^{\gamma}$ to rewrite (3.2) solely in terms of the
pressure variable $p_{\gamma}$, which yields the equation

(3.4) $\partial_{t}p_{\gamma}-|\nabla p_{\gamma}|^{2}-\gamma p_{\gamma}(\Delta
p_{\gamma}+n)=0.$

The main difficulty in obtaining (3.1) is to show that as $\gamma\to\infty$,
$u_{\gamma}=-\gamma(\Delta p_{\gamma}+n)$ converges to a meaningful limit
object $u$, whose positive part can be controlled. For the classic PME without
a source term, bounds on the negative part of $\Delta p_{\gamma}$ are known
through the celebrated Aronson-Benilan estimate [AB79]. In the presence of a
source term, AB-type bounds on quantities taking a similar form to
$\gamma(\Delta p_{\gamma}+n)$ have been studied in the literature [PQV14,
GPŚG19, DP21, Jac23], however except for [Jac23], these bounds do not scale
well with respect to $\gamma$. We adapt the arguments from [Jac23] to show
that $[u_{\gamma}]_{+}$ can be bounded uniformly with respect to $\gamma$ in
BMO-type spaces. Note that we are unable to get $L^{\infty}$ bounds on
$u_{\gamma}$ essentially because two key quantities in the estimate
$\partial_{t}n$ and $\nabla n\cdot\nabla p$ are not in general bounded in
$L^{\infty}$. It would be interesting to see whether equation (3.1) could be
obtained directly from the original system without going through PME, but we
leave this question to a future work.

Once we have obtained equation (3.1), there is still significant work required
to obtain a Hopf-Lax type control for $p$. Here the difficulty is that $u_{+}$
is not bounded in $L^{\infty}$. Noting that $p$ should satisfy

$\partial_{t}p-|\nabla p|^{2}+pu_{+}\geq 0,$

the derivative of $p$ along an arbitrary path $x(t)$ gives

$\frac{d}{dt}p(x(t),t)=\partial_{t}p(x(t),t)+x^{\prime}(t)\cdot\nabla
p(x(t),t)\geq$ $|\nabla
p(x(t),t)|^{2}-p(x(t),t)u_{+}(x(t),t)+x^{\prime}(t)\cdot\nabla
p(x(t),t)\geq-p(x(t),t)u_{+}(x(t),t)-\frac{1}{4}|x^{\prime}(t)|^{2}.$

Unfortunately, without $L^{\infty}$ control on $u_{+}$ it is not clear that
time integrals of the final quantity will be well-defined. This prevents the
usual approach to proving Hopf-Lax type bounds.

To overcome this, we adapt the approach developed in [CG15], which handles
unbounded coefficients by instead considering an average over paths indexed by
the unit ball. Our calculation is somewhat different however, as we can
exploit the specific structure of $pu_{+}$ to write $pu_{+}=\lambda
p+p(u-\lambda)_{+}$ for some scalar $\lambda\geq 0$. By choosing $\lambda$
appropriately we can force $p(u-\lambda)_{+}$ to be small while using a
Gronwall argument to handle $\lambda p$. This allows us to obtain a much more
favorable error term in our Hopf-Lax formula compared to [CG15] (c.f.
Proposition 3.7).

We begin with the aforementioned result (well-known) that says we can
approximate our system (1.1-1.2) with a sequence of smooth solutions to PME.

######  Proposition 3.1 (see e.g. [PQV14, GPŚG19, Jac23]).

There exists a sequence of smooth solutions
$(\rho_{\gamma},p_{\gamma},n_{\gamma})$ to the PME-nutrient system (3.2-3.3)
with initial data $(\rho_{0,\gamma},n_{0,\gamma})$ such that for any $\tau>0$
we have that $\rho_{\gamma}$ converges strongly in
$L^{1}([0,\tau]\times\mathbb{R}^{d})$, $p_{\gamma},n_{\gamma}$ converge
strongly in $L^{2}([0,\tau];H^{1}(\mathbb{R}^{d}))$ to the unique solution
$(\rho,p,n)$ to the system (1.1-1.2) with initial data $(\rho_{0},n_{0})$ as
$\gamma\to\infty$. Furthermore, one may choose $\rho_{0,\gamma}$ such that
$u_{\gamma,+}(\cdot,0)$ is bounded in $L^{\infty}(\mathbb{R}^{d})$ uniformly
in $\gamma$.

Next, we record the following simple result for solutions to the heat equation
with $L^{\infty}$ source.

######  Lemma 3.2.

There exists some $b_{0}>0$ such that $n$ satisfies the bound

$\exp(b_{0}\frac{|\partial_{t}n|+|\Delta n|}{n})-1\in
L^{1}([0,\tau]\times\mathbb{R}^{d})$

###### Proof.

Thanks to Lemma 2.2 we know that $\partial_{t}n$ and $D^{2}n$ are bounded in
BMO. Since we also have

$\frac{1}{2}\lVert\partial_{t}n\rVert_{L^{2}([0,\tau]\times\mathbb{R}^{d})}^{2}+\lVert\nabla
n\rVert_{L^{2}(\\{\tau\\}\times\mathbb{R}^{d})}^{2}\leq\lVert\nabla
n\rVert_{L^{2}(\\{0\\}\times\mathbb{R}^{d})}^{2}+\frac{1}{2}\lVert
n\rVert_{L^{\infty}([0,\tau]\times\mathbb{R}^{d})}\lVert\rho\rVert_{L^{2}([0,\tau]\times\mathbb{R}^{d})},$

the BMO bound implies the existence of a constant $c>0$ such that
$\exp(c|\partial_{t}n|)-1,\exp(c|\Delta n|)-1\in
L^{1}([0,T]\times\mathbb{R}^{d})$. Following the logic of Lemma 2.4, it
follows that $n$ is uniformly bounded from below on any time interval. Hence,
there must be an appropriate choice of $b_{0}$ where the result holds. ∎

######  Proposition 3.3.

If $(p_{\gamma},n_{\gamma})$ is a smooth solution to the system (3.3-3.4) for
some $\gamma\in(1,\infty)$, then for any $\tau>0$ there exists $b>0$ that only
depends on $\tau$ such that $(b[u_{\gamma}]_{+}-1)\exp(b[u_{\gamma}]_{+})+1$
is uniformly bounded in $L^{1}([0,\tau]\times\mathbb{R}^{d})$ with respect to
$\gamma$, where $u_{\gamma}:=-\gamma(\Delta p_{\gamma}+n_{\gamma})$.

###### Proof.

If we differentiate $\frac{1}{\gamma}u_{\gamma}$ with respect to time, we get

$\partial_{t}\frac{1}{\gamma}u_{\gamma}=-\partial_{t}n_{\gamma}-\Delta\partial_{t}p_{\gamma}=-\partial_{t}n_{\gamma}-\Delta(|\nabla
p|^{2}-p_{\gamma}u_{\gamma})$

Expanding the Laplacian, we see that

$\partial_{t}\frac{1}{\gamma}u_{\gamma}=-\partial_{t}n_{\gamma}-2|D^{2}p_{\gamma}|^{2}-2\nabla\Delta
p_{\gamma}\cdot\nabla p_{\gamma}+2\nabla p_{\gamma}\cdot\nabla
u_{\gamma}+p_{\gamma}\Delta u_{\gamma}+u_{\gamma}\Delta p_{\gamma}.$

Noting that $-\Delta p_{\gamma}=n_{\gamma}+\frac{1}{\gamma}u_{\gamma}$ we can
rewrite the previous line as

$\partial_{t}\frac{1}{\gamma}u_{\gamma}=2\nabla n_{\gamma}\cdot\nabla
p_{\gamma}-\partial_{t}n_{\gamma}-2|D^{2}p_{\gamma}|^{2}+2(1+\frac{1}{\gamma})\nabla
p_{\gamma}\cdot\nabla u_{\gamma}+p_{\gamma}\Delta
u_{\gamma}-n_{\gamma}u_{\gamma}-\frac{1}{\gamma}u_{\gamma}^{2}.$

Hence, eliminating $|D^{2}p|^{2}$, we can conclude that

(3.5)
$n_{\gamma}u_{\gamma}+\frac{1}{\gamma}(\partial_{t}u_{\gamma}+u_{\gamma}^{2})\leq
2\nabla n_{\gamma}\cdot\nabla
p_{\gamma}-\partial_{t}n_{\gamma}+2(1+\frac{1}{\gamma})\nabla
p_{\gamma}\cdot\nabla u_{\gamma}+p_{\gamma}\Delta u_{\gamma}.$

Now let $f:\mathbb{R}\to\mathbb{R}$ be a $C^{2}$ convex function such that
$f^{\prime}\geq 0$ everywhere and $f=0$ on $(-\infty,0]$. If we integrate
(3.5) against $f^{\prime}(u_{\gamma})$ on $[0,\tau]\times\mathbb{R}^{d}$ we
find that

(3.6)
$\int_{\mathbb{R}^{d}\times\\{\tau\\}}\frac{1}{\gamma}f(u_{\gamma})+\int_{\mathbb{R}^{d}\times[0,\tau]}n_{\gamma}u_{\gamma}f^{\prime}(u_{\gamma})+\frac{1}{\gamma}u_{\gamma}^{2}f^{\prime}(u_{\gamma})\leq\\\
\int_{\mathbb{R}^{d}\times\\{0\\}}\frac{1}{\gamma}f(u_{\gamma})+\int_{\mathbb{R}^{d}\times[0,\tau]}f^{\prime}(u_{\gamma})\big{(}2\nabla
n_{\gamma}\cdot\nabla
p_{\gamma}-\partial_{t}n_{\gamma}\big{)}+2(1+\frac{1}{\gamma})\nabla
p_{\gamma}\cdot\nabla\big{(}f(u_{\gamma})\big{)}+p_{\gamma}f^{\prime}(u_{\gamma})\Delta
u_{\gamma}.$

Noting that $\Delta\big{(}f(u_{\gamma})\big{)}=f^{\prime}(u_{\gamma})\Delta
u_{\gamma}+f^{\prime\prime}(u_{\gamma})|\nabla u_{\gamma}|^{2}$, we can
integrate by parts in (3.6) to obtain

(3.7)
$\int_{\mathbb{R}^{d}\times\\{\tau\\}}\frac{1}{\gamma}f(u_{\gamma})+\int_{\mathbb{R}^{d}\times[0,\tau]}n_{\gamma}u_{\gamma}f^{\prime}(u_{\gamma})+\frac{1}{\gamma}u_{\gamma}^{2}f^{\prime}(u_{\gamma})+p_{\gamma}f^{\prime\prime}(u_{\gamma})|\nabla
u_{\gamma}|^{2}\leq\\\
\int_{\mathbb{R}^{d}\times\\{0\\}}\frac{1}{\gamma}f(u_{\gamma})+\int_{\mathbb{R}^{d}\times[0,\tau]}f^{\prime}(u_{\gamma})\big{(}2\nabla
n_{\gamma}\cdot\nabla
p_{\gamma}-\partial_{t}n_{\gamma}\big{)}-(1+\frac{2}{\gamma})f(u_{\gamma})\Delta
p_{\gamma}$

We then integrate by parts in $\nabla n_{\gamma}\cdot\nabla p_{\gamma}$ to get

(3.8)
$\int_{\mathbb{R}^{d}\times\\{\tau\\}}\frac{1}{\gamma}f(u_{\gamma})+\int_{\mathbb{R}^{d}\times[0,\tau]}n_{\gamma}u_{\gamma}f^{\prime}(u_{\gamma})+\frac{1}{\gamma}u_{\gamma}^{2}f^{\prime}(u_{\gamma})+p_{\gamma}f^{\prime\prime}(u_{\gamma})|\nabla
u_{\gamma}|^{2}\leq\\\
\int_{\mathbb{R}^{d}\times\\{0\\}}\frac{1}{\gamma}f(u_{\gamma})-\int_{\mathbb{R}^{d}\times[0,\tau]}f^{\prime}(u_{\gamma})\big{(}2p_{\gamma}\Delta
n_{\gamma}+\partial_{t}n_{\gamma}\big{)}+p_{\gamma}f^{\prime\prime}(u_{\gamma})\nabla
u_{\gamma}\cdot\nabla n_{\gamma}+(1+\frac{2}{\gamma})f(u_{\gamma})\Delta
p_{\gamma}.$

Once again using $-\Delta p_{\gamma}=n_{\gamma}+\frac{1}{\gamma}u_{\gamma}$
and using the quadratic Young’s inequality on $\nabla u_{\gamma}\cdot\nabla
n_{\gamma}$ we find that

(3.9)
$\int_{\mathbb{R}^{d}\times\\{\tau\\}}\frac{1}{\gamma}f(u_{\gamma})+\int_{\mathbb{R}^{d}\times[0,\tau]}n_{\gamma}u_{\gamma}f^{\prime}(u_{\gamma})+\frac{1}{\gamma}u_{\gamma}^{2}f^{\prime}(u_{\gamma})+\frac{1}{2}p_{\gamma}f^{\prime\prime}(u_{\gamma})|\nabla
u_{\gamma}|^{2}\leq\\\
\int_{\mathbb{R}^{d}\times\\{0\\}}\frac{1}{\gamma}f(u_{\gamma})+\int_{\mathbb{R}^{d}\times[0,\tau]}(1+\frac{2}{\gamma})f(u_{\gamma})(n_{\gamma}+\frac{1}{\gamma}u_{\gamma})-f^{\prime}(u_{\gamma})\big{(}2p_{\gamma}\Delta
n_{\gamma}+\partial_{t}n_{\gamma}\big{)}+\frac{1}{2}p_{\gamma}f^{\prime\prime}(u_{\gamma})|\nabla
n_{\gamma}|^{2}.$

Next, to help compare the left and right hand sides, we divide and multiply by
multiples of $n_{\gamma}$ to get

(3.10)
$\int_{\mathbb{R}^{d}\times\\{\tau\\}}\frac{1}{\gamma}f(u_{\gamma})+\int_{\mathbb{R}^{d}\times[0,\tau]}n_{\gamma}u_{\gamma}f^{\prime}(u_{\gamma})+\frac{1}{\gamma}u_{\gamma}^{2}f^{\prime}(u_{\gamma})+\frac{1}{2}p_{\gamma}f^{\prime\prime}(u_{\gamma})|\nabla
u_{\gamma}|^{2}\leq\\\
\int_{\mathbb{R}^{d}\times\\{0\\}}\frac{1}{\gamma}f(u_{\gamma})+\int_{\mathbb{R}^{d}\times[0,\tau]}(1+\frac{2}{\gamma})f(u_{\gamma})(n_{\gamma}+\frac{1}{\gamma}u_{\gamma})-\frac{n_{\gamma}}{2}f^{\prime}(u_{\gamma})\big{(}\frac{4p_{\gamma}\Delta
n_{\gamma}+2\partial_{t}n_{\gamma}}{n_{\gamma}}\big{)}+n_{\gamma}f^{\prime\prime}(u_{\gamma})\frac{p_{\gamma}|\nabla
n_{\gamma}|^{2}}{2n_{\gamma}}.$

Using the identity $uf^{\prime}(u)-f(u)=f^{*}(f^{\prime}(u))$ and applying
Young’s inequality to
$-\frac{n_{\gamma}}{2}f^{\prime}(u_{\gamma})\big{(}\frac{2p_{\gamma}\Delta
n_{\gamma}+\partial_{t}n_{\gamma}}{n_{\gamma}}\big{)}$, we get

(3.11)
$\int_{\mathbb{R}^{d}\times\\{\tau\\}}\frac{1}{\gamma}f(u_{\gamma})+\int_{\mathbb{R}^{d}\times[0,\tau]}(n_{\gamma}+\frac{1}{\gamma}u_{\gamma})f^{*}(f^{\prime}(u_{\gamma}))+\frac{1}{2}p_{\gamma}f^{\prime\prime}(u_{\gamma})|\nabla
u_{\gamma}|^{2}\leq\\\
\int_{\mathbb{R}^{d}\times\\{0\\}}\frac{1}{\gamma}f(u_{\gamma})+\int_{\mathbb{R}^{d}\times[0,\tau]}\frac{2}{\gamma}f(u_{\gamma})(n_{\gamma}+\frac{1}{\gamma}u_{\gamma})+\frac{n_{\gamma}}{2}f^{*}(f^{\prime}(u_{\gamma}))+\frac{n_{\gamma}}{2}f\big{(}\frac{4p_{\gamma}\Delta
n_{\gamma}+2\partial_{t}n_{\gamma}}{n_{\gamma}}\big{)}+n_{\gamma}f^{\prime\prime}(u_{\gamma})\frac{p_{\gamma}|\nabla
n_{\gamma}|^{2}}{2n_{\gamma}},$

which is finally in a form that will allow us to estimate.

Fix some $b\leq\frac{b_{0}}{4\max\big{(}1,\sup_{\gamma}\lVert
p_{\gamma}\rVert_{L^{\infty}([0,\tau]\times\mathbb{R}^{d})}\big{)}}$ where
$b_{0}$ is the constant from Lemma 3.2. If we choose $f$ such that $f$ grows
like $\exp(bu)$ at infinity, then $f^{*}(f^{\prime}(u))$ grows like
$bu\exp(bu)$ at infinity, and hence $f^{*}(f^{\prime}(u))$ dominates both
$f(u)$ and $f^{\prime\prime}(u)$ at infinity. PME has finite propagation in
time (uniform in $\gamma$) [Váz07], thus, there exists a radius $R=R_{\tau}>0$
sufficiently large such that $(\rho_{\gamma},p_{\gamma})$ is supported in
$B_{R}$ independently of $\gamma$. Recalling that $u_{\gamma}=-\gamma(\Delta
p_{\gamma}+n_{\gamma})$ and $n_{\gamma}\geq 0$, it follows that
$f(u_{\gamma}),f^{\prime}(u_{\gamma}),f^{\prime\prime}(u_{\gamma})$ are all
supported on $B_{R}$ independently of $\gamma$ as well. Since we are
integrating functions with uniformly bounded support and
$f\big{(}\frac{4p_{\gamma}\Delta
n_{\gamma}+2\partial_{t}n_{\gamma}}{n_{\gamma}}\big{)}$ is bounded by Lemma
3.2 for our choice of $b$, it follows that the left-hand side of (3.11)
dominates the right-hand side and so the result follows.

∎

It essentially immediately follows that $p$ is a weak supersolution to the
appropriate HJB equation.

######  Corollary 3.4.

Given any
$L^{2}_{\operatorname{\textup{loc}}}([0,\infty);L^{2}(\mathbb{R}^{d}))$ weak
limit point $u_{+}$ of the family $u_{\gamma,+}$ $p$ solves, in the sense of
weak solutions,

(3.12) $\partial_{t}p-|\nabla p|^{2}+u_{+}p\geq 0,$

where for any $\tau>0$ there exists $b=b(\tau,d)>0$ such that
$(bu_{+}-1)e^{bu_{+}}+1\in L^{1}([0,\tau];{\mathbb{R}}^{d})$.

Although we now know that $p$ is a supersolution to an HJB equation, it is
somewhat annoying to directly obtain the Hopf-Lax formula from (3.12), due to
the fact that $p$ is not continuous. Instead, we will work towards the Hopf-
Lax formula by once again going through the $\gamma$ limit. Here, we will
still need to deal with the difficulty that $u_{\gamma,+}$ is not uniformly
bounded in $L^{\infty}$. We proceed by adapting an argument from [CG15], which
provides a method to obtain Hopf-Lax type formulas for Hamilton-Jacobi
equations with unbounded coefficients. A key difference in our setting is that
the right-hand side has the specific form $p_{\gamma}u_{\gamma,+}$. This
structure allows us to combine their approach with Gronwall-type estimates to
obtain much stronger bounds.

######  Lemma 3.5.

Choose a decreasing nonnegative function $\lambda\in L^{1}([0,t_{1}-t_{0}])$.
Given any $\gamma\in(1,\infty)$ and any points $(x_{1},t_{1}),(x_{0},t_{0})$
with $t_{0}<t_{1}$ there exists a constant $C=C(t_{1},d)$ such that

(3.13) $p_{\gamma}(x_{0},t_{0})\leq
e^{\Lambda_{\gamma}(t_{1}-t_{0})}\Big{(}p_{\gamma}(x_{1},t_{1})+\frac{|x_{1}-x_{0}|^{2}}{4\int_{0}^{t_{1}-t_{0}}e^{\Lambda_{\gamma}(s)}\,ds}+C(t_{1}-t_{0})^{7/10}e^{-\lambda(t_{1}-t_{0})}\Big{)}$

where $b$ is the constant from Proposition 3.3 and

$\Lambda_{\gamma}(t):=\frac{5}{4b}\int_{0}^{t}\lambda(a)\,da+\frac{1}{b}\int_{0}^{t}\log(1+\lVert\exp(bu_{\gamma,+})-1\rVert_{L^{1}(\\{t_{1}-a\\}\times\mathbb{R}^{d})})\,da.$

###### Proof.

Define $\varphi_{\gamma}(x,t):=p_{\gamma}(x,t_{1}-t)$. It then follows that
$\varphi_{\gamma}$ satisfies the differential inequality

$\partial_{t}\varphi_{\gamma}(x,t)+|\nabla\varphi_{\gamma}(x,t)|^{2}\leq\varphi_{\gamma}(x,t)u_{\gamma,+}(x,t_{1}-t)$

almost everywhere. Define

$\bar{\lambda}_{\gamma}(s):=\lambda(s)+\frac{1}{b}\log(\lVert\exp(bu_{\gamma,+})-1\rVert_{L^{1}(\\{a\\}\times\mathbb{R}^{d})})$

and split

$\varphi_{\gamma}(x,t)u_{\gamma,+}(x,t_{1}-t)\leq\varphi_{\gamma}(x,t)\bar{\lambda}_{\gamma}(t)+\varphi_{\gamma}(x,t)(u_{\gamma}(x,t_{1}-t)-\bar{\lambda}_{\gamma}(t))_{+}$

Multiplying both sides of the differential inequality by
$e^{-\Lambda_{\gamma}(t)}$ we see that

$\partial_{t}(e^{-\Lambda_{\gamma}(t)}\varphi_{\gamma})+e^{-\Lambda_{\gamma}(t)}|\nabla\varphi_{\gamma}|^{2}\leq
e^{-\Lambda_{\gamma}(t)}\varphi_{\gamma}(u_{\gamma}-\bar{\lambda}_{\gamma}(t))_{+}.$

Let $q_{\gamma}(t,x):=e^{-\Lambda_{\gamma}(t)}\varphi_{\gamma}(t,x)$, we then
have

$\partial_{t}q_{\gamma}+e^{\Lambda_{\gamma}(t)}|\nabla q_{\gamma}|^{2}\leq
q_{\gamma}(u_{\gamma}-\bar{\lambda}_{\gamma}(t))_{+}.$

Fix any two points $x_{1},x_{0}\in\mathbb{R}^{d}$. We now introduce a family
of paths $x_{\sigma}$ in the spirit of the path optimization argument
introduced in [CG15]. For each $\sigma$ in the unit ball $B_{1}$ let
$x_{\sigma}:[0,t_{1}-t_{0}]\to\mathbb{R}^{d}$ be a path such that
$x_{\sigma}(0)=x_{1}$ and $x_{\sigma}(t_{1}-t_{0})=x_{0}$. Consider

(3.14)
$\frac{d}{dt}\big{[}q_{\gamma}(x_{\sigma}(t),t)-\frac{1}{4}\int_{t_{0}}^{t}e^{-\Lambda_{\gamma}(s)}|x_{\sigma}^{\prime}(s)|^{2}\,ds\big{]}=\\\
\partial_{t}q_{\gamma}(x_{\sigma}(t),t)+\nabla
q_{\gamma}(x_{\sigma}(t),t)\cdot
x_{\sigma}^{\prime}(t)-\frac{e^{-\Lambda_{\gamma}(t)}}{4}|x_{\sigma}^{\prime}(t)|^{2}\leq
q_{\gamma}(x_{\sigma}(t),t)(u_{\gamma}(,x_{\sigma}(t)t_{1}-t)-\bar{\lambda}_{\gamma}(t))_{+}$

Thus,

$q_{\gamma}(x_{0},t_{1}-t_{0})\leq
q_{\gamma}(x_{1},0)+\frac{1}{4}\int_{0}^{t_{1}-t_{0}}e^{-\Lambda_{\gamma}(s)}|x_{\sigma}^{\prime}(s)|^{2}+q_{\gamma}(x_{\sigma}(s),s)(u_{\gamma}(x_{\sigma}(s),t_{1}-s)-\bar{\lambda}_{\gamma}(s))_{+}\,ds,$

It then follows that

(3.15)
$\varphi_{\gamma}(x_{0},t_{1}-t_{0})e^{-\Lambda_{\gamma}(t_{1}-t_{0})}\leq\varphi_{\gamma}(x_{1},0)+\frac{1}{4}\int_{0}^{t_{1}-t_{0}}e^{-\Lambda_{\gamma}(s)}\Big{(}|x_{\sigma}^{\prime}(s)|^{2}+\varphi_{\gamma}(s,x_{\sigma}(s))(u_{\gamma}(t_{1}-s,x_{\sigma}(s))-\bar{\lambda}_{\gamma}(s))_{+}\Big{)}\,ds.$

We now assume that $x_{\sigma}$ has the form

$x_{\sigma}(s)=\sigma\xi(s)+x_{0}+z(s)(x_{1}-x_{0}),$

where $\xi:[0,t_{1}-t_{0}]\to[0,1]$ satisfies $\xi(0)=\xi(t_{1}-t_{0})=0$ and
$z:[0,t_{1}-t_{0}]\to[0,1]$ is an increasing function such that $z(0)=0$ and
$z(t_{1}-t_{0})=1$. For notational simplicity, we will write
$\alpha_{\gamma}=\varphi_{\gamma}(u_{\gamma}-\bar{\lambda}_{\gamma})_{+}$.
Averaging (3.15) over $B_{1}$ we see that

$\varphi_{\gamma}(x_{0},t_{1}-t_{0})e^{-\Lambda_{\gamma}(t_{1}-t_{0})}\leq\varphi_{\gamma}(x_{1},0)+\frac{1}{4|B_{1}|}\int_{0}^{t_{1}-t_{0}}\int_{B_{1}}e^{-\Lambda_{\gamma}(s)}\Big{(}|\sigma|^{2}|\xi^{\prime}(s)|^{2}+|x_{1}-x_{0}|^{2}z^{\prime}(s)^{2}+\alpha_{\gamma}(s,x_{\sigma}(s))\Big{)}\,d\sigma\,ds.$

The optimality condition for $z$ implies that
$(z^{\prime}(s)e^{-\Lambda_{\gamma}(s)})^{\prime}=0$, therefore
$z^{\prime}(s)=\frac{e^{\Lambda_{\gamma}(s)}}{\int_{0}^{t_{1}-t_{0}}e^{\Lambda_{\gamma}(s)}\,ds}$.
Thus, making this choice we see that

$\varphi_{\gamma}(x_{0},t_{1}-t_{0})e^{-\Lambda_{\gamma}(t_{1}-t_{0})}\leq\varphi_{\gamma}(x_{1},0)+\frac{|x_{1}-x_{0}|^{2}}{4\int_{0}^{t_{1}-t_{0}}e^{\Lambda_{\gamma}(s)}\,ds}+\frac{1}{4|B_{1}|}\int_{0}^{t_{1}-t_{0}}\int_{B_{1}}e^{-\Lambda_{\gamma}(s)}\Big{(}|\sigma|^{2}|\xi^{\prime}(s)|^{2}+\alpha_{\gamma}(s,x_{\sigma}(s))\Big{)}\,d\sigma\,ds.$

Changing variables $y=x_{\sigma}$, it follows that

$\frac{1}{|B_{1}|}\int_{B_{1}}\alpha_{\gamma}(s,x_{\sigma}(s))d\sigma=\frac{\xi(s)^{-d}}{|B_{1}|}\int_{B_{\xi(s)}(x_{0}+z(s)(x_{1}-x_{0}))}\alpha_{\gamma}(s,y)dy$

where $B_{\xi(s)}(x_{0}+z(s)(x_{1}-x_{0}))$ is the ball of radius $\xi(s)$
centered at $x_{0}+z(s)(x_{1}-x_{0})$.

Using Hölder’s inequality with exponent $2d$, it follows that the above
quantity is bounded above by
$\xi(s)^{-1/2}\lVert\alpha_{\gamma}\rVert_{L^{2d}(\\{s\\}\times\mathbb{R}^{d})}.$
Hence, after dropping the good term $e^{-\Lambda_{\gamma}(s)}$ in the last
integral we see that

$\varphi_{\gamma}(x_{0},t_{1}-t_{0})e^{-\Lambda_{\gamma}(t_{1}-t_{0})}\leq\varphi_{\gamma}(x_{1},0)+\frac{|x_{1}-x_{0}|^{2}}{4\int_{0}^{t_{1}-t_{0}}e^{\Lambda_{\gamma}(s)}\,ds}+\frac{1}{4}\int_{0}^{t_{1}-t_{0}}\big{(}|\xi^{\prime}(s)|^{2}+\xi^{-1/2}(s)\lVert\alpha_{\gamma}(s,\cdot)\rVert_{L^{2d}(\mathbb{R}^{d})}\big{)}\,ds.$

Fix some $a>0$ and set

$\xi(s):=\begin{cases}as^{3/4}&\textup{if}\;\;t_{0}\leq s<(t_{1}-t_{0})/2,\\\
a(t_{1}-s)^{3/4}&\textup{if}\;\;(t_{1}-t_{0})/2\leq s\leq t_{1}.\end{cases}$

Using Hölder’s inequality with exponent $2$ on
$\int_{0}^{t}\xi^{-1/2}(s)\lVert\alpha_{\gamma}(s,\cdot)\rVert_{L^{2d}(\mathbb{R}^{d})}\,ds$,
we see that

$\lVert\xi^{\prime}\rVert_{L^{2}([t_{0},t_{1}])}^{2}\leq
C(t_{1}-t_{0})^{1/2}a^{2},\quad\lVert\xi^{-1/2}\rVert_{L^{2}([0,t_{1}-t_{0}])}\leq
Ca^{-1/2}(t_{1}-t_{0})^{1/8}$

thus,

$\varphi_{\gamma}(x_{0},t_{1}-t_{0})e^{-\Lambda_{\gamma}(t_{1}-t_{0})}\leq\varphi_{\gamma}(x_{1},0)+\frac{|x_{1}-x_{0}|^{2}}{4\int_{0}^{t_{1}-t_{0}}e^{\Lambda_{\gamma}(s)}\,ds}+C((t_{1}-t_{0})^{1/8}a^{-1/2}\lVert\alpha_{\gamma}\rVert_{L^{2}([0,t_{1}-t_{0}];L^{2d}(\mathbb{R}^{d}))}+(t_{1}-t_{0})^{1/2}a^{2})$

Optimizing over $a>0$, we obtain

(3.16)
$\varphi_{\gamma}(x_{0},t_{1}-t_{0})e^{-\Lambda_{\gamma}(t_{1}-t_{0})}\leq\varphi_{\gamma}(x_{1},0)+\frac{|x_{1}-x_{0}|^{2}}{4\int_{0}^{t_{1}-t_{0}}e^{\Lambda_{\gamma}(s)}\,ds}+C(t_{1}-t_{0})^{\frac{3}{10}}\lVert\alpha_{\gamma}\rVert_{L^{2}([0,t_{1}-t_{0}];L^{2d}(\mathbb{R}^{d}))}^{4/5}$

for a potentially different constant $C>0$.

Finally, it remains to estimate
$\lVert\alpha_{\gamma}\rVert_{L^{2}([0,t];L^{2d}(\mathbb{R}^{d}))}$. Recalling
that $\alpha_{\gamma}=\varphi_{\gamma}(u_{\gamma}-\lambda)_{+}$, we may write

$\lVert\alpha_{\gamma}\rVert_{L^{2d}(\\{s\\}\times\mathbb{R}^{d}))}^{2d}\leq\lVert\varphi_{\gamma}\rVert_{L^{\infty}([t_{0},t_{1}]\times\mathbb{R}^{d})}^{2}\int_{0}^{\infty}2dv^{2d-1}|\\{x\in\mathbb{R}^{d}:u_{\gamma,+}(t_{1}-s,x)>v+\bar{\lambda}_{\gamma}(s)\\}|\,dv$

By Chebyshev’s inequality, for any strictly increasing function
$f:\mathbb{R}\to\mathbb{R}$,

$\leq\lVert
p_{\gamma}\rVert_{L^{\infty}([t_{0},t_{1}]\times\mathbb{R}^{d})}^{2}\int_{0}^{\infty}2dv^{2d-1}\frac{\lVert
f(u_{\gamma,+})-f(0)\rVert_{L^{1}(\\{t_{1}-s\\}\times\mathbb{R}^{d})}}{f(\bar{\lambda}_{\gamma}(s)+v)-f(0)}\,dv$

If we choose $f(a)=\exp(ba)-1$, then we see that

$\lVert\alpha_{\gamma}\rVert_{L^{2d}(\\{s\\}\times\mathbb{R}^{d}))}\leq
Ce^{-b\bar{\lambda}_{\gamma}(s)}\lVert\exp(bu_{\gamma,+})-1\rVert_{L^{1}(\\{t_{1}-s\\}\times\mathbb{R}^{d})}=Ce^{-b\lambda(s)},$

where one should note carefully that now $\bar{\lambda}_{\gamma}$ has been
replaced by $\lambda$ in the last right-hand term. Thus, we have the estimate

$\lVert\alpha_{\gamma}\rVert_{L^{2}([t_{0},t_{1}];L^{2d}(\mathbb{R}^{d}))}^{2}\lesssim\int_{0}^{t_{1}-t_{1}}e^{-\frac{5}{2}\lambda(s)}\leq(t_{1}-t_{0})e^{-\frac{5}{2}\lambda(t_{1}-t_{0})},$

hence,

$\lVert\alpha_{\gamma}\rVert_{L^{2}([t_{0},t_{1}];L^{2d}(\mathbb{R}^{d}))}^{4/5}\lesssim(t_{1}-t_{0})^{2/5}e^{-\lambda(t_{1}-t_{0})}.$

Combining our work, we now have

(3.17)
$\varphi_{\gamma}(x_{0},t_{1}-t_{0})e^{-\Lambda_{\gamma}(t_{1}-t_{0})}\leq\varphi_{\gamma}(x_{1},0)+\frac{|x_{1}-x_{0}|^{2}}{4\int_{0}^{t_{1}-t_{0}}e^{\Lambda_{\gamma}(s)}\,ds}+C(t_{1}-t_{0})^{7/10}e^{-\lambda(t_{1}-t_{0})},$

for some possibly new constant $C$. The result follows after replacing
$\varphi_{\gamma}$ with $p_{\gamma}$ and multiplying both sides by
$e^{\Lambda_{\gamma}(t_{1}-t_{0})}$

∎

Before we can show that the Hopf-Lax formula also holds for the limiting
pressure, we first need a Lemma that gives us a pointwise well-defined
representative of our weak solution $p$. The argument is a simple adaptation
of a result from [MPQ17].

######  Lemma 3.6.

Suppose that $(\rho,p,n)$ is a weak solution to (1.1-1.2). $p$ can be
redefined on a set of measure zero so that

(3.18) $p(x,t)=\lim_{r\to
0}\frac{1}{r^{2}|B_{r}|}\int_{B_{r}(x)}\int_{t}^{t+r^{2}}p(y,s)\,ds\,dy$

for all $(x,t)$. With this definition, $p$ is spacetime upper semicontinuous
and for all $x\in\mathbb{R}^{d}$ the mapping $t\mapsto p(x,t)$ is continuous
from the right.

###### Proof.

From our control on $u_{\gamma}$ and the relation $\Delta
p_{\gamma}=-n_{\gamma}-\frac{1}{\gamma}u_{\gamma}$ it follows that after
taking limits, we have

$\Delta p\geq-n\geq-n_{0}$

in the sense of spacetime distributions. Therefore, for any $\epsilon>0$,

$\Delta\Big{(}\frac{1}{\epsilon^{2}}\int_{0}^{\epsilon^{2}}p(\cdot,t+s)\,ds\Big{)}\geq-
n_{0}$

in the sense of space distributions. Hence, the mean value property for
Laplace’s equation implies that for all $x\in\mathbb{R}^{d}$ and $t>0$ the
function

$\phi(r,\epsilon):=\frac{1}{\epsilon^{2}|B_{r}|}\int_{B_{r}(x)}\int_{0}^{\epsilon^{2}}p(y,t+s)+\frac{n_{0}}{2d}|y-x|^{2}\,ds\,dy$

is non-decreasing with respect to $r$. We also note that for $0\leq
r^{\prime}\leq r$ we have

$\phi(r,r)-\phi(r,r^{\prime})=\frac{1}{r^{2}|B_{r}|}\int_{B_{r}(x)}\int_{0}^{r^{2}}p(y,t+s)-p(y,t+(\frac{r^{\prime}}{r})^{2}s)\,ds\,dy$
$\geq\frac{1}{r^{2}|B_{r}|}\int_{B_{r}(x)}\int_{0}^{r^{2}}\int_{(\frac{r^{\prime}}{r})^{2}s}^{s}[\partial_{t}p(y,a)]_{-}\,da\,ds\,dy.$

From Corollary 3.4, it follows that $[\partial_{t}p]_{-}$ is bounded in
$L^{q}(\mathbb{R}^{d}\times[0,\tau])$ for any $q\in[1,\infty)$. Therefore,

$\phi(r,r)-\phi(r,r^{\prime})\geq-|B_{r}|^{-1/q}(r^{2}-r^{\prime
2})^{1-1/q}\lVert[\partial_{t}p]_{-}\rVert_{L^{q}([0,\tau]\times\mathbb{R}^{d})}\geq-
Cr^{1-(d+1)/q}(r-r^{\prime})^{1-1/q}$

for some constant $C>0$.

By choosing $q>d+1$, we can conclude that there exists a Holder continuous
function $g$ such that $r\mapsto\phi(r,r)+g(r)$ is nondecreasing and $g(0)=0$.
As a result, $\lim_{r\to 0^{+}}\phi(r,r)$ must exist for all $(x,t)$. Hence,
(3.18) is well defined everywhere. The Lebesgue differentiation theorem also
implies that our redefinition only changes $p$ on a set of measure zero.

Finally, to see that $p$ is upper semicontinuous, we note that $\lim_{r\to
0}\phi(r,r)=\lim_{r\to 0}\phi(r,r)+g(r)=\inf_{r>0}\phi(r,r)+g(r)$. Thus, we
may write

$p(x,t)=\inf_{r>0}g(r)+\frac{1}{r^{2}|B_{r}|}\int_{B_{r}(x)}\int_{t}^{t+r^{2}}p(y,s)\,ds\,dy.$

The infimum over a family of functions always produces an upper semicontinuous
function, hence, $p$ is upper semicontinuous.

∎

At last we obtain the main result of this Section, the Hopf-Lax formula for
our limit pressure $p$.

######  Proposition 3.7.

Given any points $(x_{1},t_{1}),(x_{0},t_{0})$ with $t_{0}<t_{1}$ and a
decreasing function $\lambda\in L^{1}([0,t_{1}-t_{0}])$, there exists a
constant $C=C(t_{1},d)$ such that

(3.19) $p(x_{0},t_{0})\leq
e^{\Lambda(t_{1}-t_{0})}\Big{(}p(x_{1},t_{1})+\frac{|x_{1}-x_{0}|^{2}}{4\int_{0}^{t_{1}-t_{0}}e^{\Lambda(s)}\,ds}+C(t_{1}-t_{0})^{7/10}e^{-\lambda(t_{1}-t_{0})}\Big{)}$

where $b$ is the constant from Proposition 3.3 and

$\Lambda(t):=\frac{5}{4b}\int_{0}^{t}\lambda(s)\,ds+\frac{t}{b}\log(1+\frac{C}{t})$

###### Proof.

Using the formula from Lemma 3.6, we have

$p(x_{0},t_{0})=\lim_{r\to
0^{+}}\frac{1}{r^{2}|B_{r}|}\int_{B_{r}(x)}\int_{t}^{t+r^{2}}p(y,s)\,ds\,dy.$

Choose a point $x_{2}\in\mathbb{R}^{d}$ and $t_{2}>t_{0}$ such that
$p_{\gamma}(x_{2},t_{2})$ converges to $p(x_{2},t_{2})$ along some subsequence
$\gamma_{k}$. Using the $L^{2}_{t}H^{1}_{x}$ strong convergence of
$p_{\gamma}$ to $p$ and then applying Lemma 3.5, we have for any
$x_{2}\in\mathbb{R}^{d}$ and $t_{2}>t_{0}$

$p(x_{0},t_{0})=\lim_{r\to
0^{+}}\lim_{k\to\infty}\frac{1}{r^{2}|B_{r}|}\int_{B_{r}(x)}\int_{t_{0}}^{t_{0}+r^{2}}p_{\gamma_{k}}(y,s)\,ds\,dy\leq$
$\lim_{r\to
0^{+}}\lim_{k\to\infty}\frac{1}{r^{2}|B_{r}|}\int_{B_{r}(x)}\int_{t_{0}}^{t_{0}+r^{2}}e^{\Lambda_{\gamma_{k}}(t_{2}-s)}\Big{(}p_{\gamma_{k}}(x_{2},t_{2})+\frac{|x_{2}-x_{0}|^{2}}{4\int_{0}^{t_{2}-s}e^{\Lambda_{\gamma_{k}}(a)}\,da}+C(t_{2}-t_{0})^{7/10}e^{-\lambda(t_{2}-s)}\Big{)}\,dy\,ds.$

Recall that

$\Lambda_{\gamma}(t):=\frac{5}{4b}\int_{0}^{t}\lambda(a)\,da+\frac{1}{b}\int_{0}^{t}\log(1+\lVert\exp(bu_{\gamma,+})-1\rVert_{L^{1}(\\{t_{1}-a\\}\times\mathbb{R}^{d})})\,da.$

Applying Jensen’s inequality, we have the bound

$\Lambda_{\gamma}(t)\leq\frac{5}{4b}\int_{0}^{t}\lambda(a)\,da+\frac{t}{b}\log(1+\frac{1}{t}\lVert\exp(bu_{\gamma,+})-1\rVert_{L^{1}([t_{1}-t,t_{1}]\times\mathbb{R}^{d})}).$

Hence, we can find a potentially new constant $C=C(\tau,d)>0$ such that

$\Lambda_{\gamma}(t)\leq\frac{5}{4b}\int_{0}^{t}\lambda(a)\,da+\frac{t}{b}\log(1+\frac{C}{t})=\Lambda(t).$

for all $\gamma$. Therefore,

$p(x_{0},t_{0})\leq
e^{\Lambda(t_{2}-t_{0})}\Big{(}p(x_{2},t_{2})+\frac{|x_{2}-x_{0}|^{2}}{4\int_{0}^{t_{2}-t_{0}}e^{\Lambda(s)}\,ds}+C(t_{2}-t_{0})^{7/10}e^{-\lambda(t_{2}-t_{0})}\Big{)}.$

Since $p_{\gamma}$ converges pointwise almost everywhere to $p$ along
appropriate subsequences, it follows that

$p(x_{0},t_{0})\leq
e^{\Lambda(t_{2}-t_{0})}\Big{(}p(x_{2},t_{2})+\frac{|x_{2}-x_{0}|^{2}}{4\int_{0}^{t_{2}-t_{0}}e^{\Lambda(a)}\,da}+C(t_{2}-t_{0})^{7/10}e^{-\lambda(t_{2}-t_{0})}\Big{)}.$

for a dense set of $(x_{2},t_{2})$ with $t_{2}>t_{0}$. The result now follows
from the upper semicontinuity of $p$. ∎

##  4\. Hölder continuity of the Hitting time

We are now going to construct a radial supersolution which will give an upper
bound for the rate of expansion for the tumor, and thus, will yield a lower
bound on the arrival times. Given a point of interest
$x_{0}\in\mathbb{R}^{d}$, our supersolution will be defined on the time-
dependent annulus

$A(t):=\\{t\\}\times\\{x:r(t)\leq|x-x_{0}|\leq mr(t)\\},\quad
A=\bigcup_{t\in[0,\epsilon]}A(t),$

for some $m>1$ and a function $r(t)\geq 0$ that we will define shortly. Given
some starting time $t_{0}$, for each $t,r\geq 0$ we define

$\bar{p}(t,r):=\sup_{x\in B_{r}(x_{0})}p(t+t_{0},x),$

where the sup is well defined since $p$ is upper semicontinuous in space.

Now we will construct our supersolution $\psi(t,x)$ by solving

$\begin{cases}-\Delta\psi(t,x)=\bar{n}_{0}&\textup{if}\;r(t)<|x-x_{0}|<mr(t),\\\
\psi(t,x)=0&\textup{if}\;|x-x_{0}|\leq r(t),\\\
\psi(t,x)=\bar{p}(t,|x-x_{0}|)&\textup{if}\;|x-x_{0}|\geq mr(t).\\\
\end{cases}$

On $A(t)$, the equation admits the explicit radial solution

(4.1)
$\psi(t,x)=h(t)\Gamma_{d}(|x-x_{0}|)-\frac{\bar{n}(0)}{2d}|x-x_{0}|^{2}+g(t),$

where $\Gamma_{d}$ is the fundamental solution of the Laplace equation in
dimension $d$, i.e. $\Gamma_{d}^{\prime}(r)=r^{1-d}$,

(4.2)
$h(t):=\frac{\bar{p}(t,mr(t))+(2d)^{-1}\bar{n}(0)(m^{2}-1)r(t)^{2}}{\Gamma_{d}(mr(t))-\Gamma_{d}(r(t))},$

and

(4.3) $g(t):=\frac{\bar{n}_{0}}{2d}-h(t)\Gamma_{d}(r(t)).$

Finally, we define $r(t)$ by choosing some initial data $r(0)$ and then
solving the ODE

(4.4) $r^{\prime}(t)=-|\nabla\psi(t,y)|$

where the right hand side is evaluated at any point $y$ such that
$|y-x_{0}|=r(t)$.

We now show that $\psi$ is indeed a supersolution as long as comparison holds
at initial time. Due to the lack of regularity for the pressure variable, we
establish comparison using the time integrated versions of $\psi$ and $p$.
This creates an annoying issue where it is difficult to establish that the
boundary data stays ordered as the annulus moves. To avoid this problem, we
establish comparison by first going through a sequence of supersolutions
$\psi_{k}$, where the $\psi_{k}$ are defined on modified annuli whose outer
radii are taken to be piecewise constant in time.

######  Lemma 4.1.

Let $\mu(t,x)$ be the characteristic function of the set
$\\{x\in\mathbb{R}^{d}:|x-x_{0}|\geq r(t)\\}$. If $\mu(0,x)\leq\rho(t_{0},x)$
for almost every $x\in\mathbb{R}^{d}$, then $p(t_{0}+t,x)\leq\psi(t,x)$ for
almost every $x\in\mathbb{R}^{d}$ and almost every time $t\geq 0$

###### Proof.

As we noted above, we will first prove the comparison for a modified sequence
of supersolutions $\psi_{k}$. The $\psi_{k}$ will be defined in precisely the
same way as $\psi$, except that we will modify the construction of the moving
annulus. Hence, given radii $r_{k}(t)<R_{k}(t)$, we define $\psi_{k}$ by
solving

$\begin{cases}-\Delta\psi_{k}(t,x)=\bar{n}_{0}&\textup{if}\;r_{k}(t)<|x-x_{0}|<R_{k}(t),\\\
\psi_{k}(t,x)=0&\textup{if}\;|x-x_{0}|\leq r_{k}(t),\\\
\psi_{k}(t,x)=\bar{p}(t,|x-x_{0}|)&\textup{if}\;|x-x_{0}|\geq R_{k}(t).\\\
\end{cases}$

$r_{k}(t)$ will be defined as before via the ODE
$r_{k}^{\prime}(t)=-|\nabla\psi_{k}(t,y)|$ where $y$ is any point satisfying
$|y-x_{0}|=r_{k}(t)$. We then define $R_{k}$ by setting

$R_{k}(t):=mr_{k}(t_{k,j}),\quad\textup{if}\;t\in[t_{k,j},t_{k,j+1}),$

where we inductively define the points $t_{k,j}$ by setting $t_{k,0}=0$ and
then taking

$t_{k,j+1}:=\inf\\{t\geq t_{k,j}:r_{k}(t)<(1-\frac{1}{k+1})r_{k}(t_{k,j})\\}.$

As before, on the annulus $r_{k}(t)\leq|x-x_{0}|\leq R_{k}(t)$, the $\psi_{k}$
will admit the explicit radial solutions

(4.5)
$\psi_{k}(t,x)=h_{k}(t)\Gamma_{d}(|x-x_{0}|)-\frac{\bar{n}(0)}{2d}|x-x_{0}|^{2}+g_{k}(t),$

where

(4.6)
$h_{k}(t):=\frac{\bar{p}(t,R_{k}(t))+(2d)^{-1}\bar{n}(0)(R_{k}(t)^{2}-r_{k}(t)^{2})}{\Gamma_{d}(R_{k}(t))-\Gamma_{d}(r_{k}(t))},$

and

(4.7) $g_{k}(t):=\frac{\bar{n}_{0}}{2d}-h_{k}(t)\Gamma_{d}(r_{k}(t)).$

Let $\Psi_{k}(t,x)=\int_{0}^{t}\psi_{k}(s,x)\,ds$. Since $\psi_{k}$ is clearly
Lipschitz in space on $|x-x_{0}|\leq R_{k}(t)$, it follows that $\Psi_{k}$ is
Lipschitz in space on $|x-x_{0}|\leq R_{k}(t)$ and

$\nabla\Psi_{k}(t,x)=\int_{0}^{t}\nabla\psi_{k}(s,x)\,ds$

almost everywhere on $|x-x_{0}|\leq R_{k}(t)$. Define $\tilde{t}_{k}(r)$ to be
the inverse function of $r_{k}(t)$. From the definition of $\psi_{k}$, it
follows that

$\nabla\Psi_{k}(t,x)=\int_{\min(t,\tilde{t}_{k}(|x-x_{0}|))}^{t}\nabla\psi_{k}(s,x)\,ds.$

Now if $x$ is a point such that $|x-x_{0}|<R_{k}(t)$ and $|x-x_{0}|<r(0)$,
then for each fixed $s\in(\tilde{t}(|x-x_{0}|),t]$, there exists a
neighborhood of $x$ such that $\nabla\psi_{k}(s,x)$ is differentiable and
$-\Delta\psi_{k}(s,x)=\bar{n}_{0}$. Thus, it follows that

$-\Delta\Psi_{k}(t,x)=\operatorname{sgn}_{+}\big{(}t-\tilde{t}(|x-x_{0}|)\big{)}\tilde{t}^{\prime}(|x-x_{0}|)|\nabla\psi_{k}(\tilde{t}(|x-x_{0}|),x)|+\int_{\min(t,\tilde{t}_{k}(|x-x_{0}|))}^{t}\bar{n}_{0}\,ds.$

Since
$r^{\prime}_{k}(\tilde{t}_{k}(|x-x_{0}|))=-|\nabla\psi(\tilde{t}(|x-x_{0}|),x)|$
and $\tilde{t}_{k}(r)$ is the inverse of $r_{k}(t)$, we see that

$-\Delta\Psi_{k}(t,x)=-\operatorname{sgn}_{+}\big{(}t-\tilde{t}(|x-x_{0}|)\big{)}+(t-\tilde{t}(|x-x_{0}|)_{+}\bar{n}_{0}.$

On the other hand, if $x$ is a point such that $r(0)<|x-x_{0}|<R_{k}(t)$, then

$-\Delta\Psi_{k}(t,x)=t\bar{n}_{0}.$

Let $\mu_{k}(t,x)$ be the characteristic function of the set
$\\{(t,x):|x-x_{0}|\geq r_{k}(t)\\}$ and note that
$\mu_{k}(t,x)=\operatorname{sgn}_{+}(t-\tilde{t}_{k}(|x-x_{0}|))$ and
$\int_{0}^{t}\mu_{k}(s,x)\,ds=(t-\tilde{t}_{k}(|x-x_{0}|)_{+}.$ Combining our
work from above, we can conclude that for almost every $x$ satisfying
$|x-x_{0}|<R_{k}(t)$ we have

$-\Delta\Psi_{k}(t,x)=\mu_{k}(0,x)-\mu_{k}(t,x)+\int_{0}^{t}\mu_{k}(s,x)\bar{n}_{0},$

and for almost all $x\in\mathbb{R}^{d}$ we have
$\Psi_{k}(t,x)(1-\mu_{k}(t,x))=0,$ as well as $\psi_{k}(1-\mu_{k}(t,x))=0$.

Now let us define the time shifted variables
$\tilde{w}(t,x):=\int_{0}^{t}p(s+t_{0},x)\,ds$,
$\tilde{\rho}(t,x):=\rho(t_{0}+t,x)$, and $\tilde{n}(t,x):=n(t_{0}+t,x)$. It
then follows that
$-\Delta\tilde{w}(t,x)=\tilde{\rho}(t,x)-\tilde{\rho}(0,x)+\int_{t_{0}}^{t}\tilde{\rho}(s,x)\tilde{n}(s,x)$
and $(1-\rho(t+t_{0},x))w_{0}=0$ almost everywhere. For any time
$t\in[0,t_{k,1})$, the definition of $\psi_{k}$ guarantees that
$\Psi_{k}(t,x)\geq w(t,x)$ for all $x$ satisfying $|x-x_{0}|=R_{k}(t)=mr(0)$.
Hence, for any $t\in[0,t_{1})$ and any increasing $C^{1}$ function
$\eta:\mathbb{R}\to\mathbb{R}$ such that $\eta(a)=0$ if $a\leq 0$, we have

$\int_{\\{|x-x_{0}|\leq
R_{k}(0)\\}}(\tilde{\rho}-\mu_{k})\eta(\tilde{w}-\Psi_{k})+\eta^{\prime}(\tilde{w}-\Psi_{k})|\nabla(\tilde{w}-\Psi_{k})|^{2}\leq\int_{\\{|x-x_{0}|\leq
R_{k}(0)\\}}\int_{0}^{t}\eta(\tilde{w}-\Psi_{k})(\tilde{\rho}\tilde{n}-\mu_{k}\bar{n}_{0})$

Letting $\eta$ approach $\operatorname{sgn}_{+}$ and using the fact that
$\operatorname{sgn}_{+}(\tilde{w}-\Psi_{k})=\operatorname{sgn}_{+}(\tilde{\rho}-\mu_{k})$,
we can conclude that

$\int_{\\{|x-x_{0}|\leq
R_{k}(0)\\}}(\tilde{\rho}-\mu_{k})_{+}\leq\int_{\\{|x-x_{0}|\leq
R_{k}(0)\\}}\bar{n}_{0}\int_{0}^{t}(\tilde{\rho}-\mu_{k})_{+}.$

Hence, Gronwall’s inequality now implies that
$\tilde{\rho}(t,x)\leq\mu_{k}(t,x)$ for all $t\in[0,t_{k,1})$ and almost all
$x\in\mathbb{R}^{d}$ (recall it is immediate that $\tilde{\rho}\leq\mu_{k}$ on
$|x-x_{0}|\geq R(0)$ from the definition of $\mu_{k}$). The masses of the
differences $\mu_{k}(t,x)-\mu_{k}(0,x)$ and
$\tilde{\rho}(t,x)-\tilde{\rho}(0,x)$ are continuous functions of time,
therefore, the ordering $\tilde{\rho}\leq\mu_{k}$ must hold at time $t_{1}$.
This allows us to run the above argument on $[t_{k,1},t_{k,2})$. Iterating, we
conclude that the ordering $\tilde{\rho}\leq\mu_{k}$ must hold for all times
$t$ when $r_{k}(t)>0$.

Now we wish to argue that $\liminf_{k\to\infty}r_{k}(t)\geq r(t)$. Let

$t_{*}=\inf\\{t>0:\liminf_{k\to\infty}r_{k}(t)<r(t)\\},$

and note that $\liminf_{k\to\infty}r_{k}(t_{*})=r(t_{*})$. Using the explicit
formulas (4.1) and (4.5), as well as the upper semicontinuity of
$r\mapsto\bar{p}(t,r)$, it follows that

$\liminf_{k\to\infty}r_{k}^{\prime}(t_{*})\geq r^{\prime}(t_{*})$

whenever $r(t_{*})>0$. Hence, $r(t)\leq\liminf_{k\to\infty}r_{k}(t)$ for all
times where $r(t)>0$. This implies that $\tilde{\rho}(t,x)\leq\mu(t,x)$ for
all $t$ and almost all $x$.

Finally, we note that the ordering $\tilde{\rho}(t,x)\leq\mu(t,x)$ implies
that for almost every time $t$

$(p-\psi)_{+}(\Delta p+n)=0,\quad(p-\psi)_{+}(\Delta\psi+\bar{n}_{0})=0.$

distributionally. Thus, for any $T>0$

$\int_{Q_{T}}|\nabla(p-\psi)_{+}|^{2}=\int_{Q_{T}}(p-\psi)_{+}(n-\bar{n}_{0})$

which is only possible if $(p-\psi)_{+}=0$ almost everywhere.

∎

We can now use this barrier supersolution to get bounds on the Hölder
continuity of the hitting time. The key is to use our Hopf-Lax estimate from
3.7 to ensure that the supersolution arrives at the point of interest at the
correct time.

######  Theorem 4.2.

$T$ is locally Hölder continuous on the set
$\\{x\in\mathbb{R}^{d}:0<T(x)<\infty\\}$. In particular, for any
$x_{1}\in\mathbb{R}^{d}$ such that $T(x_{1})\in(0,\infty)$, we have

(4.8) $\sup_{y\in B_{R}(x_{1})}T(x_{1})-T(y)\lesssim R^{\alpha_{d}}$

for all $R>0$ sufficiently small, where

(4.9) $\alpha_{d}:=\begin{cases}\frac{2}{e}&\textup{if}\;d=2,\\\
2(\frac{2}{d})^{\frac{d}{d-2}}&\textup{if}\;d>2.\\\ \end{cases}$

###### Proof.

Let $\epsilon>0$ be a small value that we will choose later. Let

$\delta=\delta(\epsilon):=\inf\\{R>0:\sup_{y\in
B_{R}(x_{1})}T(x_{1})-T(y)\geq\epsilon\\}.$

Since $T$ is continuous at $x_{1}$, it follows that $\lim_{\epsilon\to
0}\delta(\epsilon)=0$. Let $t_{0}=T(x_{1})-\epsilon$ and $t_{1}=T(x_{1})$.
Thanks to the super solution that we have constructed above, we know that

$\inf_{y\in B_{r(t)}(x_{1})}T(y)\geq t_{0}+t,$

which implies

$\sup_{y\in B_{r(t)}(x_{1})}T(x_{1})-T(y)\leq(t_{1}-t_{0}-t).$

Hence, if we can provide lower bounds on $r(t)$ in terms of $t$, we can get a
Hölder estimate for $T$ at $x_{1}$. In particular, a bound of the form
$(t_{1}-t_{0}-t)^{1/\alpha}\lesssim r(t)$ will imply that $\sup_{y\in
B_{R}(x_{1})}T(x_{1})-T(y)\lesssim R^{\alpha}$.

To bound $r(t)$ from below, we must consider the ODE (4.4), which can be
simplified to

$r^{\prime}(t)=-|h(t)||\Gamma_{d}^{\prime}(r(t))|-\frac{\bar{n}(0)}{d}r(t).$

Noting that in any dimension there exists a function $\xi_{d}(m)$ such that
$\frac{|\Gamma_{d}^{\prime}(r(t))|}{|\Gamma_{d}(m(r(t))-\Gamma_{d}(r(t))|}=r(t)^{-1}\xi_{d}(m)$,
it follows from the structure of $h$ and the ODE that there exists some
constant $K>0$ such that

(4.10) $r^{\prime}(t)+Kr(t)\geq-\frac{\bar{p}(t,mr(t))\xi_{d}(m)}{r(t)}.$

Now we want to estimate $\bar{p}(t,mr(t))$. To do so, we will apply the bounds
from Lemma 3.5, choosing to evaluate $p$ at $(x_{1},t_{1})$ and leaving the
choice of $\lambda\in L^{1}([0,t_{1}-t_{0}])$ until later. With these choices,
we see that

$\begin{array}[]{lll}\bar{p}(t,mr(t))&=&\sup_{x\in
B_{mr(t))}(x_{1})}p(t+t_{0},x)\\\ &\leq&\sup_{x\in
B_{mr(t))}(x_{1})}H(t)|x-x_{1}|^{2}+F(t)=m^{2}r^{2}H(t)+F(t),\end{array}$

where we have defined

(4.11)
$H(t):=e^{\Lambda(t_{1}-t_{0}-t)}(4\int_{0}^{t_{1}-t_{0}-t}e^{\Lambda(s)}\,ds)^{-1},\quad
F(t):=C(t_{1}-t_{0}-t)^{7/10}e^{-\lambda(t_{1}-t_{0}-t)+\Lambda(t_{1}-t_{0}-t)}$

for notational convenience.

Returning to equation (4.10) and applying the upper bound on $\bar{p}$
obtained above, we have

$r^{\prime}(t)+r(t)(K+m^{2}\xi_{d}(m)H(t))\geq-\xi_{d}(m)\frac{F(t)}{r(t)}$

Multiplying both sides by $2r(t)$ and defining $z(t)=r(t)^{2}$, we get

(4.12) $z^{\prime}(t)+z(t)(2K+2m^{2}\xi_{d}(m)H(t))\geq-\xi_{d}(m)F(t).$

Now we choose $m$ by optimizing $m^{2}\xi_{d}(m)$. Define

$\xi_{d}:=\inf_{m>1}\frac{m^{2}}{2}\xi_{d}(m).$

One can then check that $\xi_{d}=(\frac{d}{2})^{\frac{d}{d-2}}$ and
$\operatorname*{argmin}\frac{m^{2}}{2}\xi_{d}(m)=(\frac{d}{2})^{\frac{1}{d-2}}$
(where these should be understood in a limiting sense when $d=2$). Thus, we
have

(4.13) $z^{\prime}(t)+z(t)(2K+4\xi_{d}H(t))\geq-dF(t).$

Let $\bar{H}(t)=\int_{0}^{t}4H(s)\,ds$. Multiplying both sides of (4.13) by
$\exp(2Kt+\xi_{d}\bar{H}(t))$ and integrating in time, we can conclude that

(4.14) $z(t)e^{2Kt+\xi_{d}\bar{H}(t)}\geq
z(0)-d\int_{0}^{t}F(s)e^{2Ks+\xi_{d}\bar{H}(s)}\,ds.$

Now we need to provide upper bounds on $\exp(\xi_{d}\bar{H}(t))$. To do so, we
will need to make a choice for $\lambda$. Fix some $\theta>0$ and set

$\lambda(s)=\theta+s^{-1/2}.$

We then have

$\Lambda(t)=\frac{5}{4b}(\theta t+2t^{1/2})+\frac{t}{b}\log(1+\frac{C}{t}).$

Using the above estimates, we see that

$4H(t)\leq\frac{\exp\Big{(}\frac{5}{4b}(\theta(t_{1}-t_{0}-t)+2(t_{1}-t_{0}-t)^{1/2})+(t_{1}-t_{0}-t)\log(1+C/(t_{1}-t_{0}-t))\Big{)}}{\int_{0}^{t_{1}-t_{0}-t}e^{\frac{5}{4b}\theta
s}\,ds}=$
$\frac{\frac{5\theta}{4b}\exp\Big{(}2(t_{1}-t_{0}-t)^{1/2})+(t_{1}-t_{0}-t)\log(1+C/(t_{1}-t_{0}-t))\Big{)}}{1-e^{-\frac{5}{4b}\theta(t_{1}-t_{0}-t)}}.$

Since $(t_{1}-t_{0}-t)\leq(t_{1}-t_{0})=\epsilon$, we can assume that
$\epsilon$ is sufficiently small that

$4H(t)\leq\frac{5\theta}{4b(1-e^{-\frac{5}{4b}\theta(t_{1}-t_{0}-t)})}+\frac{20\theta(t_{1}-t_{0}-t)^{1/2}}{4b(1-e^{-\frac{5}{4b}\theta(t_{1}-t_{0}-t)})}$

Hence, for some possibly new constant $C>0$ independent of $\epsilon$ and
$\theta$ we get

$\bar{H}(t)\leq\log\big{(}\frac{e^{\frac{5\theta}{4b}\epsilon}-1}{e^{\frac{5\theta}{4b}(t_{1}-t_{0}-t)}-1}\big{)}+C(1+\epsilon^{3/2}\theta).$

Thus,

$\exp(\xi_{d}\bar{H}(t))\lesssim\big{(}\frac{e^{\frac{5\theta}{4b}\epsilon}-1}{e^{\frac{5\theta}{4b}(t_{1}-t_{0}-t)}-1}\big{)}^{\xi_{d}}\exp(C\epsilon^{3/2}\theta)$

Now we move to estimating $F(s)\exp(2Ks+\xi_{d}\bar{H}(s))$. From our choice
of $\lambda$, it is clear that $\lambda(t)\geq 2\Lambda(t)$ for all $t$
sufficiently small. Thus, it follows that

$F(s)\lesssim(t_{1}-t_{0}-s)^{7/10}e^{-\lambda(t_{1}-t_{0}-s)/2},$

hence, we have the bound

$F(s)\exp(2Ks+\xi_{d}\bar{H}(s))\lesssim\big{(}\frac{e^{\frac{5\theta}{4b}\epsilon}-1}{e^{\frac{5\theta}{4b}(t_{1}-t_{0}-s)}-1}\big{)}^{\xi_{d}}(t_{1}-t_{0}-s)^{7/10}\exp\Big{(}\frac{1}{2}\big{(}4Ks+(2C\epsilon^{3/2}-1)\theta-(t_{1}-t_{0}-s)^{-1/2})\big{)}\Big{)}$

Therefore, once $\epsilon$ is small enough that
$2C\epsilon^{3/2}+\frac{5}{4b}\epsilon<1$ we can choose $\theta$ large enough
that

$z(0)-d\int_{0}^{\epsilon}F(s)e^{2Ks+\xi_{d}\bar{H}(s)}\,ds\geq z(0)/2,$

and from there we can conclude that

$z(0)\lesssim z(t)e^{2Kt+\xi_{d}\bar{H}(t)}$

for all $t\in[0,\epsilon]$. This implies that

$z(0)(\theta(t_{1}-t_{0}-t))^{\xi_{d}}\lesssim z(t)=r(t)^{2}.$

Hence, $r(0)(t_{1}-t_{0}-t)^{\xi_{d}/2}\lesssim r(t).$ The result now follows
from the fact that $\alpha_{d}=2/\xi_{d}$.

∎

##  5\. Results from Obstacle Problem Theory

In this section, we use techniques from the theory of the obstacle problem to
study the local behavior of the interface. The main technique here is the
quadratic blowup, which classifies free boundary points into regular points,
where the zero set is asymptotically a half-space, and singular points, where
the zero set is asymptotically lower dimensional. With sufficiently regular
source term, the blowup limit approximates the solution at a uniform scale,
and we can use this to extract information on the local geometry of the
positive set. The regularity of the source term in the equation satisfied by
$w$ is governed by the regularity of the nutrient and the regularity of the
hitting time. Since the nutrient enjoys parabolic regularity as in Lemma 2.2,
Hölder continuity of the hitting time leads to Hölder continuous source, which
is enough to control the blowup limit at both types of free boundary points.

Thus, using that $T\in
C^{0,\alpha}_{\operatorname{\textup{loc}}}(\mathcal{O})$ for the
$\alpha\in(0,1)$ from Theorem 4.2, we show that the regular points form an
open set of full measure in the spacetime interface, on which $T$ improves to
locally Lipschitz and the spatial interface evolves as a locally $C^{1,1-}$
graph. The scale and bounds for which this regularity is achieved can be
quantified in terms of the Hölder seminorm of $T$ and the scale at which the
zero set achieves sufficiently large density near the regular point. We also
show Hölder regularity of the unit normal to the interface in spacetime, using
the spatial regularity of the interface and the monotonicity of its expansion.
Under the stronger assumption that no singular points occur for some time
interval, this lets us improve $T$ to
$C^{1,1/2-}_{\operatorname{\textup{loc}}}$ on the corresponding region.

As for singular points, we show that they form a relatively closed set in
$\mathcal{O}$ contained in a $C^{1}$ manifold of dimension $d-1$. This
improves on the standard obstacle problem result that the singular points at a
fixed time are contained in a $C^{1}$ manifold of dimension $d-1$ and implies
that the worse case, where the singular points have positive $d-1$ Hausdorff
measure for some time, only occurs for at most countably many times. Under the
additional assumption that $T$ is Lipschitz up to the singular set, we show a
stronger generic regularity result which gives that for a.e. time the singular
points have $d-2$ Hausdorff measure 0. In dimension 2, this would imply that
the times with singular points have zero measure, as a relatively closed
subset of $(0,\infty)$.

We note that we are not currently able to prove that $T$ is Lipschitz up to
singular points. It is not clear to what extent the obstacle problem can be
leveraged to understand the geometry of the patch at times just before a
singular point occurs, in order to prove nondegeneracy of the pressure.
Nondegeneracy at later times is also uncertain, but appears more tractable
since the blowup is available. For example, suppose one knew, for a singular
point $x_{0}$ and all sufficiently small $r$, that the set
$\\{w(\cdot,T(x_{0}))=0\\}\cap B_{r}(x_{0})$ is contained in a strip of width
$Cr^{1+\alpha}$, for some $C,\alpha$ depending on the source term. Assuming
such a strip condition, then the Hopf lemma could be applied to establish
nondegeneracy of $p$ near $x_{0}$ at times $t\geq T(x_{0})$. This strip
condition has been proven by [FS19] for singular points in the
$(d-1)$-dimensional stratum, albeit using methods which require much stronger
regularity than $C^{0,\alpha}$ source. A related result on the rate of
convergence of the quadratic blowup at singular points in dimension 2 has been
proven by [CSV18]. As far as we are aware, it is not currently known whether
or in what sense this strip condition may hold for the obstacle problem with
$C^{0,\alpha}$ source.

Nevertheless, we expect that $T$ is indeed Lipschitz, as it seems unlikely
that the pressure sometimes becomes degenerate, but only at instances of
merging or topological change. It is clear that we cannot hope for better than
Lipschitz, since $T$ cannot be differentiable when two pieces of the boundary
collide while traveling at different speeds. We also discuss examples in
Remark 5.7 which show that Lipschitz continuity of $T$ is sharp at regular
points without the additional assumption to give global control over singular
points.

For the obstacle problem, the $C^{1,\alpha}$ regularity of the free boundary
at regular points and the $C^{1}$ manifold covering singular points are well-
known for Hölder source (see appendix for more details). Thus, the main
challenge in lies in the analyzing the time-indexed family of obstacle
problems satisfied by $w(\cdot,t)$ to control these properties in time. We
note that such parameterized families have now been studied extensively for
the constant source obstacle problem with varying fixed boundary data, mainly
with the goal of understanding generic behavior of singular points ([Mon03],
[FROS20]). Our problem differs in that we must contend with a varying low
regularity source and no fixed boundary data, which rules out many of the
techniques typically used. In particular, the results of [FROS20], including
that the singular set has $(d-4)$-Hausdorff measure zero, do not appear to be
in reach with even Lipschitz source. Our approach draws from arguments in
[Mon03] to establish the $C^{1}$ manifold property for the singular set.
However, whereas comparison arguments with the fixed boundary data allow
Monneau to prove directly that the hitting time is Lipschitz, we must work
harder to get lower regularity for $T$. Finally, the analysis of the regular
set for the time-parameterized family, to our knowledge, is new. The main
facts we make use of are the Hölder continuity of $T$, the spacetime
continuity of $w$, the $L^{1}$ time-continuity of $\rho$, and the monotonicity
of $\rho$ and $w$ in time.

We remark that Proposition 2.12 shows that the interface strictly expands, and
in space-time is exactly the graph of $T$. Therefore, regularity improvements
to $T$ correspond exactly to regularity of the space-time interface as a
$d$-dimensional manifold. As a result, we will generally not consider the
space-time perspective directly, preferring to work with the subsets of
${\mathbb{R}}^{d}$ traced out by the moving interface.

Now, we proceed to study the local situation at the free boundary. As we noted
above, the new regularity of $T$ from Theorem 4.2 feeds back into the obstacle
problem satisfied by $w$ through the dependence of $\eta$, as defined in
(1.8), on $T$. We state this precisely below:

######  Lemma 5.1.

Up to $C^{1,1-}$, $\eta(\cdot,t)$ has the same spatial regularity as $T$ on
$\overline{\\{w(\cdot,t)>0\\}}$. In particular, for any $\tau>0$, we have
$\eta\in L^{\infty}_{t}C^{0,\alpha}_{x}([0,\tau];{\mathbb{R}}^{d})$.

###### Proof.

The first part follows immediately from the $L^{\infty}_{t}C^{1,1-}$
regularity of $n$, from Lemma 2.2. The second part follows from the
$C^{0,\alpha}$ regularity of $T$. ∎

The exact regularity of $\eta$ is relevant for determining the spatial
regularity of the free boundary near regular points. However, for most results
in this section, we only require Hölder continuity to give uniqueness of the
quadratic blowup limit introduced in Lemma 2.10, and to give spatial
equicontinuity of $\eta(\cdot,t)$. The uniqueness of the blowup limit and its
subsequent characterization is best expressed as the following dichotomy,
originally due to Caffarelli:

######  Lemma 5.2.

Let $u$ be a solution of the obstacle problem $\Delta u=f\chi_{\\{u>0\\}}$ in
${\mathbb{R}}^{d}$ with $f$ positive and $C^{0,\alpha}$ near 0. If
$0\in\partial\\{u>0\\}$, then one of the following holds:

  1. (1)

$\\{u=0\\}$ has density $\frac{1}{2}$ at 0, and the quadratic rescalings
$r^{-2}u(rx)$ converge in $C^{1,1-}(B_{1})$ to $\frac{f(0)}{2}(x\cdot
e)^{2}_{+}$ for some unit vector $e$.

  2. (2)

$\\{u=0\\}$ has density 0 at 0, and the quadratic rescalings $r^{-2}u(rx)$
converge in $C^{1,1-}(B_{1})$ to $\frac{f(0)}{2}x\cdot D^{2}u(0)x$, where
$D^{2}u(0)$ exists in the classical sense and is a positive semidefinite
matrix with trace 1.

Points of the first type are called regular points, and points of the second
type are called singular points.

This dichotomy was proven in [Caf98] for the constant source obstacle problem,
with the note that minor modifications could extend the proof to the Hölder
continuous case. An energetic criterion for the dichotomy appears in [Wei99].
Careful proofs for the uniqueness of the blowup limit in the Hölder continuous
case are given in [Bla00] for regular points and [Mon03] for singular points.

The dichotomy applies to the free boundary of $w(\cdot,t)$ for each $t$. We
let $R_{t}$ denote the regular points of $\partial\Omega_{t}$, and
$\Sigma_{t}$ denote the singular points of $\partial\Omega_{t}$, for the
obstacle problem solved by $w(\cdot,t)$ at each time. Subsequently, we take

$R:=\bigcup_{t>0}R_{t}\hbox{ and }\Sigma:=\bigcup_{t>0}\Sigma_{t},$

so that

$\mathcal{O}=\\{0<T(x)<\infty\\}=R\cup\Sigma.$

Let us also mention that $\\{R_{t}\\}_{t>0}$ is a foliation of $R$, and so is
$\\{\Sigma_{t}\\}_{t>0}$ for $\Sigma$, due to Proposition 2.12. We further
subdivide singular points into strata by the dimensionality of the zero set;
specifically, for $0\leq k\leq d-1$ we denote

$\Sigma^{k}_{t}:=\\{x\in\partial\Omega_{t}:\dim\ker D^{2}w=k\\}\hbox{ and
}\Sigma^{k}:=\bigcup_{t}\Sigma^{k}_{t}.$

For the obstacle problem, regular points are relatively open in the free
boundary ([Bla00], Corollary 4.8), and thus singular points form a closed set.
This topological control is lost in the union over all times, so a first step
is to reestablish that control for our $R$ and $\Sigma$. For this, we use a
lemma due to Blank, which allows us to identify regular points by finite-scale
behavior.

######  Lemma 5.3 ([Bla00] Theorem 4.5).

Let $u\geq 0$ solve $\Delta u=f\chi_{\\{u>0\\}}$ in $B_{1}(0)$, with
$0\in\partial\\{u>0\\}$ and $0<f\leq 1$. Then there exist universal parameters
$\lambda_{0},r_{0},\tau\in(0,1)$ such that if $\lambda_{0}<f$ in $B_{1}$ and

$\frac{|\\{x:u(x)=0\\}\cap B_{r}(x_{0})|}{|B_{r}|}\geq\frac{1}{8}\hbox{ for
some }r<r_{0},$

then

$\frac{|\\{x:u(x)=0\\}\cap B_{s}(x_{0})|}{|B_{s}|}\geq\frac{3}{8}\hbox{ for
all }s<\tau r.$

In particular, if the hypothesis of the lemma holds, then $\\{u=0\\}$ has
positive density at 0, so 0 is a regular point. This lemma is applicable even
when $f$ is less regular than Hölder, and results in a modified regular-
singular dichotomy in that case. Essentially, one may take away that
nonuniqueness of the blowup limit in the low regularity setting can occur due
to infinite rotation, but not due to any sort of mixing of regular and
singular point behavior at different scales. Since the following result only
relies on the previous lemma and continuity of $T$, it also holds when $T$ is
less regular than Hölder.

######  Proposition 5.4.

$R$ is open. Thus, $\Sigma$ is relatively closed in $\mathcal{O}$.

###### Proof.

Let $\lambda_{0},r_{0},\tau\in(0,1)$ be the parameters given by Lemma 5.3.
Suppose $x_{0}\in R$, so that $\Omega_{T(x_{0})}$ has density $\frac{1}{2}$ at
$x_{0}$, and thus there exists $r>0$ such that

$\frac{|\\{x:w(x,T(x_{0}))=0\\}\cap B_{r}(x_{0})|}{|B_{r}|}\geq\frac{1}{8}$

Then the result of the lemma is that for all $s<\tau r$,

$\frac{|\\{x:T(x)\geq T(x_{0})\\}\cap B_{s}(x_{0})|}{|B_{s}|}\geq\frac{3}{8}$

Fix $s_{0}=\tau r/2$. By Lemma 2.2, $|\Omega_{t}|$ is continuous in $t$, so we
can choose $\delta$ for $s_{0}$ such that

$\frac{|\\{x:T(x)>T(x_{0})+\delta\\}\cap
B_{s_{0}}(x_{0})|}{|B_{s_{0}}|}\geq\frac{5}{16}$

Let $s_{1}<s_{0}$ be such that if $|x-x_{0}|<s_{1}$, then
$|T(x)-T(x_{0})|<\delta$. Then for $x_{1}\in B_{s_{1}}(x_{0})$, we compute

$\displaystyle\frac{|\\{x:T(x)>T(x_{1})\\}\cap
B_{s_{0}+|x_{1}-x_{0}|}(x_{1})|}{|B_{s_{0}+|x_{1}-x_{0}|}|}$
$\displaystyle\geq\frac{|\\{x:T(x)>T(x_{0})+\delta\\}\cap
B_{s_{0}}(x_{0})|}{|B_{s_{0}+|x_{1}-x_{0}|}|}$
$\displaystyle\geq\left(\frac{5}{16}\right)\left(1-\frac{s_{0}}{s_{0}+s_{1}}\right)^{d}$

Thus, if we take $s_{1}$ sufficiently small, this last quantity is greater
than $\frac{1}{8}$, and all points in $B_{s_{1}}(x_{0})$ are regular. ∎

###  5.1. Regular points

We now turn toward understanding the behavior of the interface near regular
points. Standard obstacle problem theory ([Caf98], [Bla00]) gives that for
$C^{0,\alpha}$ source, the interface is locally $C^{1,\alpha}$ at regular
points, with the scale at which the regularity is achieved depending on the
scale at which the zero set is sufficiently large. We discuss the dependence
of this regularity in greater detail in the appendix. In particular, for this
problem we have:

######  Proposition 5.5.

$R$ can be covered by open neighborhoods $V$, each with the property that
there exist constants $C,r>0$ such that for each $x\in V$,
$B_{r}(x)\cap\Omega_{T(x)}$ is the intersection of $B_{r}(x)$ with the lower
graph of a $C^{1,\alpha}$ function in some coordinate system (depending on
$x$) with seminorm bounded by $C$.

###### Proof.

By Lemma 6.9, we need only show that if the zero set reaches density
sufficiently close to $\frac{1}{2}$ at scale $r$ near $x$, then it does so at
the same scale at all points near $x$. This is essentially immediate from the
fact that $\Omega_{t}$ expands monotonically in $t$, with the measure
$|\Omega_{t}|$ Lipschitz as a function of $t$ by Lemma 2.2. ∎

The dependence of the coordinate system in Proposition 5.5 is only a minor
inconvenience, and we will eventually remove it in Proposition 5.12. To better
understand this dependence, we introduce $\nu(x)$, defined for $x\in R$ as the
outward unit normal to $\Omega_{T(x)}$ at $x$. We have spatial regularity of
$\nu$ from the obstacle problem; namely, $\nu\in
C^{0,\alpha}_{\operatorname{\textup{loc}}}(R_{t})$ for each $t$. Our goal will
be to improve this to regularity of $\nu$ on $R$.

The key ingredients will be the regularity of $\Omega_{t}$ near regular
points, and the strictly monotonic expansion of the $\Omega_{t}$. The
essential idea will be that if the tangent planes to $\Omega_{T(x)}$ at $x$
and to $\Omega_{T(y)}$ at $y$ intersect for some points $x,y$ with different
hitting times, they must intersect well away from $x$ and $y$ or else we will
be able to use the regularity of the interfaces to show that
$\partial\Omega_{T(x)}$ and $\partial\Omega_{T(y)}$ intersect, which
contradicts monotonicity. This then gives control over the angle at which the
tangent planes may intersect in terms of the distance between $x$ and $y$.

######  Proposition 5.6.

The outward unit normal vector $\nu$ to $\Omega_{T(x)}$ at $x$ satisfies
$\nu\in C^{0,\alpha/(1+\alpha)}_{\operatorname{\textup{loc}}}(R)$.

###### Proof.

By Proposition 5.5, we may cover $R$ with neighborhoods $V$ such that for each
$x\in V$ uniformly, $\partial\Omega_{T(x)}\cap V$ is the lower graph in some
coordinate system of a function $f_{T(x)}$, with the $f_{t}$ uniformly bounded
in $C^{1,\alpha}$. We will restrict to such a $V$ for the remainder of the
proof.

Then as a preliminary step, we can observe continuity of $\nu$ from the
regularity and monotonicity of the interface by a purely geometrical argument.
Namely, a $C^{1,\alpha}$ domain entertains a uniform interior and exterior
cone condition, where the angle of the cone improves toward $\pi$ as we allow
its height to approach 0; specifically, the cone in $B_{r}$ can be taken with
angle $2\arccos(Cr^{\alpha})$, when the $C^{1,\alpha}$ seminorm is $C$. Thus,
if $\nu$ were discontinuous at some $x\in R$, we could use compactness to find
a sequence $(y_{n})$ converging to $x$ with $T(y_{n})$ either increasing or
decreasing to $T(x)$ and $\nu(y_{n})$ converging to some unit vector distinct
from $\nu(x)$. Then for $n$ sufficiently large, at a sufficiently small scale,
the interior cone at $x$ will intersect with the exterior cone of a $y_{n}$,
or vice versa, and we draw a contradiction with the monotonic expansion of the
$\Omega_{t}$ depending on whether the $T(y_{n})$ are decreasing or increasing.

Then, we have checked that $x\to\nu(x)$ is continuous. Now, to obtain a
quantitative local continuity estimate in view of the cone regularity we
described above, we may restrict attention to $x,y\in V$ with
$\frac{1}{2}<\nu(x)\cdot\nu(y)<1$. Moreover, since the case $T(x)=T(y)$ is
managed by the spatial regularity of the interface, we may assume that
$T(x)>T(y)$. For notation, we let $r=\nu(x)\cdot\nu(y)$ and use $P_{x},P_{y}$
to refer to the tangent planes to $\Omega_{T(x)}$ at $x$ and $\Omega_{T(y)}$
at $y$ respectively.

Let $v:=\frac{\nu(y)-|r|\nu(x)}{\sqrt{1-r^{2}}}$ be the projection of $\nu(y)$
into $\nu(x)^{\perp}$ , scaled to unit norm. Considering the point $x-hv$ for
$h>0$, we compute that its $\nu(y)$ component is
$x\cdot\nu(y)-h\sqrt{1-r^{2}}$. Thus $x-hv$ reaches $P_{y}$ precisely when
$h=\frac{(x-y)\cdot\nu(y)}{\sqrt{1-r^{2}}}$, and in general we have

(5.1) $d(x-hv,P_{y})\geq h\sqrt{1-r^{2}}-\delta\quad\hbox{ where
}\delta:=|x-y|$

Now, we apply the $C^{1,\alpha}$ regularity of the interface in $V$. For all
$h$ sufficiently small, this regularity implies that $\partial\Omega_{T(x)}$
in $B_{h}(x)$ is contained in a $Ch^{1+\alpha}$-neighborhood of $P_{x}$; in
other words, $\Omega_{T(x)}$ and its exterior contain the following
halfspaces:

(5.2) $\\{z\in B_{h}(x):(z-x-Ch^{1+\alpha})\cdot\nu(x)\leq
0\\}\subset\Omega_{T(x)}\cap B_{h}(x)$ (5.3) $\\{z\in
B_{h}(x):(z-x+Ch^{1+\alpha})\cdot\nu(x)\geq
0\\}\subset({\mathbb{R}}^{d}\setminus\Omega_{T(x)})\cap B_{h}(x)$

In particular, since $x-hv\in P_{x}$, it follows that there is a point
$\tilde{x}\in\partial\Omega_{T(x)}$ with

(5.4) $|\tilde{x}-(x-hv)|<Ch^{1+\alpha}.$

Let $y_{1}$ be the nearest point in $P_{y}$ to $\tilde{x}$. We illustrate this
with the figure below.

Figure 1.

As the figure may suggest, $\tilde{x}$ cannot be too far below $y_{1}$, or
else it falls into $\Omega_{T(y)}$, contradicting that $T(x)>T(y)$.
Specifically, as in (5.2), we can apply the $C^{1,\alpha}$ regularity to to
get that for all sufficiently small $r$, $\Omega_{T(y)}\cap B_{r}(y)$ contains
the halfspace $\\{z:(z-y-Cr^{1+\alpha}\nu(y))\cdot\nu(y)\leq 0\\}\cap
B_{r}(y)$. Since $\tilde{x}\notin\overline{\Omega_{T(y)}}$, $\tilde{x}$ must
not be contained in that halfspace, and we have

(5.5) $(\tilde{x}-y_{1})\cdot\nu(y)>-C|y_{1}-y|^{1+\alpha}$

The left side here is the signed distance of $\tilde{x}$ to $P_{y}$. From
(5.1), the signed distance of $x-hv$ to $P_{y}$ is bounded above by
$-(h\sqrt{1-r^{2}}-\delta)$, so using (5.4), we conclude that the left side
above is bounded above by $Ch^{1+\alpha}-(h\sqrt{1-r^{2}}-\delta)$. On the
other hand, we have

$|y_{1}-y|\leq|\tilde{x}-y|\leq|\tilde{x}-(x-hv)|+|(x-hv)-x|+|x-y|\leq
C_{0}h^{1+\alpha}+h+\delta$

When $\delta<h<1$, this is $O(h)$, and so $|y_{1}-y|^{1+\alpha}\geq-
Ch^{1+\alpha}$ for some $C$. Thus, the inequality (5.5) becomes

$Ch^{1+\alpha}-(h\sqrt{1-r^{2}}-\delta)>-Ch^{1+\alpha}$

Rearranging and absorbing constants, this means

$h\sqrt{1-r^{2}}-\delta\leq Ch^{1+\alpha}$

so that

$|\nu(x)-\nu(y)|\leq\sqrt{1-r^{2}}\leq Ch^{\alpha}+\delta h^{-1}$

Optimizing $h$, we get

$|\nu(x)-\nu(y)|\leq C\delta^{\frac{\alpha}{1+\alpha}}$

where the constant depends only on the uniform bound for the $C^{1,\alpha}$
seminorms of the graphs, and on $\alpha$. Hence we conclude. ∎

As a result of the regularity of $\nu$, we can improve Proposition 5.5 to also
have the coordinate system chosen locally uniformly. In other words, near
regular points, one can fix a local coordinate system in which the free
boundary evolves as a $C^{1,\alpha}$ graph over some time interval.

We now turn toward applying the improved geometry of the patch at regular
points to the pressure. Elliptic regularity for $C^{1,\alpha}$ domains implies
that the pressure $p(\cdot,t)$ has a well-defined gradient on $R_{t}$, and the
Hopf lemma for $C^{1,\alpha}$ domains implies that $\nabla p(\cdot,t)$ is
nonvanishing on $R_{t}$. However, there is an important limitation here:
$x\mapsto\nabla p(x,T(x))$ is not necessarily continuous on $R$, complicating
our analysis. This is illustrated with the following example:

######  Remark 5.7.

Singular points can exert a nonlocal effect on the pressure gradient.

For example, if we consider a pressure supported on a strip of width $h$ with
zero boundary conditions and constant Laplacian $-1$, then we can see that
$|\nabla p|=\frac{h}{2}$ on the boundary, since the solution to the one-
dimensional problem with $p(0)=p(h)=0$ is $p(x)=\frac{1}{2}x(h-x)$. In
particular, it follows that if we have a patch which consists of two strips,
and those strips merge along a hyperplane at some time, then $\nabla p$ has a
jump discontinuity in time at every regular point at the time those strips
merge.

Similar examples can be considered for singular points in each stratum
$\Sigma^{k}$ by examining a cylindrical patch with a cylindrical hole as the
radius of the cylinder shrinks to 0. Here, by cylinder we mean the product of
${\mathbb{R}}^{k}$ and a $(d-k)$-dimensional ball, for $0\leq k\leq d-1$.

We note that the discontinuity in the gradient in Remark 5.7 is a jump in
magnitude, not in direction. Indeed, the zero boundary condition implies that
$\frac{\nabla p(x,T(x))}{|\nabla p(x,T(x))|}=-\nu(x)$ on $R_{t}$, so the
regularity of $\nu$ from Proposition 5.6 rules out such discontinuities.

This example also proves to be an obstacle to higher regularity of $T$; we
will later see in Lemma 5.18 that the derivative of $T$ closely depends on
$\nabla p(x,T(x))$ when it exists. As a result, we can hope for $T$ to be at
best Lipschitz on $R$.

The key in establishing this regularity will be to obtain a quantitative
estimate from the Hopf lemma, to get a locally uniform lower bound for
$|\nabla p(x,T(x))|$ on $R$. For this, we require an a priori estimate for the
growth of the solution, in addition to control over the geometry. We can
obtain this growth estimate from the strict superharmonicity of $p(\cdot,t)$,
and so we have the following statement:

######  Lemma 5.8.

Let $r_{0},c_{0},C_{0}>0$ and $\alpha\in(0,1)$. Suppose $u$ is a positive
$C^{2}$ solution to $\Delta u\leq-c_{0}<0$ on
$B_{r_{0}}(0)\cap\\{x:x_{d}>C_{0}|x^{\prime}|^{1+\alpha}\\}$ with $u(0)=0$,
where we write $x=(x^{\prime},x_{d})$. Then there exist $\varepsilon,\delta>0$
depending only on $d,\alpha,r_{0},c_{0},C_{0}$ such that for all
$h\in(0,\delta)$, we have $u(he_{n})\geq\varepsilon h$.

###### Proof.

First, we let $C_{1}=C_{0}+1$, to give additional separation from the boundary
when we are away from $0$, and let
$U_{r}=B_{r}(0)\cap\\{x:x_{d}>C_{1}|x^{\prime}|^{1+\alpha}\\}$ for
$0<r<r_{0}$. For a given $r$, we will decompose $\partial U_{r}$ as
$\Gamma_{1}\cup\Gamma_{2}$, where $\Gamma_{2}:=\partial
B_{r}(0)\cap\\{x:x_{d}\geq C_{1}|x^{\prime}|^{1+\alpha}\\}$ is a spherical
cap, and $\Gamma_{1}=B_{r}(0)\cap\\{x:x_{d}=C_{1}|x^{\prime}|^{1+\alpha}\\}$.

Figure 2. the region $U_{r}$ (shaded)

The proof of the Hopf lemma proceeds by perturbing $u$ by a function $v$
constructed specially for the domain and applying the comparison to the
result. In particular, if for some $v,\varepsilon>0,\delta>0$ we have

$\begin{cases}\Delta(u-\varepsilon v)\leq 0\hbox{ on }U_{\delta}\\\
u-\varepsilon v\geq 0\hbox{ on }\partial U_{\delta}\\\
\partial_{d}v(0)=1\end{cases}$

then the comparison principle implies that $u-\varepsilon v\geq 0$ on
$U_{\delta}$, and the control over the derivative of $v$ at 0 implies the
result.

Borrowing from the proof of the Hopf lemma for $C^{1,\alpha}$ domains in
[LN07], we let

$v(x)=x_{d}+\frac{2C_{1}}{\alpha}(\alpha+d-1)x_{d}^{1+\alpha}-2C_{1}|x|^{1+\alpha}$

We note by direct computation that $\partial_{d}v(0)=1$ and

$\Delta
v(x)=2C_{1}(1+\alpha)(\alpha+d-1)\left(x_{d}^{\alpha-1}-|x|^{\alpha-1}\right)$

In particular, $\Delta v(x)\geq 0$ on $U_{r}$ for any $r$.

Thus, we reduce to verifying the boundary condition, which requires using the
particular behavior of $v$ on $\Gamma_{1}$, and choosing $\delta$ and
$\varepsilon$ appropriately to control $v$ on $\Gamma_{1}$ and $\Gamma_{2}$.
Following this plan, we first check the boundary inequality on $\Gamma_{1}$,
defined above, where we have $x_{d}=C_{1}|x^{\prime}|^{1+\alpha}$. Along a
curve $x(t)=(te^{\prime},C_{0}t^{1+\alpha})$ for a unit vector
$e^{\prime}\in{\mathbb{R}}^{d-1}$, we have

$v(x(t))=C_{1}t^{1+\alpha}+\frac{2C_{0}^{2}}{\alpha}(\alpha+d-1)t^{(1+\alpha)^{2}}-2C_{0}(t^{2}+C_{1}t^{2(1+\alpha)})^{(1+\alpha)/2}$

For small $t>0$, we can drop the higher order terms and see that $v$ grows
like $C_{1}t^{1+\alpha}-2C_{1}t^{1+\alpha}=-C_{1}t^{1+\alpha}$. In particular,
there exists $\delta_{0}=\delta_{0}(d,\alpha,C_{0})$ such that if
$\delta\leq\delta_{0}$, then $v(x)\leq 0$ on
$B_{\delta}(0)\cap\\{x:x_{d}=C_{1}|x^{\prime}|^{1+\alpha}\\}$. We will set
$\delta=\min(\delta_{0},r_{0}/2)$ for the rest of the proof.

Next, we handle the boundary inequality on $\Gamma_{2}$, the spherical cap
defined by $\Gamma_{2}=\partial B_{\delta}(0)\cap\\{x:x_{d}\geq
C_{1}|x^{\prime}|^{1+\alpha}\\}$. By compactness, $\Gamma_{2}$ has positive
distance to $\\{x:x_{d}\geq C_{0}|x^{\prime}|^{1+\alpha}\\}$. Let
$D=D(d,\alpha,\delta,C_{0})$ denote this distance. By shrinking $D$ to
$\frac{r_{0}}{2}$ if necessary, we get that $u$ is defined on a ball of radius
$D$ at each point of $\Gamma_{2}$. The paraboloid on the ball of radius $D$
with zero boundary data and Laplacian $-c_{0}$ is a subsolution to $u$, so we
conclude that $u(x)\geq\frac{c_{0}}{2}D^{2}$ on $\Gamma_{2}$. Let
$M=\max(1,\max_{\Gamma_{2}}v)$, and then we can take

$\varepsilon=\frac{c_{0}D^{2}}{2M}$

to get $u-\varepsilon v\geq 0$ on $\Gamma_{2}$, with $\varepsilon$ depending
on all of the parameters in the statement of the lemma.

On $\Gamma_{1}$, we have $u\geq 0$ and $v\leq 0$, so $u-\varepsilon v\geq 0$.
On $\Gamma_{2}$, we chose $\varepsilon$ so that $u-\varepsilon v\geq 0$. Thus,
$u-\varepsilon v\geq 0$ on $\partial U_{\delta}$, so we conclude.

∎

With the Hopf lemma estimate and the regularity of the boundary near regular
points, we can conclude that $p(\cdot,T(x))$ has linear growth at $x$, locally
uniformly on $R$.

######  Proposition 5.9.

$R$ can be covered with neighborhoods $V$ with the following property: there
exist parameters $C,c,r_{0}>0$ such that for any $x\in V$,

$c\leq r^{-1}\sup_{B_{r}(x)}p(\cdot,T(x))\leq C$

In particular, $|\nabla p(x,T(x))|\sim 1$ on $V$, for implicit constants
depending on $V$.

###### Proof.

To obtain the lower bound, we apply Lemma 5.8. Thus, we must show that the
parameters $r_{0},c_{0},C_{0}$ from the statement of the lemma can be chosen
locally uniformly on $R$. By Lemma 2.4, we can choose $c_{0}>0$ locally
uniformly in time so that $\Delta p=-n\leq-c_{0}<0$, using the assumption that
$n_{0}$ is bounded away from 0. By Proposition 5.5, for any $x_{0}\in R$, we
can find a neighborhood $V$ of $x_{0}$ in which we have a uniform $C_{0}$ so
that each $R_{t}$ which intersects the neighborhood does so as a graph with
$C^{1,\alpha}$ seminorm controlled by $C_{0}$. By taking $r_{0}$ so that
$B_{2r_{0}}(x_{0})\subset V$, we can use $r_{0},C_{0}$ as the parameters for
all points in $B_{r_{0}}(x_{0})$.

To obtain the upper bound, we first use Proposition 5.5. This gives us a
neighborhood $V$ and a parameter $r_{1}>0$ for which
$B_{r_{1}}(x)\cap\partial\Omega_{T(x)}$ is a $C^{1,\alpha}$ graph, with
uniform control over the $C^{1,\alpha}$ seminorm. In particular, this
regularity implies that there is an $r_{2}>0$ such that for all $x\in V$ and
all $r\in(0,r_{2}]$, we have $x-r\nu(x)\in\Omega_{T(x)}$. Then the mean value
theorem along the path $x-r\nu(x)$ implies that there exists some $r$ for
which

$\nabla p(x,x-r\nu(x))\cdot-\nu(x)=\frac{p(x-r_{2}\nu(x),T(x))}{r_{2}}$

We recall from the proof of Lemma 2.3 that $p\in
L^{\infty}({\mathbb{R}}^{d}\times[0,\tau])$, for any $\tau\in(0,\infty)$,
using that the patch has bounded support and comparing to a sufficiently large
paraboloid supersolution. Thus, we have

$|\nabla p(x,x-r\nu(x))\cdot\nu(x)|\leq Cr_{2}^{-1}$

Using the regularity of the boundary, the $L^{\infty}$ bound on the pressure,
and the $L^{\infty}$ bound on the nutrient from Lemma 2.2, we can invoke
boundary Schauder estimates to have $p(\cdot,T(x))$ uniformly $C^{1,\alpha}$
on $B_{r_{1}}(x)\cap\Omega_{T(x)}$ for each $x\in V$. Thus, we can transfer
our bound to the boundary:

$|\nabla p(x,T(x))|=|\nabla p(x,T(x))\cdot\nu(x)|\leq|(\nabla p(x,T(x))-\nabla
p(x-r\nu(x)))\cdot\nu(x)|+|\nabla p(x-r\nu(x))\cdot\nu(x)|\leq Cr+Cr_{2}^{-1}$

This also gives the upper bound on the linear growth of $p(\cdot,T(x))$ near
$x$, so we conclude. ∎

Having established nondegeneracy of the pressure, we get improved regularity
of $T$ on $R$.

######  Corollary 5.10.

$T\in C^{0,1}_{\operatorname{\textup{loc}}}(R)$. In other words, $T$ attains
its optimal regularity on $R$ in light of Remark 5.7, barring additional
assumptions on $\Sigma$.

###### Proof.

This follows from the boundary regularity from Proposition 5.5 and the linear
nondegeneracy of the pressure from Proposition 5.9, via a method similar to
the comparison arguments of Section 4. Using these properties, we can
construct a radial subsolution initially supported on an annulus in
$\Omega_{t}$ which expands at a constant rate, near any $x_{0}\in R$. From
this argument, we get

(5.6) $(T(x)-T(x_{0}))_{+}\leq C|x-x_{0}|$

for $x_{0}\in R$ and $x$ sufficiently close to $x_{0}$. Since the boundary
regularity and linear growth rate are uniform for $x_{0}$ restricted to a
compact $K\subset R$, the one-sided bound (5.6) holds uniformly for
$x,x_{0}\in K$, and so we conclude that $T$ is locally Lipschitz on $R$. ∎

Finally, we turn to refining the statement of Proposition 5.5. The first step
will be to use the linear growth of the pressure gradient to control the
boundary in the Hausdorff metric.

######  Lemma 5.11.

For any $x\in R$, there exists parameters $r,\delta>0$ such that for
$t_{1},t_{2}\in(T(x)-\delta,T(x)+\delta)$, we have

$\sup_{y_{1}\in R_{t_{1}}\cap B_{r}(x)}\inf_{y_{2}\in R_{t_{2}}\cap
B_{r}(x)}|y_{1}-y_{2}|\sim|t_{1}-t_{2}|$

In other words, $D(R_{t_{1}}\cap B_{r}(x),R_{t_{2}}\cap
B_{r}(x))\sim|t_{1}-t_{2}|$, where $D$ denotes Hausdorff distance. Here, all
parameters and implicit constants depend on $x$.

###### Proof.

Let $B$ be a ball compactly contained in $R$, and let $y_{1},y_{2}\in B$ with
$T(y_{1})<T(y_{2})$. Since $T\in C^{0,1}_{\operatorname{\textup{loc}}}(R)$ by
Corollary 5.10, we have $|T(y_{1})-T(y_{2})|\leq C(B)|y_{1}-y_{2}|$. We get
the reverse bound by an analogous argument using the boundary regularity of
Proposition 5.5 and using Proposition 5.9’s upper bound on the linear growth
of $p(\cdot,t_{1})$ away from $y_{1}$. That is, following the approach of
Section 4, we can construct a supersolution supported outside a ball in the
exterior of $\Omega_{T(y_{1})}$ near $y_{1}$, such that the supersolution
expands at a constant rate. From that argument, we get that there exists a
point $\tilde{y_{2}}\in R_{t_{2}}$ with $|y_{1}-\tilde{y_{2}}|\leq
C(B)|T(y_{1})-T(y_{2})|$. Letting $t_{1}=T(y_{1}),t_{2}=T(y_{2})$, and
restricting to the case where $|t_{1}-t_{2}|$ is sufficiently small to
guarantee that the $\tilde{y_{2}}$ from before is in $B$, we have shown that
$D(R_{t_{1}}\cap B,R_{t_{2}}\cap B)\sim|t_{1}-t_{2}|$ for implicit constants
depending on $B$, from which we can obtain the original statement. ∎

Using the previous result, we can now state our final improved form of
Proposition 5.5.

######  Proposition 5.12.

$R$ is covered by neighborhoods $V$ with the following property: there exists
$r>0$, a coordinate system $(x^{\prime},x_{n})$, and a locally defined
function $f(x^{\prime},t)$ such that for each $x\in V$, $\Omega_{T(x)}\cap
B_{r}(x)$ is the lower graph $\\{y\in B_{r}(x):y_{n}\leq
f(y^{\prime},T(x))\\}$. Moreover, $f$ is uniformly $C^{1,1}$ in space and
$C^{0,1}$ in time.

###### Proof.

First, by Proposition 5.6, we note that the coordinate system in Proposition
5.5 can be chosen locally uniformly, so that we get $r>0$ and the family
$f(x^{\prime},t)$ which are uniformly $C^{1,1}$ in space by the Lipschitz
regularity of $T$.