# Regularity properties of passive scalars with rough divergence-free drifts

Dallas Albritton  Courant Institute of Mathematical Sciences, New York
University, New York, NY 10012<EMAIL_ADDRESS>and  Hongjie Dong
Division of Applied Mathematics, Brown University, 182 George Street,
Providence, RI 02912<EMAIL_ADDRESS>
###### Abstract.

We present sharp conditions on divergence-free drifts in Lebesgue spaces for
the passive scalar advection-diffusion equation

$\partial_{t}\theta-\Delta\theta+b\cdot\nabla\theta=0$

to satisfy local boundedness, a single-scale Harnack inequality, and upper
bounds on fundamental solutions. We demonstrate these properties for drifts
$b$ belonging to $L^{q}_{t}L^{p}_{x}$, where $\frac{2}{q}+\frac{n}{p}<2$, or
$L^{p}_{x}L^{q}_{t}$, where $\frac{3}{q}+\frac{n-1}{p}<2$. For steady drifts,
the condition reduces to $b\in L^{\frac{n-1}{2}+}$. The space
$L^{1}_{t}L^{\infty}_{x}$ of drifts with ‘bounded total speed’ is a borderline
case and plays a special role in the theory. To demonstrate sharpness, we
construct counterexamples whose goal is to transport anomalous singularities
into the domain ‘before’ they can be dissipated.

###### Contents

  1. 1 Introduction
  2. 2 Local boundedness and Harnack’s inequality
  3. 3 Bounded total speed
  4. 4 Counterexamples
  5. 5 Upper bounds on fundamental solutions

##  1\. Introduction

We consider the linear advection-diffusion equation

$\partial_{t}\theta-\Delta\theta+b\cdot\nabla\theta=0.$ (A-D)

The solution $\theta=\theta(x,t)$ is known as a _passive scalar_ , and the
prescribed divergence-free velocity field $b=b(x,t)$ is known as the _drift_.

Our goal is to understand, in detail, the regularity properties of (A-D) when
the drift $b$ is _rough_. Rough divergence-free drifts arises naturally in the
context of nonlinear PDEs in fluid dynamics.

To understand what is ‘rough’, we recall the scaling symmetry

$u\to u(\lambda x,\lambda^{2}t),\quad b\to\lambda b(\lambda
x,\lambda^{2}t),\quad\lambda>0.$ (1.1)

In dimensional analysis, one writes $[x]=L$, $[t]=L^{2}$, and $[b]=L^{-1}$.
The scaling (1.1) identifies the Lebesgue spaces $L^{q}_{t}L^{p}_{x}$, where
$\frac{2}{q}+\frac{n}{p}\leq 1$, as _(sub)critical_ spaces for the drift,
meaning spaces whose norms do not grow upon ‘zooming in’ with the scaling
symmetry. For example,

$X=L^{\infty}_{t}L^{n}_{x},\,L^{2}_{t}L^{\infty}_{x},\,L^{n+2}_{t,x}$ (1.2)

are _critical spaces_ , whose norms are dimensionless, i.e., invariant under
the symmetry (1.1). Here, $n\geq 2$ is the spatial dimension.

When $b$ belongs to one of the above critical Lebesgue spaces, it is not
difficult to adapt the work of De Giorgi, Nash, and Moser [DG57, Nas58, Mos64]
to demonstrate that weak solutions of (A-D) are Hölder continuous and satisfy
Harnack’s inequality. The above threshold is known to be _sharp_ for Hölder
continuity within the scale of Lebesgue spaces, see the illuminating
counterexamples in [SVZ13, Wu21]. The divergence-free condition moreover
allows access to drifts in the critical spaces
$X=L^{\infty}_{t}L^{-1,\infty}_{x},L^{\infty}_{t}{\rm BMO}^{-1}_{x}$,
considered by [Osa87, SSvZ12, QX19b]. In these spaces, it is furthermore
possible to prove Gaussian upper and lower bounds on fundamental solutions in
the spirit of Aronson [Aro67].

In this paper, we are concerned with _supercritical_ drifts, for which
continuity may fail. Nonetheless, much of the regularity theory may be
salvaged. The divergence-free structure plays a crucial role here,111The
divergence-free structure plays a more subtle role in the critical case.
Without this structure, the drift is required to be small in a critical
Lebesgue space or Kato class, and local boundedness may depend on the
‘profile’ of the drift. as is already visible from the computation

$\int(b\cdot\nabla\theta)\theta\phi^{2}\,dx\,dt=\int
b\cdot\nabla\left(\frac{\theta^{2}}{2}\right)\phi^{2}\,dx\,dt=-\int\theta^{2}(b\cdot\nabla\phi)\phi.$
(1.3)

With this well known observation, one may apply Moser’s iteration scheme to
demonstrate that, when $b\in L^{q}_{t}L^{p}_{x}$ and
$\frac{2}{q}+\frac{n}{p}<2$, solutions are _locally bounded_ , see [NUt11].
Typical examples are

$X=L^{\infty}_{t}L^{\frac{n}{2}+}_{x},L^{\frac{n+2}{2}+}_{t,x}.$ (1.4)

Under these conditions (and a weak background assumption), the Harnack
inequality persists as a _single-scale Harnack inequality_ [IKR14, IKR16]: In
the steady case $\theta=\theta(x)$,

$\sup_{B_{R}}\theta\leq C_{R}\inf_{B_{R}}\theta,$ (1.5)

where $C_{R}$ may become unbounded as $R\to 0^{+}$. Whereas a scale-invariant
Harnack inequality implies Hölder continuity, it is perhaps less well known
that a single-scale Harnack may hold in the absence of Hölder continuity.
Finally, in this setting, pointwise upper bounds on fundamental solutions
continue to hold, although they have ‘fat tails’ compared to their Gaussian
counterparts [Zha04, QX19a].

Our first contribution is to understand the sharpness of the condition $b\in
L^{q}_{t}L^{p}_{x}$, $\frac{2}{q}+\frac{n}{p}<2$. We find that
$L^{1}_{t}L^{\infty}_{x}$, the space of drifts with ‘bounded total speed’ (in
the terminology of [Tao13]), plays a special role and informs the
counterexamples we construct in Section 4. We summarize the results pertaining
to this condition in Theorem 1.1.

In the steady case, there is an additional subtle feature, which is not well
known and that we find surprising: Local boundedness continues to hold when
$b\in L^{\frac{n-1}{2}+}$. To our knowledge, this ‘dimension reduction’ was
first observed in this context by Kontovourkis [Kon07] in his thesis.222 This
‘dimension reduction’ itself goes back at least to work of Frehse and
R$\overset{\circ}{\text{u}}$žička [FR96] on the steady Navier-Stokes equations
in $n=6$. The ‘slicing’ is also exploited by Struwe in [Str88]. Heuristically,
Kontovourkis’ key observation is as follows. Consider the basic $L^{2}$ energy
estimate, in a ball $B_{r}$ without smooth cut-off. The drift contributes the
boundary term

$\int_{B_{r}}(b\cdot\nabla\theta)\theta\,dx=\int_{\partial
B_{r}}\frac{\theta^{2}}{2}b\cdot n\,d\sigma,$ (1.6)

where $d\sigma$ is the surface area measure. Since $\nabla\theta\in
L^{2}(B_{R})$, on ‘many slices’ $r\in(R/2,R)$, we have $\nabla\theta\in
L^{2}(\partial B_{r})$, with a quantitative bound. Similarly, $b$ belongs to
$L^{\frac{n-1}{2}+}(\partial B_{r})$ on ‘many slices’. Thus, one may exploit
Sobolev embedding on the sphere $\partial B_{r}$ to estimate the boundary
term. This dimension reduction was recently rediscovered by Bella and
Schäffner [BS21], who proved local boundedness and a single-scale Harnack
inequality in the context of certain degenerate elliptic PDEs, which we review
below.

Since the work of Kontovourkis, it has been an interesting question, what
dimension reduction holds in the parabolic setting? In particular, is $b\in
L^{\frac{n+1}{2}+}_{t,x}$ enough for local boundedness? Very recently, X.
Zhang [Zha20] generalized the work [BS21] of Bella and Schäffner to the
parabolic setting, and among other things, demonstrated local boundedness
under the condition $b\in L^{p}_{x}L^{q}_{t}$, $\frac{3}{q}+\frac{n-1}{p}<2$,
$p\leq q$, see Corollary 1.5 therein. Crucially, the order of integration is
reversed. The condition $b\in L^{\frac{n-1}{2}+}_{x}L^{\infty}_{t}$ implies
the elliptic case of Kontovourkis. From this condition, we see that, perhaps,
one dimension is not ‘reduced’, but rather hidden into the time variable.

In Theorem 1.2, we present our second contribution, namely, (i) the parabolic
Harnack inequality and pointwise upper bounds on fundamental solutions in this
setting, and (ii) counterexamples which illustrate the meaning and sharpness
of the ‘dimension reduction’.

We now present the main results, which constitute a detailed picture of the
local regularity theory for the passive scalar advection-diffusion equation
(A-D) with supercritical drifts.

Let $n\geq 2$, $\Omega\subset\mathbb{R}^{n}$ be a bounded domain, and
$\Omega^{\prime}\subset\subset\Omega$ be a subdomain. Let $I=(S,T]$ and
$I^{\prime}=(S^{\prime},T^{\prime}]\subset I$ be finite intervals such that
$S<S^{\prime}$. Let $Q_{I}=\Omega\times I$ and
$Q_{I}^{\prime}=\Omega^{\prime}\times I^{\prime}$. Let $p,q\in[1,+\infty]$.

######  Theorem 1.1 ($b\in L^{q}_{t}L^{p}_{x}$).

_(Local boundedness)_ [NUt11] If

$\zeta_{0}:=\frac{2}{q}+\frac{n}{p}<2,$ (1.7)

then we have the following _quantitative local boundedness_ property: If
$\theta\in L^{1}(Q_{I})\cap C^{\infty}(Q_{I})$ satisfies the drift-diffusion
equation (A-D) in $Q_{I}$ with divergence-free drift $b\in C^{\infty}(Q_{I})$
and

$b\in L^{q}_{t}L^{p}_{x}(Q_{I}),$ (1.8)

then

$\sup_{Q_{I}^{\prime}}|\theta|\lesssim\lVert\theta\rVert_{L^{1}(Q_{I})},$
(1.9)

where the implied constant depends on $n$, $\Omega$, $\Omega^{\prime}$, $I$,
$I^{\prime}$, $p$, $q$, and $\lVert b\rVert_{L^{q}_{t}L^{p}_{x}(Q_{I})}$.

_(Single-scale Harnack)_ [IKR16] If, additionally, $b\in
L^{2}_{t}H^{-1}_{x}(Q_{I})$ and $\theta>0$, then we have the following
_quantitative Harnack inequality_ : If $I_{1},I_{2}\subset\subset I$ are
intervals satisfying $\sup I_{1}<\inf I_{2}$, then

$\sup_{\Omega^{\prime}\times I_{1}}\theta\lesssim\inf_{\Omega^{\prime}\times
I_{2}}\theta,$ (1.10)

where the implied constant depends on $n$, $\Omega$, $\Omega^{\prime}$, $I$,
$I_{1}$, $I_{2}$, $p$, $q$, $\lVert b\rVert_{L^{q}_{t}L^{p}_{x}(Q_{I})}$, and
$\lVert b\rVert_{L^{2}_{t}H^{-1}_{x}(Q_{I})}$.

_(Bounded total speed)_ If $(p,q)=(\infty,1)$, then the above quantitative
local boundedness property holds with constants depending on $b$ itself rather
than $\|b\|_{L^{1}_{t}L^{\infty}_{x}(Q_{I})}$. (The property is false without
this adjustment.)

_(Sharpness)_ Let $Q=B_{1}\times(0,1)$. There exist a smooth divergence-free
drift $b\in C^{\infty}(Q)$ belonging to $L^{q}_{t}L^{p}_{x}(Q)$ for all
$(p,q)\in[1,+\infty]^{2}$ with $2/q+n/p=2$, $(p,q)\neq(\infty,1)$, and
satisfying the following property. There exists a smooth solution $\theta\in
L^{\infty}_{t}L^{1}_{x}\cap C^{\infty}(Q)$ to the advection-diffusion equation
(A-D) in $Q$ with

$\sup_{B_{1/2}\times(0,T)}|\theta|\to+\infty\text{ as }T\to 1_{-}.$

In particular, the above quantitative local boundedness property fails when
$2/q+n/p=2$ and $q>1$.

_(Upper bounds on fundamental solutions)_ [QX19a] If the divergence-free drift
$b\in C^{\infty}_{0}(\mathbb{R}^{n}\times[0,+\infty))$ belongs to
$L^{q}_{t}L^{p}_{x}(\mathbb{R}^{n}\times\mathbb{R}_{+})$ and
$1\leq\zeta_{0}<2$, then the fundamental solution $\Gamma=\Gamma(x,t;y,s)$ to
the parabolic operator $L=\partial_{t}-\Delta+b\cdot\nabla$ satisfies

$\displaystyle\Gamma(x,t;y,s)$ $\displaystyle\leq
C(t-s)\left(\frac{1}{t-s}+\frac{M_{0}}{(t-s)^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}$
(1.11)
$\displaystyle\quad\times\left[\exp\left(-\frac{c|x-y|^{2}}{t-s}\right)+\exp\left(-\frac{c|x-y|^{1+\frac{1}{1+\alpha_{0}}}}{((t-s)M_{0})^{\frac{1}{1+\alpha_{0}}}}\right)\right]$

for all $x,y\in\mathbb{R}^{n}$ and $0\leq s<t<+\infty$. Here,

$\alpha_{0}=\frac{\zeta_{0}-1}{\theta_{2}},\quad
M_{0}=C\|b\|_{L^{q}_{t}L^{p}_{x}(\mathbb{R}^{n}\times\mathbb{R}_{+})}^{\frac{1}{\theta_{2}}}+\frac{1}{4},\quad\theta_{2}=1-\frac{\zeta_{0}}{2}.$
(1.12)

See Figure 1 for an illustration of Theorem 1.1.

Figure 1. Divergence-free drift $b\in L^{q}_{t}L^{p}_{x}$ in dimension $n\geq
2$ (dimension $n=3$ illustrated above). _Region A_ ($2/q+n/p\leq 1$): Local
boundedness, Harnack inequality, and Hölder continuity. _Region B_
($1<2/q+n/p<2$ or $(p,q)=(\infty,1)$): Local boundedness and single-scale
Harnack inequality. _Dashed line_ : Local boundedness is false.

######  Theorem 1.2 ($b\in L^{p}_{x}L^{q}_{t}$).

_(Local boundedness)_ [Zha20] If

$\zeta_{0}:=\frac{3}{q}+\frac{n-1}{p}<2,\quad p\leq q,$ (1.13)

then we have the following _quantitative local boundedness_ property: If
$\theta\in L^{1}(Q_{I})\cap C^{\infty}(Q_{I})$ satisfies the drift-diffusion
equation (A-D) in $Q_{I}$ with divergence-free drift $b\in C^{\infty}(Q_{I})$
and

$b\in L^{p}_{x}L^{q}_{t}(Q_{I}),$ (1.14)

then

$\sup_{Q_{I}^{\prime}}|\theta|\lesssim\lVert\theta\rVert_{L^{1}(Q_{I})},$
(1.15)

where the implied constant depends on $n$, $\Omega$, $\Omega^{\prime}$, $I$,
$I^{\prime}$, $p$, $q$, and $\lVert b\rVert_{L^{p}_{x}L^{q}_{t}(Q_{I})}$.

_(Single-scale Harnack)_ If, additionally, $b\in L^{2}_{t}H^{-1}_{x}(Q_{I})$
and $\theta>0$, then we have the following _quantitative Harnack inequality_ :
If $I_{1},I_{2}\subset\subset I$ are intervals satisfying $\sup I_{1}<\inf
I_{2}$, then

$\sup_{\Omega^{\prime}\times I_{1}}\theta\lesssim\inf_{\Omega^{\prime}\times
I_{2}}\theta,$ (1.16)

where the implied constant depends on $n$, $\Omega$, $\Omega^{\prime}$, $I$,
$I_{1}$, $I_{2}$, $p$, $q$, $\lVert b\rVert_{L^{p}_{x}L^{q}_{t}(Q_{I})}$, and
$\lVert b\rVert_{L^{2}_{t}H^{-1}_{x}(Q_{I})}$.

_(Sharpness, steady case)_. Let $n\geq 3$. The quantitative local boundedness
property fails for steady drifts $b\in L^{\frac{n-1}{2}}(B)$ and steady
solutions $\theta$ in the ball $B$.

_(Sharpness, time-dependent case)_. Let $n\geq 2$ and $Q=B_{1}\times(0,1)$.
There exist a smooth divergence-free drift $b\in C^{\infty}(Q)$ belonging to
$L^{p}_{x}L^{q}_{t}(Q)$ for all $p,q\in[1,+\infty]$ with $p\leq q$ and
$3/q+(n-1)/p>2$ and satisfying the following property. There exists a smooth
solution $\theta\in L^{\infty}_{t}L^{1}_{x}\cap C^{\infty}(Q)$ to the
advection-diffusion equation (A-D) in $Q$ with

$\sup_{B_{1/2}\times(0,T)}|\theta|\to+\infty\text{ as }T\to 1_{-}.$

In particular, the above quantitative local boundedness property fails when
$3/q+(n-1)/p>2$ and $p\leq q$. Finally, the drift additionally belongs to
$L^{q}_{t}L^{p}_{x}(Q)$ for all $(p,q)\in[1,+\infty]^{2}$ with $2/q+n/p>2$.

_(Upper bounds on fundamental solutions)_ If the divergence-free drift $b\in
C^{\infty}_{0}(\mathbb{R}^{n}\times[0,+\infty))$ belongs to
$L^{p}_{x}L^{q}_{t}(\mathbb{R}^{n}\times\mathbb{R}_{+})$ and
$1\leq\zeta_{0}<2$, then the fundamental solution $\Gamma=\Gamma(x,t;y,s)$ to
the parabolic operator $L=\partial_{t}-\Delta+b\cdot\nabla$ satisfies

$\displaystyle\Gamma(x,t;y,s)$ $\displaystyle\leq
C(t-s)\left(\frac{1}{t-s}+\frac{M_{0}}{(t-s)^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}$
(1.17)
$\displaystyle\quad\times\left[\exp\left(-\frac{c|x-y|^{2}}{t-s}\right)+\exp\left(-\frac{c|x-y|^{1+\frac{1}{1+\alpha_{0}}}}{((t-s)M_{0})^{\frac{1}{1+\alpha_{0}}}}\right)\right]$

for all $x,y\in\mathbb{R}^{n}$ and $0\leq s<t<+\infty$. Here,

$\alpha_{0}=\frac{\zeta_{0}-1}{\theta_{2}},\quad
M_{0}=C\|b\|_{L^{p}_{x}L^{q}_{t}(\mathbb{R}^{n}\times\mathbb{R}_{+})}^{\frac{1}{\theta_{2}}}+\frac{1}{4},\quad\theta_{2}=1-\frac{\zeta_{0}}{2}.$
(1.18)

See Figures 2 and 3 for an illustration of Theorem 1.2.

Figure 2. Divergence-free drift $b\in L^{p}_{x}L^{q}_{t}$, dimension $n\geq 3$
(dimension $n=4$ illustrated above). Local boundedness and single-scale
Harnack inequality. Figure 3. Divergence-free drift $b\in L^{p}_{x}L^{q}_{t}$,
dimension $n=2$. Local boundedness and single-scale Harnack inequality.

We prove Theorems 1.1 and 1.2, including the known results, in a self-
contained way below.

######  Remark 1.3.

The local boundedness property and Harnack inequality in Theorems 1.1 and 1.2
can be easily extended to accomodate drifts satisfying
$\operatorname{div}b\leq 0$ (with the background assumption $b\in
L^{2}_{t,x}(Q_{I})$ in the Harnack inequality). These properties and the
fundamental solution estimates can also be extended to divergence-form
elliptic operators $\operatorname{div}a\nabla\cdot$ with bounded, uniformly
elliptic $a$.

######  Remark 1.4.

Our condition $b\in L^{2}_{t}H^{-1}_{x}(Q_{I})$ appears to be new and is used
to connect the forward- and backward-in-time regions in the Harnack
inequality, see Lemma 2.11. In contrast, the background condition in the
single-scale Harnack inequality in [IKR16] is $b\in
L^{\infty}_{t}L^{2}_{x}(Q_{I})$.

### Discussion of dimension reduction principle

The ‘slicing’ described above in the steady setting is more subtle in the
time-dependent setting because the anisotropic condition $\theta\in
L^{\infty}_{t}L^{2}_{x}$ does not restrict well to slices in the radial
variable $r$; compare this to the isotropic condition $\nabla\theta\in
L^{2}_{t,x}$. Indeed, to ‘slice’ in a variable, it seems necessary for that
variable to be summed ‘last’ (that is, on the outside) in the norm. The
condition $b\in L^{p}_{x}L^{q}_{t}$, $p\leq q$, $\frac{3}{q}+\frac{n-1}{p}<2$,
in Theorem 1.2 comes, roughly speaking, from interpolating between the
isotropic condition $b\in L^{\frac{n+2}{2}+}_{t,x}$, in which the order of
integration may be changed freely, and the dimensionally reduced condition
$b\in L^{\frac{n-1}{2}+}_{x}L^{\infty}_{t}$, which implies that $b\in
L^{\infty}_{t}L^{\frac{n-1}{2}+}_{\sigma}(\partial B_{r}\times I)$ on ‘many
slices’, say, a set of $r\in A\subset(1/2,1)$ with measure $|A|>1/4$. Local
boundedness under this condition was already observed by X. Zhang in [Zha20,
Corollary 1.5], and the counterexamples we construct answer an open question
in Remark 1.6 therein.

Our proof of local boundedness and the Harnack inequality is built on a
certain _form boundedness condition_ (FBC), see Section 2, which subsumes a
wide variety of possible assumptions on $b$. For example, in Proposition 2.3,
we verify (FBC) not only in the context of Theorems 1.1 and 1.2 but also under
the more general conditions

$b\in L^{q}_{t}L^{\beta}_{r}L^{\gamma}_{\sigma}((B_{R}\setminus B_{R/2})\times
I),\quad\beta\geq\frac{n}{2},\quad\frac{2}{q}+\frac{1}{\beta}+\frac{n-1}{\gamma}<2$
(1.19)

and

$b\in L^{\kappa}_{r}L^{q}_{t}L^{p}_{\sigma}((B_{R}\setminus B_{R/2})\times
I),\quad p\leq q,\quad\frac{3}{q}+\frac{n-1}{p}<2.$ (1.20)

Furthermore, we allow arbitrarily low integrability $\kappa>0$ in the radial
variable; the slicing method does not require high integrability. Our proof of
upper bounds on fundamental solutions is centered on a variant
($\widetilde{\text{FBC}}$) of the form boundedness condition, see Section 5,
partially inspired by the work of Qi S. Zhang [Zha04].

We now describe the work [BS21], which was generalized to the parabolic
setting in [Zha20]. The conditions in [BS21] are on the ellipticity matrix
$a$, which is allowed to be degenerate: Define

$\lambda(x):=\inf_{|\xi|=1}\xi\cdot
a(x)\xi,\quad\mu(x):=\sup_{|\xi|=1}\frac{|a(x)\xi|^{2}}{\xi\cdot a(x)\xi}.$
(1.21)

If $n\geq 2$, $p,q\in(1,+\infty]$, and

$\lambda^{-1}\in L^{q}(B),\quad\mu\in
L^{p}(B),\quad\frac{1}{p}+\frac{1}{q}<\frac{2}{n-1},$ (1.22)

then weak solutions of $-\operatorname{div}a\nabla u=0$ are locally bounded
and satisfy a single-scale Harnack inequality. The analogous condition with
$\frac{2}{n}$ on the right-hand side is due to Trudinger in [Tru71]. By
examples in [FSSC98], the right-hand side cannot be improved to
$\frac{2}{n-1}+\varepsilon$. Divergence-free drifts $b$ belong to the above
framework: Under general conditions, it is possible to realize $b$ as the
divergence of an antisymmetric _stream matrix_ : $b_{i}=d_{ij,i}$. Then we
have
$-\Delta\theta+b\cdot\nabla\theta=-\operatorname{div}[(I+d)\nabla\theta]$, and
$\mu$ captures the antisymmetric part $d$. The steady examples we construct in
Section 4 handle the equality case in (1.22). We mention also the works
[BS20b, BS20a].

Earlier, it was hoped that the dimension reduction could be further adapted to
treat the case $b\in L^{\frac{n+1}{2}+}_{t,x}$ in the parabolic setting by
estimating a half-derivative in time: $|\partial_{t}|^{1/2}\theta\in
L^{2}_{t,x}$, since this condition is better adapted to slicing than
$\theta\in L^{\infty}_{t}L^{2}_{x}$. On the other hand, our counterexamples
rule out this possibility. Half time derivatives in parabolic PDE go back, at
least, to [LSUt68, Chapter III, Section 4], see [ABES19] for further
discussion.

### Discussion of counterexamples and ‘bounded total speed’

Solutions of (A-D) in the whole space and evolving from initial data
$\theta_{0}\in L^{1}(\mathbb{R}^{n})$ become bounded instantaneously. This is
captured by the famous Nash estimate [Nas58]:

$\|\theta(\cdot,t)\|_{L^{\infty}(\mathbb{R}^{n})}\lesssim
t^{-\frac{n}{2}}\|\theta_{0}\|_{L^{1}(\mathbb{R}^{n})},$ (1.23)

where the implied constant is _independent_ of the divergence-free drift $b$.
The Nash estimate indicates that a divergence-free drift does not impede
smoothing, in the sense of boundedness, of a density, even if the density is
initially a Dirac mass. Therefore, for rough drifts, local boundedness must be
violated in a different way: The danger is that the drift can ‘drag’ an
anomolous singularity into the domain of observation from outside. There is a
competition between the drift, which transports the singularity with some
speed, and the diffusion, which smooths the singularity at some rate. Will the
singularity, entering from outside, be smoothed before it can be observed
inside the domain?

Consider a Dirac mass $\delta_{x=-\vec{e_{1}}}$, which we seek to transport
inside the domain. If one can transport the Dirac mass inside $B_{1/2}$
_instantaneously_ , one can violate local boundedness. This can be done easily
via the drift $b(x,t)=\delta_{t=0}\vec{e}_{1}$, which is singular in time.
This example already demonstrates the importance of the space
$L^{1}_{t}L^{\infty}_{x}$, whose drifts cannot transport the mass inside
arbitrarily quickly.

To improve this example, we seek the most efficient way to transport the Dirac
mass. Heuristically, the evolution of the Dirac mass is mostly supported in a
ball of radius $R(t)\sim\sqrt{t}$. Therefore, we define our drift $b$ to be
$S(t)\vec{e}_{1}$ restricted to this support. That is, the drift lives on a
ball of radius $R(t)$ moving in the $x_{1}$-direction at speed $S(t)$. Since
we wish to move the Dirac mass instantaneously, we guess that $S(t)\sim 1/t$.
A back-of-the-envelope calculation gives

$\|b\|_{L^{q}_{t}L^{p}_{x}}^{q}\sim\int_{0}^{1}S(t)^{q}R(t)^{\frac{nq}{p}}\,dt\sim\int_{0}^{1}t^{-q+\frac{nq}{2p}}\,dt.$
(1.24)

The above quantity is finite when $2/q+n/p>2$; more care is required to get
the borderline cases in Theorem 1.1, see Section 4. This heuristic is the
basis for our time-dependent counterexamples in Section 4, except that we use
appropriate subsolutions to keep the compact support property, we glue
together many of these Dirac masses, and $S(t)$ must be chosen more carefully.

The above transport-vs.-smoothing phenomenon is also responsible for the ‘fat
tails’ in the upper bounds (1.11) and (1.17) on the fundamental solutions.
These upper bounds do not align with the Nash estimate (1.23) because the Nash
estimate does not effectively capture spatial localization.

The elliptic counterexample with $b\in L^{\frac{n-1}{2}}$ is achieved by
introducing an ansatz which reduces the problem to counterexamples for the
steady Schrödinger equation $-\Delta u+Vu=0$ in dimension $n-1$. These steady
counterexamples are singular on a line through the domain, as they must be to
respect the maximum principle.

The time-dependent counterexamples in $L^{p}_{x}L^{q}_{t}$ seem to be more
subtle, and we only exhibit them in the non-borderline cases
$\zeta_{0}=\frac{3}{q}+\frac{n-1}{p}>2$ and $p\leq q$. When $\zeta_{0}=2$, we
have counterexamples in the cases $p=q=\frac{n+2}{2}$ and
$(p,q)=(\frac{n-1}{2},+\infty)$ (the steady example). We believe that local
boundedness fails also between these two points, but the counterexamples are
yet to be exhibited, see Remark 4.4.

### Review of existing literature

Following the seminal works of De Giorgi [DG57] and Nash [Nas58], Moser
introduced his parabolic Harnack inequality [Mos64, Mos67] (see [Mos61] for
the elliptic case), whose original proof relied on a parabolic generalization
of the John-Nirenberg theorem concerning exponential integrability of ${\rm
BMO}$ functions. Later, Moser published a simplified proof [Mos71], whose
basic methods we follow. In [SSvZ12], Seregin, Silvestre, Šverák, and Zlatoš
generalized Moser’s methods to accomodate drifts in $L^{\infty}_{t}{\rm
BMO}^{-1}_{x}$. For recent work on boundary behavior in this setting, see
[LP19, HLMP21]. Generalizations to critical Morrey spaces and the
supercritical Lebesgue spaces are due to [NUt11, IKR14, IKR16, Ign14].

The Gaussian estimates on fundamental solutions were discovered by Aronson
[Aro67] and were generalized to divergence-free drifts by Osada in [Osa87]
($L^{\infty}_{t}L^{-1,\infty}_{x}$) and Qian and Xi ($L^{\infty}_{t}{\rm
BMO}^{-1}_{x}$) in [QX19b, QX19a]. Important contributions are due to [Zha04],
who developed Gaussian-like upper bounds in the supercritical case $b\in
L^{\frac{n}{2}+}(\mathbb{R}^{n})$, $n\geq 4$, and [LZ04, MS04, Sem06], among
others. For recent progress on Green’s function estimates with sharp
conditions on lower order terms, see [DK18, KS19, Sak20, Mou19].

The primary examples concerning the regularity of solutions to (A-D) can be
found in [SSvZ12, SVZ13, Wu21]. Counterexamples to continuity with time-
dependent drifts can be constructed by colliding two discs of $+1$
(subsolution) and $-1$ (supersolution) with radii $R(t)\sim\sqrt{1-t}$ and
speeds $S(t)\sim 1/\sqrt{1-t}$. The parabolic counterexamples with _steady_
velocity fields constructed therein are more challenging. See [Fil13, FS18]
for examples in the elliptic setting. We also mention Zhikov’s counterexamples
[Zhi04] to uniqueness when $b$ does not belong to $L^{2}$, whereas weak
solutions with zero Dirichlet conditions are known to be unique when $b\in
L^{2}$ [Zha04].

For recent counterexamples in the regularity theory of parabolic systems based
on self-similarity, see [Moo17].

######  Remark 1.5.

At a technical level, there is a small but, perhaps, non-trivial gap in the
proof of the weak Harnack inequality in [IKR16], see (3.22) therein, where it
is claimed that $\log_{+}(\theta/\mathbf{K})$ is a supersolution. This seems
related to a step in the proof of Lemma 6.20, p. 124, in Lieberman’s book
[Lie96], which we had difficulty following, see the first inequality therein.
Both of these are related to improving the weak $L^{1}$ inequality.

We opt to follow Moser’s proof in [Mos71] more directly and skip the weak
Harnack inequality. In principle, one could directly apply the parabolic John-
Nirenberg inequality in [Mos64, FG85] to obtain the weak Harnack inequality.

##  2\. Local boundedness and Harnack’s inequality

Let $b$ be a smooth, divergence-free vector field defined on $B_{R_{0}}\times
I_{0}$, where $R_{0}>0$ and $I_{0}$ is an open interval. In the sequel, we
will use a _form boundedness condition_ , which we denote by (FBC):

> There exist constants $M,N,\alpha>0$, $\varepsilon\in[0,1/2)$, and
> $\delta\in(0,1]$ satisfying the following property. For every
> $R\in[R_{0}/2,R_{0}]$, subinterval $I\subset I_{0}$, and Lipschitz $u\in
> W^{1,\infty}(B_{R}\times I)$, there exists a measurable set
> $A=A(\varrho,R,I,u)\subset(\varrho,R)$ with $|A|\geq\delta(R-\varrho)$ and
> satisfying
>
> $\displaystyle-\frac{1}{|A|}\int\kern-3.39996pt\int_{B_{A}\times
> I}\frac{|u|^{2}}{2}(b\cdot
> n)\,dx\,dt\leq\frac{MR_{0}^{\alpha}}{\delta^{\alpha}R_{0}^{2}(R-\varrho)^{\alpha}}\int\kern-3.39996pt\int_{(B_{R}\setminus
> B_{\varrho})\times I}|u|^{2}\,dx\,dt$ (FBC)
> $\displaystyle\quad\quad+N\int\kern-3.39996pt\int_{(B_{R}\setminus
> B_{\varrho})\times I}|\nabla u|^{2}\,dx\,dt+\varepsilon\sup_{t\in
> I}\int_{B_{R}}|u(x,t)|^{2}\,dx,$
>
> where $B_{A}=\cup_{r\in A}\partial B_{r}$ and $n$ is the outer unit normal.

The LHS of (FBC) appears on the RHS of the energy estimates.

In the situations we consider, $M$ may depend on $R_{0}$, and we can predict
its dependence based on dimensional analysis. For example, since $b$ has
dimensions of $L^{-1}$, the quantity

$R_{0}^{1-\frac{2}{q}-\frac{n}{p}}\|b\|_{L^{q}_{t}L^{p}_{x}(B_{R_{0}}\times
R_{0}^{2}I)}$

is dimensionless.

In Proposition 2.3, we show that (FBC) is satisfied under the hypotheses of
Theorems 1.1 and 1.2.

_Notation_. In this section, $R_{0}/2\leq\varrho<R\leq R_{0}$ and
$-\infty<T<\tau<0$. Let us introduce the backward parabolic cylinders
$Q_{R,T}=B_{R}\times(T,0)$. Our working assumptions are that $\theta$ is a
non-negative Lipschitz function and $b$ is a smooth, divergence-free vector
field. To give precise constants, we will frequently use the notation

$\mathbf{C}(\varrho,\tau,R,T,M,\delta,\alpha)=\frac{1}{\delta^{2}(R-r)^{2}}+\frac{MR_{0}^{\alpha}}{\delta^{\alpha}R_{0}^{2}(R-r)^{\alpha}}+\frac{1}{\tau-T}$
(2.1)

involving the various parameters from (FBC). Our convention throughout the
paper is that all implied constants may depend on $n$.

######  Theorem 2.1 (Local boundedness).

Let $\theta$ be a non-negative Lipschitz subsolution and $b$ satisfy (FBC) on
$Q_{R,T}$. Then, for all $\gamma\in(0,2]$,

$\sup_{Q_{\varrho,\tau}}\theta\lesssim_{N,\alpha,\varepsilon,\gamma}\mathbf{C}^{\frac{n+2}{2\gamma}}\lVert\theta\rVert_{L^{\gamma}(Q_{R,T}\setminus
Q_{\varrho,\tau})}.$ (2.2)

######  Theorem 2.2 (Harnack inequality).

Let $\theta$ be a non-negative Lipschitz solution on
$Q^{*}=B\times(-T^{*},T^{*})$. Let $b\in
L^{2}_{t}H^{-1}_{x}(Q^{*};\mathbb{R}^{n})$ satisfying (FBC) on $Q^{*}$. Let
$0<\ell<T^{*}$ be the time lag. Then

$\sup_{B_{1/2}\times(-T^{*}+\ell,0)}\theta\lesssim_{N,M,A,T^{*},\delta,\alpha,\varepsilon,\ell}\inf_{B_{1/2}\times(\ell,T^{*})}\theta,$
(2.3)

where $A=\lVert b\rVert_{L^{2}_{t}H^{-1}_{x}(Q^{*};\mathbb{R}^{n})}^{2}$.

###  2.1. Verifying (FBC)

We verify that (FBC) is satisfied in the setting of the main theorems:

######  Proposition 2.3 (Verifying FBC).

Let $p,q,\beta,\gamma\in[1,+\infty]$, $\kappa\in(0,+\infty]$, and $b$ be a
smooth, divergence-free vector field defined on $B_{R_{0}}\times I_{0}$.

1\. If

$b\in L^{q}_{t}L^{\beta}_{r}L^{\gamma}_{\sigma}((B_{R_{0}}\setminus
B_{R_{0}/2})\times
I_{0}),\quad\beta\geq\frac{n}{2},\quad\zeta_{0}:=\frac{2}{q}+\frac{1}{\beta}+\frac{n-1}{\gamma}<2,$
(2.4)

then $b$ satisfies (FBC) with $\alpha=1/\theta_{2}$, $A=(\varrho,R)$,
$\delta=1$, $N=\varepsilon=1/4$, and

$M=C\left(R_{0}^{1-\zeta_{0}}\|b\|_{L^{q}_{t}L^{\beta}_{r}L^{\gamma}_{\sigma}((B_{R_{0}}\setminus
B_{R_{0}/2})\times I_{0})}\right)^{\frac{1}{\theta_{2}}}+\frac{1}{4},$ (2.5)

where $\theta_{2}=1-\zeta_{0}/2$.

2\. If

$b\in L^{\kappa}_{r}L^{q}_{t}L^{p}_{\sigma}((B_{R_{0}}\setminus
B_{R_{0}/2})\times I_{0}),\quad p\leq
q,\quad\zeta_{0}:=\frac{3}{q}+\frac{n-1}{p}<2,$ (2.6)

then $b$ satisfies (FBC) with $\alpha=(1/\kappa+(q-1)/q)/\theta_{2}$,
$\delta=1/2$, $N=\varepsilon=1/4$, and

$M=C\left(R_{0}^{1-\zeta_{0}}\|b\|_{L^{\kappa}_{r}L^{q}_{t}L^{p}_{\sigma}((B_{R_{0}}\setminus
B_{R_{0}/2})\times I_{0})}\right)^{\frac{1}{\theta_{2}}}+\frac{1}{4},$ (2.7)

where $\theta_{2}=1-\zeta_{0}/2$.

######  Corollary 2.4 (FBC in $L^{p}_{x}L^{q}_{t}$).

Let $p,q\in[1,+\infty]$ and $b$ be as above. If

$b\in L^{p}_{x}L^{q}_{t}((B_{R_{0}}\setminus B_{R_{0}/2})\times I_{0}),\quad
p\leq q,\quad\zeta_{0}:=\frac{3}{q}+\frac{n-1}{p}<2,$ (2.8)

then $b$ satisfies (FBC) with $\alpha=(1/p+(q-1)/q)/\theta_{2}$, $\delta=1/2$,
$N=\varepsilon=1/4$, and

$M=C\left(R_{0}^{1-\zeta_{0}}\|b\|_{L^{p}_{x}L^{q}_{t}((B_{R_{0}}\setminus
B_{R_{0}/2})\times I_{0})}\right)^{\frac{1}{\theta_{2}}}+\frac{1}{4},$ (2.9)

where $\theta_{2}=1-\zeta_{0}/2$.

By Minkowski’s inequality, (2.8) is a special case of (2.6) with $\kappa=p$.

######  Remark 2.5.

Condition (2.4) automatically enforces $q\in(1,+\infty]$ and
$\gamma\in\left(\frac{n-1}{2},+\infty\right]$. Condition (2.6) automatically
enforces $q\in\left(\frac{n+2}{2},+\infty\right]$ and
$p\in\left(\frac{n-1}{2},+\infty\right]$.

###### Proof of Proposition 2.3.

First, we rescale $R_{0}=1$. Let $1/2\leq\varrho<R\leq 1$, $I\subset I_{0}$.
We restrict to $n\geq 3$ and summarize $n=2$ afterward. Unless stated
otherwise, the norms below are on $(B_{R}\setminus B_{\varrho})\times I$.

1\. Summary of embeddings for $u$. By the Sobolev embedding theorem, we have

$\|u\|_{L^{2}_{t}L^{2^{*}_{n}}_{x}}\lesssim\|u\|_{L^{2}}+\|\nabla
u\|_{L^{2}},\quad 2^{*}_{n}=\frac{2n}{n-2}.$ (2.10)

After interpolating with $\|u\|_{L^{\infty}_{t}L^{2}_{x}}$ and
$\|u\|_{L^{2}_{t,x}}$, we have

$\|u\|_{L^{q_{1}}_{t}L^{p_{1}}_{x}}\lesssim\|u\|_{L^{2}}^{\theta_{1}}(\|u\|_{L^{2}}+\|\nabla
u\|_{L^{2}}+\|u\|_{L^{\infty}_{t}L^{2}_{x}})^{1-\theta_{1}}$ (2.11)

for suitable $\theta_{1}=\theta_{1}(q_{1},p_{1})\in(0,1)$ whenever

$\frac{2}{q_{1}}+\frac{n}{p_{1}}>\frac{n}{2},\quad q_{1}\in[2,+\infty),\quad
p_{1}\in[2,2^{*}_{n}).$ (2.12)

Next, we employ the Sobolev embedding theorem on the spheres $\partial B_{r}$,
$r\in(\varrho,R)$:

$\|u\|_{L^{2}_{t}L^{2}_{r}L^{2^{*}_{n-1}}_{\sigma}}\lesssim\|u\|_{L^{2}}+\|\nabla
u\|_{L^{2}}.$ (2.13)

Interpolating with (2.11), we have

$\|u\|_{L^{q_{2}}_{t}L^{\beta_{2}}_{r}L^{\gamma_{2}}_{\sigma}}\lesssim\|u\|_{L^{2}}^{\theta_{2}}(\|u\|_{L^{2}}+\|\nabla
u\|_{L^{2}}+\|u\|_{L^{\infty}_{t}L^{2}_{x}})^{1-\theta_{2}}$ (2.14)

for suitable $\theta_{2}=\theta_{2}(q_{2},\beta_{2},\gamma_{2})\in(0,1)$
whenever

$\frac{2}{q_{2}}+\frac{1}{\beta_{2}}+\frac{n-1}{\gamma_{2}}>\frac{n}{2},\quad
q_{2}\in[2,+\infty),\quad\beta_{2}\in[2,2^{*}_{n}),\quad\gamma_{2}\in[2,2^{*}_{n-1}).$
(2.15)

The condition (2.15) describes a region in the parameter space
$\\{(1/q_{2},1/\beta_{2},1/\gamma_{2})\\}\subset\mathbb{R}^{3}$ in a
tetrahedron with vertices $(1/2,1/2,1/2)$, $(1/2,1/2^{*}_{n},1/2^{*}_{n})$,
$(1/2,1/2,1/2^{*}_{n-1})$, and $(0,1/2,1/2)$. We compute $\theta_{2}$
according to

$\theta_{2}\left(\frac{n+2}{2}\right)+(1-\theta_{2})\left(\frac{n+2}{2}-1\right)=\frac{2}{q_{2}}+\frac{1}{\beta_{2}}+\frac{n-1}{\gamma_{2}}.$
(2.16)

2\. Verifying (FBC) for condition (2.4). Choose $q_{2}/2=q^{\prime}$ and
$\beta_{2}/2=\beta^{\prime}$, where ′ denotes Hölder conjugates. This is
admissible according to the restrictions on $q$, $\beta$, $q_{2}$, and
$\beta_{2}$ described in (2.4), Remark 2.5, and (2.15). We choose $\gamma_{2}$
to satisfy

$\left(\frac{2}{q}+\frac{1}{\beta}+\frac{n-1}{\gamma}\right)+2\left(\frac{2}{q_{2}}+\frac{1}{\beta_{2}}+\frac{n-1}{\gamma_{2}}\right)=n+2.$
(2.17)

This is possible according to the numerology in (2.4) and (2.15). Unfolding
(2.17), we find $\gamma_{2}/2=\gamma^{\prime}$. Finally, we choose
$A=(\varrho,R)$ and compute $\theta_{2}$ according to (2.16). Then

$\displaystyle\left|\frac{1}{|A|}\int\kern-3.39996pt\int_{B_{A}\times
I}\frac{|u|^{2}}{2}b\cdot
n\,dx\,dt\right|\leq\frac{1}{R-\varrho}\|b\|_{L^{q}_{t}L^{\beta}_{r}L^{\gamma}_{\sigma}}\|u\|_{L^{q_{2}}_{t}L^{\beta_{2}}_{r}L^{\gamma_{2}}_{\sigma}}^{2}$
(2.18)
$\displaystyle\quad\overset{\eqref{eq:uinterpineq}}{\leq}\frac{C}{R-\varrho}\|b\|_{L^{q}_{t}L^{\beta}_{r}L^{\gamma}_{\sigma}}\|u\|_{L^{2}}^{2\theta_{2}}(\|u\|_{L^{2}}+\|\nabla
u\|_{L^{2}}+\|u\|_{L^{\infty}_{t}L^{2}_{x}})^{2-2\theta_{2}}$
$\displaystyle\quad\quad\leq
C\frac{\|b\|_{L^{q}_{t}L^{\beta}_{r}L^{\gamma}_{\sigma}}^{\frac{1}{\theta_{2}}}}{(R-\varrho)^{\frac{1}{\theta_{2}}}}\|u\|_{L^{2}}^{2}+\frac{1}{4}(\|u\|_{L^{2}}^{2}+\|\nabla
u\|_{L^{2}}^{2}+\|u\|_{L^{\infty}_{t}L^{2}_{x}}^{2}),$

where we used Young’s inequality in the last step. This implies (FBC).

3\. Verifying (FBC) for condition (2.6). First, we identify _good slices_ for
$b$. Specifically, we apply Chebyshev’s inequality in $r$ to the integrable
function

$r\mapsto\big{\lVert}b|_{I\times\partial
B_{r}}\big{\rVert}_{L^{q}_{t}L^{p}_{\sigma}(I\times\partial B_{r})}^{\kappa}$
(2.19)

to obtain that, on a set $A=A(\varrho,R,I)$ of measure $|A|\geq(R-\varrho)/2$,
we have

$\lVert b\rVert_{L^{\infty}_{r}L^{q}_{t}L^{p}_{\sigma}(B_{A}\times
I)}\lesssim\frac{1}{(R-\varrho)^{\frac{1}{\kappa}}}\lVert
b\rVert_{L^{\kappa}_{r}L^{q}_{t}L^{p}_{\sigma}((B_{R}\setminus
B_{\varrho})\times I)}$ (2.20)

where $B_{A}=\cup_{r\in A}\partial B_{r}$.

Now we choose $q_{2}/2=q^{\prime}$, $\beta_{2}=q_{2}$, and
$\gamma_{2}\in[2,2^{*}_{n-1})$ satisfying (2.15) and

$\left(\frac{3}{q}+\frac{n-1}{p}\right)+2\left(\frac{3}{q_{2}}+\frac{n-1}{\gamma_{2}}\right)=n+2.$
(2.21)

Unfolding (2.21), we discover $\gamma_{2}/2=p^{\prime}$. This allows us to
compute $\theta_{2}$ according to (2.16). Moreover, since $q_{2}=\beta_{2}$ in
(2.15), we have
$L^{q_{2}}_{t}L^{\beta_{2}}_{r}L^{\gamma_{2}}_{\sigma}=L^{\beta_{2}}_{r}L^{q_{2}}_{t}L^{\gamma_{2}}_{\sigma}$.
Then

$\displaystyle\left|\frac{1}{|A|}\int\kern-3.39996pt\int_{B_{A}\times
I}\frac{|u|^{2}}{2}b\cdot
n\,dx\,dt\right|\leq\frac{C}{R-\varrho}\|b\|_{L^{\infty}_{r}L^{q}_{t}L^{p}_{\sigma}}(R-\varrho)^{1-\frac{2}{\beta_{2}}}\|u\|_{L^{\beta_{2}}_{r}L^{q_{2}}_{t}L^{\gamma_{2}}_{\sigma}}^{2}$
(2.22)
$\displaystyle\quad\overset{\eqref{eq:uinterpineq}}{\leq}\frac{C}{(R-\varrho)^{\frac{1}{\kappa}+\frac{2}{\beta_{2}}}}\lVert
b\rVert_{L^{\kappa}_{r}L^{q}_{t}L^{p}_{\sigma}((B_{R}\setminus
B_{\varrho})\times I)}\|u\|_{L^{2}}^{2\theta_{2}}(\|u\|_{L^{2}}+\|\nabla
u\|_{L^{2}}+\|u\|_{L^{\infty}_{t}L^{2}_{x}})^{2(1-\theta_{2})}$
$\displaystyle\quad\quad\leq C\frac{\lVert
b\rVert_{L^{\kappa}_{r}L^{q}_{t}L^{p}_{\sigma}((B_{1}\setminus B_{1/2})\times
I)}^{\frac{1}{\theta_{2}}}}{(R-\varrho)^{\frac{1}{\theta_{2}}(\frac{1}{\kappa}+\frac{2}{\beta_{2}})}}\|u\|_{L^{2}}^{2}+\frac{1}{4}(\|u\|_{L^{2}}^{2}+\|\nabla
u\|_{L^{2}}^{2}+\|u\|_{L^{\infty}_{t}L^{2}_{x}}^{2}).$

This completes the proof of (FBC).

4\. Dimension $n=2$. The Sobolev embedding (2.13) on the sphere bounds instead

$\|u\|_{L^{2}_{t}L^{2}_{r}C^{1/2}_{\sigma}}\lesssim\|u\|_{L^{2}}+\|\nabla
u\|_{L^{2}},$

so we must adjust the proof of the interpolation (2.14). After the initial
interpolation step between $L^{2}_{t}L^{2^{*}_{n}}_{x}$,
$L^{\infty}_{t}L^{2}_{x}$ and $L^{2}_{t,x}$ (where $2^{*}_{n}$ is now any
large finite number), we apply the following Gagliardo-Nirenberg inequality on
spheres:

$\|u\|_{L^{q}(B_{r})}\lesssim\|u\|_{L^{p}(B_{r})}^{\theta}\|u\|_{C^{1/2}(B_{r})}^{1-\theta},\quad-\frac{n-1}{q}=\theta\left(-\frac{n-1}{p}\right)+\frac{1-\theta}{2},$
(2.23)

where $p,q\in[1,+\infty]$, $\theta\in(0,1]$, and $r\in[1/2,1]$. This allows us
to recover (2.14), where now $\gamma_{2}=+\infty$ is also allowed, and
complete the proof.

To prove (2.23), we first use local coordinates on the sphere and a partition
of unity333Alternatively, since $n=2$, we could argue on the flat torus
without a partition of unity. The argument we present here is more general. to
reduce to functions $f$ on $\mathbb{R}^{n-1}$. Next, we use that $L^{p}\subset
B^{-\frac{n-1}{p}}_{\infty,\infty}$, $C^{1/2}=B^{1/2}_{\infty,\infty}$, and
real interpolation

$[B^{-\frac{n-1}{p}}_{\infty,\infty},B^{1/2}_{\infty,\infty}]_{\theta,1}=B^{0}_{q,1}\subset
L^{q}$ (2.24)

to demonstrate

$\|f\|_{L^{q}}\lesssim\|f\|_{L^{p}(B_{r})}^{\theta}\|f\|_{C^{1/2}(B_{r})}^{1-\theta}\lesssim\varepsilon^{\frac{1}{\theta}}\|f\|_{L^{p}(B_{r})}+\varepsilon^{-\frac{1}{1-\theta}}\|f\|_{C^{1/2}(B_{r})}.$
(2.25)

We piece together $u$ from the functions $f$ and optimize in $\varepsilon$ to
obtain (2.23). ∎

######  Remark 2.6 (On $\widetilde{\text{FBC}}$).

For pointwise upper bounds on fundamental solutions in Section 5, we require a
global, revised form boundedness condition ($\widetilde{\text{FBC}}$), in
which we allow $\varepsilon$ to vary in $(0,+\infty)$. We refine (2.18) by
keeping track of the $\varepsilon$ dependence in Hardy’s inequality:

$\displaystyle\left|\frac{1}{|A|}\int\kern-3.39996pt\int_{B_{A}\times
I}\frac{|u|^{2}}{2}b\cdot
n\,dx\,dt\right|\leq\frac{C\|b\|_{L^{q}_{t}L^{\beta}_{r}L^{\gamma}_{\sigma}((B_{R}\setminus
B_{\varrho})\times
I)}^{\frac{1}{\theta_{2}}}}{\varepsilon^{\frac{1-\theta_{2}}{\theta_{2}}}(R-\varrho)^{\frac{1}{\theta_{2}}}}\|u\|_{L^{2}}^{2}+\varepsilon\left(\|u\|_{L^{2}}^{2}+\|\nabla
u\|_{L^{2}}^{2}+\|u\|_{L^{\infty}_{t}L^{2}_{x}}^{2}\right),$ (2.26)

A similar refinement holds in (2.22). We require these refinements in the
justification of Lemma 5.3.

###  2.2. Proof of local boundedness

To begin, we prove Cacciopoli’s inequality:

######  Lemma 2.7 (Cacciopoli inequality).

Under the hypotheses of Theorem 2.1,

$\sup_{t\in(\tau,0)}\int_{B_{r}}|\theta(x,t)|^{2}\,dx+\int\kern-3.39996pt\int_{Q_{\varrho,\tau}}|\nabla\theta|^{2}\,dx\,dt\lesssim_{N,\alpha,\varepsilon}\mathbf{C}\int\kern-3.39996pt\int_{Q_{R,T}\setminus
Q_{\varrho,\tau}}|\theta|^{2}\,dx\,dt.$ (2.27)

###### Proof of Lemma 2.7.

Let $\eta\in C^{\infty}_{0}(T,+\infty)$ satisfying $0\leq\eta\leq 1$ on
$\mathbb{R}$, $\eta\equiv 1$ on $(\tau,+\infty)$ and $0\leq d\eta/dt\lesssim
1/(\tau-T)$. Let $r\in(\varrho,R)$ and $t\in(\tau,0)$. To begin, we multiply
by $\theta\eta^{2}$ and integrate over $B_{r}\times(T,t)$:

$\displaystyle\frac{1}{2}\int_{B_{r}}|\theta(x,t)|^{2}\,dx+\int\kern-3.39996pt\int_{B_{r}\times(\tau,t)}|\nabla\theta|^{2}\,dx\,ds$
(2.28)
$\displaystyle\quad\leq\int\kern-3.39996pt\int_{B_{R}\times(T,t)}|\theta|^{2}\frac{d\eta}{ds}\eta\,dx\,ds+\int\kern-3.39996pt\int_{\partial
B_{r}\times(T,t)}\Big{(}\frac{d\theta}{dn}\theta\eta^{2}-\frac{\theta^{2}}{2}(b\cdot
n)\eta^{2}\Big{)}\,d\sigma\,ds.$

Next, we average in the $r$ variable over the set of ‘good slices’,
$A=A(\varrho,R,(T,t),\theta\eta)$, which was defined in (FBC):

$\displaystyle\frac{1}{2}\int_{B_{\varrho}}|\theta(x,t)|^{2}\,dx+\int\kern-3.39996pt\int_{B_{\varrho}\times(\tau,t)}|\nabla\theta|^{2}\,dx\,ds$
(2.29)
$\displaystyle\quad\leq\frac{C}{\tau-T}\int\kern-3.39996pt\int_{B_{R}\times(T,\tau)}|\theta|^{2}\,dx\,ds+\frac{1}{|A|}\int\kern-3.39996pt\int_{B_{A}\times(T,t)}\Big{(}\frac{d\theta}{dn}\theta\eta^{2}-\frac{\theta^{2}}{2}(b\cdot
n)\eta^{2}\Big{)}\,dx\,ds,$

where $B_{A}=\cup_{r\in A}\partial B_{r}$. Let us estimate the term containing
$d\theta/dn$:

$\displaystyle\frac{1}{|A|}\int\kern-3.39996pt\int_{B_{A}\times(T,t)}\frac{d\theta}{dn}\theta\eta^{2}\,dx\,ds$
(2.30)
$\displaystyle\quad\leq\frac{1}{\delta^{2}(R-r)^{2}}\int\kern-3.39996pt\int_{Q_{R,T}\setminus
Q_{\varrho,\tau}}|\theta|^{2}\,dx\,ds+\int\kern-3.39996pt\int_{Q_{R,T}\setminus
Q_{\varrho,\tau}}|\nabla\theta|^{2}\,dx\,ds.$

To estimate the term containing $b$, we use (FBC) with $u=\theta\eta$:

$\displaystyle-\frac{1}{|A|}\int\kern-3.39996pt\int_{B_{A}\times(T,t)}\frac{\theta^{2}}{2}(b\cdot
n)\eta^{2}\,dx\,ds\leq\frac{M}{\delta^{\alpha}(R-r)^{\alpha}}\int\kern-3.39996pt\int_{Q_{R,T}\setminus
Q_{\varrho,\tau}}|\theta|^{2}\,dx\,ds$ (2.31)
$\displaystyle\quad+N\int\kern-3.39996pt\int_{Q_{R,T}\setminus
Q_{\varrho,\tau}}|\nabla\theta|^{2}\,dx\,ds+\varepsilon\sup_{s\in(T,t)}\int_{B_{R}}|\theta(x,s)|^{2}\,dx.$

Combining everything and applying $\sup_{t\in(T,0)}$, we obtain

$\displaystyle\frac{1}{2}\sup_{t\in(\tau,0)}\int_{B_{\varrho}}|\theta(x,t)|^{2}\,dx+\int\kern-3.39996pt\int_{Q_{\varrho,\tau}}|\nabla\theta|^{2}\,dx\,dt\leq
C\times\mathbf{C}\int\kern-3.39996pt\int_{Q_{R,T}\setminus
Q_{\varrho,\tau}}|\theta|^{2}\,dx\,dt$ (2.32)
$\displaystyle\quad\quad+(1+N)\int\kern-3.39996pt\int_{Q_{R,T}\setminus
Q_{\varrho,\tau}}|\nabla\theta|^{2}\,dx\,dt+\varepsilon\sup_{t\in(T,0)}\int_{B_{R}}|\theta(x,t)|^{2}\,dx.$

By Widman’s hole-filling trick, there exists
$\gamma:=\max\\{(N+1)/(N+2),2\varepsilon\\}\in(0,1)$ satisfying

$\displaystyle\frac{1}{2(N+2)}\sup_{t\in(\tau,0)}\int_{B_{\varrho}}|\theta(x,t)|^{2}\,dx+\int\kern-3.39996pt\int_{Q_{\varrho,\tau}}|\nabla\theta|^{2}\,dx\,dt\leq
C(N)\times\mathbf{C}\int\kern-3.39996pt\int_{Q_{R,T}\setminus
Q_{\varrho,\tau}}|\theta|^{2}\,dx\,dt$ (2.33)
$\displaystyle\quad\quad+\gamma\int\kern-3.39996pt\int_{Q_{R,T}}|\nabla\theta|^{2}\,dx\,dt+\frac{\gamma}{2(N+2)}\sup_{t\in(T,0)}\int_{B_{R}}|\theta(x,t)|^{2}\,dx.$

To remove the extra terms on the RHS, we use a standard iteration argument on
a sequence of scales (progressing ‘outward’) $\varrho_{0}=\varrho$,
$\varrho_{k+1}=\varrho+(1-\lambda^{k+1})(R-\varrho)$, $R_{k}=\varrho_{k+1}$,
$\tau_{0}=\tau$, $\tau_{k+1}=\tau+(1-\lambda^{2k+2})(T-\tau)$,
$T_{k}=\tau_{k+1}$, $k=0,1,2,\ldots$, where $0<\lambda<1$ is defined by the
relation $\lambda^{\max(\alpha,2)}=\gamma/2$. See [Giu03, p. 191, Lemma 6.1],
for example. This gives the desired Cacciopoli inequality. ∎

Next, we require a simple corollary:

######  Corollary 2.8 (Interpolation inequality).

Let $\chi=1+2/n$. Then

$\left(\int\kern-3.39996pt\int_{Q_{\varrho,\tau}}|\theta|^{2\chi}\,dx\,dt\right)^{\frac{1}{\chi}}\lesssim_{N,\alpha,\varepsilon}\mathbf{C}\int\kern-3.39996pt\int_{Q_{R,T}\setminus
Q_{\varrho,\tau}}|\theta|^{2}\,dx\,dt.$ (2.34)

###### Proof.

Let $0\leq\varphi\in C^{\infty}_{0}(B_{(R+\varrho)/2})$ satisfying
$\varphi\equiv 1$ on $B_{\varrho}$ and $|\nabla\varphi|\lesssim
1/(R-\varrho)$. Using (2.27) at an intermediate scale, we find

$\displaystyle\sup_{t\in(\tau,0)}\int_{B_{\varrho}}|\theta(x,t)\varphi|^{2}\,dx+\int\kern-3.39996pt\int_{Q_{\varrho,\tau}}|\nabla(\theta\varphi)|^{2}\,dx\,dt\lesssim_{N,\alpha,\varepsilon}\mathbf{C}\int\kern-3.39996pt\int_{Q_{R,T}\setminus
Q_{\varrho,\tau}}|\theta|^{2}\,dx\,dt.$ (2.35)

Then (2.34) follows from the Gagliardo-Nirenberg inequality on the whole
space. ∎

We are now ready to use Moser’s iteration:

###### Proof of local boundedness.

Let $\beta_{k}:=\chi^{k}$, where $k=0,1,2,\ldots$. A standard computation
implies that $\theta^{\beta_{k}}$ is also a non-negative Lipschitz
subsolution. Hence, it satisfies the Cacciopoli inequality (2.34) with
$R_{k}=\varrho+2^{-k}(R-\varrho)$, $r_{k}=R_{k+1}$,
$T_{k}=\tau-2^{-2k}(\tau-T)$, $\tau_{k}=T_{k+1}$, $k=0,1,2,\ldots$ (iterating
inward). In other words,

$\lVert\theta^{2\beta_{k+1}}\rVert_{L^{1}(Q_{\varrho_{k},\tau_{k}})}^{\frac{1}{\chi}}\lesssim_{N,\alpha,\varepsilon}\mathbf{C}(\varrho_{k},\tau_{k},R_{k},T_{k},M,\delta,\alpha)\lVert\theta^{2\beta_{k}}\rVert_{L^{1}(Q_{R_{k},T_{k}}\setminus
Q_{\varrho_{k},\tau_{k}})}.$ (2.36)

We may expand the domain of integration on the RHS as necessary. Define

$M_{0}:=\lVert\theta^{2}\rVert_{L^{1}(Q_{R_{0},T_{0}}\setminus
Q_{\varrho_{0},\tau_{0}})}$ (2.37)

and

$M_{k+1}:=\lVert\theta^{2}\rVert_{L^{\beta_{k+1}}(Q_{\varrho_{k},\tau_{k}})}=\lVert\theta^{2\beta_{k+1}}\rVert_{L^{1}(Q_{\varrho_{k},\tau_{k}})}^{\frac{1}{\beta_{k+1}}},\quad
k=0,1,2,\ldots$ (2.38)

Raising (2.36) to $1/\beta_{k}$ and using Eq. (2.1) defining $\mathbf{C}$, we
obtain

$M_{k+1}\leq
C(N,\alpha,\varepsilon)^{\frac{1}{\beta_{k}}}2^{\frac{\max(2,\alpha)k}{\beta_{k}}}\mathbf{C}(\varrho,\tau,R,T,M,\delta,\alpha)^{\frac{1}{\beta_{k}}}M_{k}.$
(2.39)

Iterating, we have

$M_{k+1}\leq
C(N,\alpha,\varepsilon)^{\sum_{j=0}^{k}\frac{1}{\chi^{j}}}2^{\sum_{j=0}^{k}\frac{\max(2,\alpha)j}{\chi^{j}}}\mathbf{C}(\varrho,\tau,R,T,M,\delta,\alpha)^{\sum_{j=0}^{k}\frac{1}{\chi_{j}}}M_{0}.$
(2.40)

Finally, we send $k\to+\infty$ and substitute $\sum_{j\geq
0}1/\chi^{j}=(n+2)/2$ to obtain

$\lVert\theta\rVert_{L^{\infty}(Q_{\varrho,\tau})}\lesssim_{N,\alpha,\varepsilon}\mathbf{C}^{\frac{n+2}{4}}\lVert\theta\rVert_{L^{2}(Q_{R,T}\setminus
Q_{\varrho,\tau})}.$ (2.41)

We now demonstrate how to replace $L^{2}$ on the RHS of (2.34) with
$L^{\gamma}$ ($0<\gamma<2$). To begin, use the interpolation inequality
$\lVert\theta\rVert_{L^{2}}\leq\lVert\theta\rVert_{L^{\gamma}}^{\gamma/2}\lVert\theta\rVert_{L^{\infty}}^{1-\gamma/2}$
in (2.41) and split the product using Young’s inequality. This gives

$\lVert\theta\rVert_{L^{\infty}(Q_{\varrho,\tau})}\leq
C(N,\alpha,\varepsilon,\gamma)\mathbf{C}^{\frac{n+2}{2\gamma}}\lVert\theta\rVert_{L^{\gamma}(Q_{R,T}\setminus
Q_{\varrho,\tau})}+\frac{1}{2}\lVert\theta\rVert_{L^{\infty}(Q_{R,T})}.$
(2.42)

The second term on the RHS is removed by iterating outward along a sequence of
scales, as in the proof of the Cacciopoli inequality in Lemma 2.7. ∎

######  Remark 2.9 (Elliptic case).

The analogous elliptic result is

$\sup_{B_{\varrho}}\theta\lesssim_{N,\alpha,\varepsilon,\gamma}\mathbf{C}^{\frac{n}{2\gamma}}\lVert\theta\rVert_{L^{\gamma}(B_{R}\setminus
B_{\varrho})},$ (2.43)

where
$\mathbf{C}(\varrho,R,M,\delta,\alpha)=1/[\delta^{2}(R-\varrho)^{2}]+MR_{0}^{\alpha-2}/[\delta^{\alpha}(R-\varrho)^{\alpha}]$.
The proof is the same except that $\chi=n/(n-2)$ and $\sum 1/\chi^{j}=n/2$.

###  2.3. Proof of Harnack inequality

In this subsection, $\theta$ is a strictly positive Lipschitz
solution.444There is no loss of generality if we replace $\theta$ by
$\theta+\kappa$ and let $\kappa\to 0^{+}$. Then $\log\theta$ is well defined.
Let $0\leq\psi\in C^{\infty}_{0}(B)$ be a radially decreasing function
satisfying $\psi\equiv 1$ on $B_{3/4}$. We use the notation

$f_{\rm avg}=\frac{1}{{\rm Vol}}\int_{\mathbb{R}^{n}}f\psi^{2}\,dx,\quad{\rm
Vol}=\int_{\mathbb{R}^{n}}\psi^{2}\,dx,$ (2.44)

whenever $f\in L^{1}_{\rm loc}(B)$. Let

$\mathbf{K}=\exp\left(\log\theta(\cdot,0)\right)_{\rm avg}.$ (2.45)

whose importance will be made clear in the proof of Lemma 2.11. Define

$v=\log\left(\frac{\theta}{\mathbf{K}}\right).$ (2.46)

Then $v(\cdot,0)_{\rm avg}=0$. A simple computation yields

$|\nabla v|^{2}=\partial_{t}v-\Delta v+b\cdot\nabla v.$ (2.47)

That is, $v$ is itself a supersolution, though it may not itself be positive.
We crucially exploit that $|\nabla v|^{2}$ appears on the LHS of (2.47).
First, we require the following decomposition of the drift:

######  Lemma 2.10 (Decomposition of drift).

We have the following decomposition on $B_{3/4}\times(-T^{*},T^{*})$:

$b=b_{1}+b_{2},\quad
b_{1}=-\operatorname{div}a,\quad\operatorname{div}b_{2}=0,$ (2.48)

where $a:B_{3/4}\times(-T^{*},T^{*})\to\mathbb{R}^{n\times n}_{\rm anti}$ is
antisymmetric and

$\|a\|_{L^{2}(B_{3/4}\times(-T^{*},T^{*}))}+\|b_{2}\|_{L^{2}(B_{3/4}\times(-T^{*},T^{*}))}\lesssim\|b\|_{L^{2}_{t}H^{-1}_{x}(Q^{*})},$
(2.49)

where $Q^{*}=B_{1}\times(-T^{*},T^{*})$.

###### Proof.

Let $\phi\in C^{\infty}_{0}(B_{1})$ with $\phi\equiv 1$ on $B_{15/16}$. Let
$\widetilde{b}=\phi b$. Hence,
$\|\widetilde{b}\|_{L^{2}_{t}H^{-1}_{x}(Q^{*})}\lesssim\|b\|_{L^{2}_{t}H^{-1}_{x}(Q^{*})}$.
We may decompose $\widetilde{b}(\cdot,t)\in H^{-1}(\mathbb{R}^{n})$ into
$\widetilde{b}_{1}(\cdot,t)\in\dot{H}^{-1}(\mathbb{R}^{n})$, whose Fourier
transform is supported outside of $B_{2}$, and $\widetilde{b}_{2}(\cdot,t)\in
L^{2}(\mathbb{R}^{n})$. Define

$a_{ij}=\Delta^{-1}(-\partial_{j}\widetilde{b}_{1i}+\partial_{i}\widetilde{b}_{1j}),\quad
g=\Delta^{-1}(-\operatorname{div}\widetilde{b}_{1}).$ (2.50)

This amounts to performing the Hodge decomposition in $\mathbb{R}^{n}$ ‘by
hand’.555We are simply exploiting the identity $\Delta=dd^{*}+d^{*}d$ on
differential $k$-forms, up to a sign convention, for differential $1$-forms
$\cong$ vector fields. Clearly, $a$ is antisymmetric, and we have the
decomposition

$-\widetilde{b}_{1}=\operatorname{div}a+\nabla g$ (2.51)

and the estimates

$\|a(\cdot,t)\|_{L^{2}(\mathbb{R}^{n})}+\|g(\cdot,t)\|_{L^{2}(\mathbb{R}^{n})}\lesssim\|\widetilde{b}_{1}(\cdot,t)\|_{\dot{H}^{-1}(\mathbb{R}^{n})}.$
(2.52)

Similarly, we decompose

$\widetilde{b}_{2}=\mathbb{P}\widetilde{b}_{2}+\mathbb{Q}\widetilde{b}_{2},$
(2.53)

where $\mathbb{P}$ is the Leray (orthogonal) projector onto divergence-free
fields, and $\mathbb{Q}=I-\mathbb{P}$ is the orthogonal projector onto
gradient fields. We denote $\mathbb{Q}\widetilde{b}_{2}=\nabla f$.

Since
$\widetilde{b}=\widetilde{b}_{1}+\widetilde{b}_{2}=-\operatorname{div}a+\mathbb{P}\widetilde{b}_{2}+\nabla(f-g)$
is divergence free in $B_{7/8}$ on time slices, we have

$\Delta(f-g)(\cdot,t)=0\text{ in }B_{7/8},$ (2.54)

and by elliptic regularity, for all $k\geq 0$,

$\|\nabla(f-g)(\cdot,t)\|_{H^{k}(B_{3/4})}\lesssim_{k}\|\widetilde{b}_{1}(\cdot,t)\|_{\dot{H}^{-1}(B_{7/8})}+\|\widetilde{b}_{2}(\cdot,t)\|_{L^{2}(B_{7/8})}.$
(2.55)

Finally, we define

$b_{2}=\mathbb{P}\widetilde{b}_{2}+\nabla(f-g)\in
L^{2}_{t,x}(B_{3/4}\times(-T^{*},T^{*})),$ (2.56)

which satisfies the claimed estimates and is divergence free in
$B_{7/8}\times(-T^{*},T^{*})$. ∎

We now proceed with the proof of Harnack’s inequality:

######  Lemma 2.11.

For all non-zero $t\in[-T^{*},T^{*}]$, we write $I_{t}=[0,t]$ if $t>0$ and
$I_{t}=[t,0]$ if $t<0$. Then

$-\operatorname{sgn}(t)v(\cdot,t)_{\rm avg}+\int_{I_{t}}\left(|\nabla
v|^{2}(\cdot,s)\right)_{\rm avg}\,ds\lesssim|t|+\lVert
b\rVert_{L^{2}_{t}H^{-1}_{x}(B\times I_{t})}^{2}.$ (2.57)

###### Proof.

We multiply (2.47) by $\psi^{2}$ and integrate over $B\times I_{t}$:

$\operatorname{sgn}(t)\int_{B}\left(v(x,0)-v(x,t)\right)\psi^{2}\,dx+\int\kern-3.39996pt\int_{B\times
I_{t}}|\nabla v|^{2}\psi^{2}\,dx\,ds\leq\int\kern-3.39996pt\int_{B\times
I_{t}}2\psi\nabla\psi\cdot\nabla v+(b\cdot\nabla v)\psi^{2}\,dx\,ds.$ (2.58)

By (2.45), $\int_{\mathbb{R}^{n}}v(x,0)\psi^{2}\,dx=0$. The first term on the
RHS is easily estimated:

$\int\kern-3.39996pt\int_{B\times I_{t}}2\psi\nabla\psi\cdot\nabla
v\,dx\,ds\leq\frac{1}{4}\int\kern-3.39996pt\int_{B\times I_{t}}|\nabla
v|^{2}\psi^{2}\,dx\,ds+C|t|.$ (2.59)

To estimate the term containing $b$, we require the drift decomposition
$b=b_{1}+b_{2}$ from Lemma 2.10. Then

$\int\kern-3.39996pt\int_{B\times I_{t}}(b_{1}\cdot\nabla
v)\psi^{2}\,dx\,ds=\int\kern-3.39996pt\int_{B\times I_{t}}2\psi
a(\nabla\psi,\nabla v)\,dx\,ds\leq\frac{1}{4}\int\kern-3.39996pt\int_{B\times
I_{t}}|\nabla v|^{2}\psi^{2}\,dx\,ds+C\lVert a\rVert_{L^{2}(B\times
I_{t})}^{2}.$ (2.60)

and

$\int\kern-3.39996pt\int_{B\times I_{t}}(b_{2}\cdot\nabla
v)\psi^{2}\,dx\,ds\leq\frac{1}{4}\int\kern-3.39996pt\int_{B\times
I_{t}}|\nabla v|^{2}\psi^{2}\,dx\,ds+C\lVert b_{2}\rVert_{L^{2}(B\times
I_{t})}^{2}.$ (2.61)

Recall the estimate (2.49) from the decomposition. Combining (2.58)–(2.61) and
dividing by ${\rm Vol}$ gives (2.57). ∎

In the following, we write $v=v_{+}-v_{-}$, where $v_{+},v_{-}\geq 0$. We also
use the notation

$A^{+}=\int_{0}^{T^{*}}\|b(\cdot,t)\|_{H^{-1}(B)}^{2}\,dt,\quad
A^{-}=\int_{-T^{*}}^{0}\|b(\cdot,t)\|_{H^{-1}(B)}^{2}\,dt.$ (2.62)

######  Lemma 2.12 (Weak–$L^{1}$ estimates).

With the above notation, we have

$\lVert v_{+}\rVert_{L^{1,\infty}_{t,x}(B_{3/4}\times(-T^{*},0))}\lesssim
1+T^{*}(T^{*}+A^{-})$ (2.63)

and

$\lVert v_{-}\rVert_{L^{1,\infty}_{t,x}(B_{3/4}\times(0,T^{*}))}\lesssim
1+T^{*}(T^{*}+A^{+}).$ (2.64)

###### Proof.

By (2.57) and a weighted Poincaré inequality [Mos64, Lemma 3, p. 120],

$-\operatorname{sgn}(t)v_{\rm
avg}(t)+\frac{1}{C_{1}}\int_{I_{t}}\left(|v-v_{\rm avg}|^{2}\right)_{\rm
avg}\,ds\leq
C_{0}\int_{I_{t}}\left(1+\|b(\cdot,t)\|_{H^{-1}(B)}^{2}\right)ds,$ (2.65)

where $C_{0}>0$ is the implied constant in (2.57). In the following, we focus
on the case $t\in[-T^{*},0]$. We use (2.65) to obtain a sub/supersolution
inequality corresponding to a quadratic ODE. First, we remove the forcing in
the ODE by defining

$p(x,t):=v(x,t)-C_{0}\underbrace{\int_{I_{t}}\left(1+\|b(\cdot,t)\|_{H^{-1}(B)}^{2}\right)ds}_{\leq
T^{*}+A^{-}}.$ (2.66)

Then (2.65) becomes

$p_{\rm avg}(t)+\frac{1}{C_{1}}\int_{I_{t}}\left(|p-p_{\rm
avg}|^{2}\right)_{\rm avg}\,ds\leq 0.$ (2.67)

Let us introduce the super-level sets, whose measures $\eta$ appear as a
coefficient in the ODE:

$\eta(\mu,t):=\left|\\{x\in B_{3/4}:p(x,t)>\mu\\}\right|,\quad\mu>0.$ (2.68)

Since $p_{\rm avg}\leq 0$, we have that $p(x,t)-p_{\rm avg}(t)>\mu-p_{\rm
avg}(t)>0$ whenever $p(x,t)>\mu$. Then

$p_{\rm avg}(t)+\frac{1}{C_{1}{\rm Vol}}\int_{I_{t}}\eta(\mu,s)(\mu-p_{\rm
avg})^{2}\,ds\leq 0.$ (2.69)

It is convenient to rephrase (2.69) in terms of a positive function evolving
forward-in-time: $\bar{p}(t)=-p_{\rm avg}(-t)$ with $t\in[0,T^{*}]$. Then
(2.69) becomes

$\bar{p}(t)\geq\frac{1}{C_{1}{\rm
Vol}}\int_{0}^{t}\eta(\mu,|s|)(\mu+\bar{p}(s))^{2}\,ds.$ (2.70)

The above inequality means that $\bar{p}$ is a supersolution of the quadratic
ODE

$\dot{q}=\frac{1}{C_{1}{\rm Vol}}\times\eta(\mu,|t|)(\mu+q)^{2}$ (2.71)

with $q(0)=0$. The above scalar ODE has a comparison principle. _A priori_ ,
since (2.71) is quadratic, its solutions may quickly blow-up depending on the
size of $\eta(\mu,\cdot)$ and $\mu$. However, because $\bar{p}$ lies above the
solution $q$, $q$ does not blow up, and we obtain a bound for the density
$\eta(\mu,\cdot)$ in the following way. After separating variables in (2.71),
we obtain

$\frac{1}{C_{1}{\rm
Vol}}\int_{0}^{T^{*}}\eta(\mu,|s|)\,ds=\frac{1}{\mu}-\frac{1}{\mu+q(T^{*})}\leq\frac{1}{\mu},$
(2.72)

since $q\geq 0$. That is,

$\lVert p_{+}\rVert_{L^{1,\infty}(B_{3/4}\times(-T^{*},0))}\lesssim 1.$ (2.73)

Finally, since $\lVert\cdot\rVert_{L^{1,\infty}(B_{3/4}\times(-T^{*},0))}$ is
a quasi-norm and $v\leq p+C_{0}(T^{*}+A^{-})$ pointwise due to (2.66), we have

$\displaystyle\lVert v_{+}\rVert_{L^{1,\infty}(B_{3/4}\times(-T^{*},0))}$
$\displaystyle\leq 2\lVert
p_{+}\rVert_{L^{1,\infty}(B_{3/4}\times(-T^{*},0))}+2C_{0}\lVert
T^{*}+A^{-}\rVert_{L^{1,\infty}(B_{3/4}\times(-T^{*},0))}$ (2.74)
$\displaystyle\lesssim 1+T^{*}(T^{*}+A^{-}).$

The proof for $t\in[0,T^{*}]$ is similar except that one uses sub-level sets
in (2.68) with $\mu<0$. ∎

We now require the following lemma of Moser [Mos71], which we quote almost
directly, and in which we denote by $Q(\varrho)$, $\varrho>0$ any family of
domains satisfying $Q(\varrho)\subset Q(r)$ for $0<\varrho<r$.

######  Lemma 2.13 (Lemma 3 in [Mos71]).

Let $m,\zeta,c_{0}$, $1/2\leq\theta_{0}<1$ be positive constants, and let
$w>0$ be a continuous function defined in a neighborhood of $Q(1)$ for which

$\sup_{Q(\varrho)}w^{p}<\frac{c_{0}}{(r-\varrho)^{m}{\rm
meas}(Q(1))}\int\kern-3.39996pt\int_{Q(r)}w^{p}\,dt\,dx$ (2.75)

for all $\varrho,r,p$ satisfying

$\frac{1}{2}\leq\theta_{0}\leq\varrho<r\leq 1,\quad 0<p<\zeta^{-1}.$ (2.76)

Moreover, let

${\rm meas}\\{(x,t)\in Q(1):\log w>\mu\\}<\frac{c_{0}\zeta}{s}{\rm
meas}(Q(1))$ (2.77)

for all $\mu>0$. Then there exists a constant function
$\gamma=\gamma(\theta_{0},m,c_{0})$ such that

$\sup_{Q(\theta_{0})}w<\gamma^{\zeta}.$ (2.78)

###### Proof of Harnack inequality.

We apply Lemma 2.13 to $w=\theta/\mathbf{K}$ with
$Q(\varrho)=B_{\varrho}\times(-T^{*}+2\ell(1-\varrho),0)$ and
$\theta_{0}=1/2$.666Technically, to satisfy the conditions in Lemma 2.13,
$w=\theta/\mathbf{K}$ should be extended arbitrarily to be continuous in a
neighborhood of $Q(1)$. Indeed, we recognize the requirement (2.75) as the
local boundedness guaranteed by Theorem 2.1 and (2.77) as the weak $L^{1}$
estimate from Lemma 2.12. This gives

$\sup_{B_{1/2}\times(-T^{*}+\ell,0)}\frac{\theta}{\mathbf{K}}\lesssim 1.$
(2.79)

Here, we suppress also the dependence on the time lag $\ell$. Meanwhile,
$v_{-}=\log_{+}(\mathbf{K}/\theta)$ is a subsolution. Hence,

$\lVert v_{-}\rVert_{L^{\infty}(B_{1/2}\times(\ell,T^{*}))}\lesssim\lVert
v_{-}\rVert_{L^{1,\infty}(B_{3/4}\times(0,T^{*}))}.$ (2.80)

On the other hand,

$\frac{\mathbf{K}}{\inf\theta}=\sup\frac{\mathbf{K}}{\theta}=\exp\left(\sup\log\frac{\mathbf{K}}{\theta}\right)\leq\exp\left(\sup
v_{-}\right)\lesssim\exp\left(\lVert v_{-}\rVert_{L^{1,\infty}}\right)\lesssim
1,$ (2.81)

where the $\inf$ and $\sup$ are taken on $B_{1/2}\times(\ell,T^{*})$.
Combining (2.79) and (2.81), we arrive at

$\sup_{B_{1/2}\times(-T^{*}+\ell,0)}\theta\lesssim\mathbf{K}\lesssim\inf_{B_{1/2}\times(\ell,T^{*})}\theta,$
(2.82)

as desired. ∎

##  3\. Bounded total speed

In this section, we prove the statements in Theorem 1.1 concerning the space
$L^{1}_{t}L^{\infty}_{x}$.

######  Proposition 3.1 (Local boundedness).

Let $T\in(-\infty,0)$ and $\tau\in(T,0)$. Let $b\colon
Q_{1,T}\to\mathbb{R}^{n}$ be a smooth divergence-free drift satisfying

$\lVert b\rVert_{L^{1}_{t}L^{\infty}_{x}(Q_{1,T})}\leq 1/8.$ (3.1)

Let $\theta$ be a non-negative Lipschitz subsolution on $Q_{1,T}$. Then, for
all $\gamma\in(0,2]$, we have

$\|\theta\|_{L^{\infty}(Q_{1/2,\tau})}\lesssim_{\gamma}\left(1+\frac{1}{\tau-T}\right)^{\frac{n+2}{2\gamma}}\lVert\theta\rVert_{L^{\gamma}(Q_{1,T})}.$
(3.2)

###### Proof.

For smooth $\lambda:[T,0]\to(0,+\infty)$, we define $x=\lambda y$ and

$\widetilde{\theta}(y,t)=\theta(\lambda y,t).$ (3.3)

That is, $\widetilde{\theta}$ is obtained by dynamically rescaling $\theta$ in
space. The new PDE is

$\partial_{t}\widetilde{\theta}-\frac{1}{\lambda^{2}}\Delta_{y}\widetilde{\theta}+\frac{1}{\lambda}\widetilde{b}\cdot\nabla_{y}\widetilde{\theta}\leq
0$ (3.4)

where

$\widetilde{b}(y,t)=b(\lambda y,t)-\dot{\lambda}y.$ (3.5)

Choose $\lambda(T)=1$, $\dot{\lambda}=-2\|b(\cdot,t)\|_{L^{\infty}}$ when
$t\in[T,0]$. Clearly, $3/4\leq\lambda\leq 1$. Our picture is that
$\widetilde{\theta}$ dynamically ‘zooms in’ on $\theta$. In particular, using
(3.1) and (3.5),

$\widetilde{b}(\cdot,t)\cdot\frac{y}{|y|}\geq-\|b(\cdot,t)\|_{L^{\infty}}+2\|b(\cdot,t)\|_{L^{\infty}}|y|\geq
0\text{ when }y\in B_{1}\setminus B_{1/2},$ (3.6)

and

$\operatorname{div}\widetilde{b}=2n\|b(\cdot,t)\|_{L^{\infty}}\geq 0.$ (3.7)

We now demonstrate Cacciopoli’s inequality in the new variables. Let
$3/4\leq\varrho<R\leq 1$. Let $\varphi\in C^{\infty}_{0}(B_{R})$ be a radially
symmetric and decreasing function satisfying $0\leq\varphi\leq 1$ on
$\mathbb{R}^{n}$, $\varphi\equiv 1$ on $B_{\varrho}$, and
$|\nabla\varphi|\lesssim 1/(R-\varrho)$. Let $\eta\in
C^{\infty}_{0}(T,+\infty)$ satisfying $0\leq\eta\leq 1$ on $\mathbb{R}$,
$\eta\equiv 1$ on $(\tau,0)$, and $|d\eta/dt|\lesssim 1/(\tau-T)$. Let
$\Phi=\varphi^{2}\eta$. We integrate Eq. (3.4) against
$\widetilde{\theta}\Phi$ on $B_{R}\times(T,t)$ for $t\in(\tau,0)$. Then

$\displaystyle\frac{1}{2}\int|\widetilde{\theta}|^{2}(y,t)\Phi\,dy+\int\kern-3.39996pt\int\lambda^{-2}|\nabla\widetilde{\theta}|^{2}\Phi\,dy\,ds$
(3.8)
$\displaystyle\quad\leq\frac{1}{2}\int\kern-3.39996pt\int|\widetilde{\theta}|^{2}\partial_{t}\Phi\,dy\,ds-\int\kern-3.39996pt\int\lambda^{-2}\widetilde{\theta}\nabla\widetilde{\theta}\cdot\nabla\Phi\,dy\,ds$
$\displaystyle\quad\quad+\frac{1}{2}\int\kern-3.39996pt\int\lambda^{-1}\underbrace{\widetilde{b}\cdot\nabla\Phi}_{\leq
0\text{ by
}\eqref{eq:tildebpointsoutward}}|\widetilde{\theta}|^{2}\,dy\,ds+\frac{1}{2}\int\kern-3.39996pt\int\lambda^{-1}\operatorname{div}\widetilde{b}|\widetilde{\theta}|^{2}\Phi\,dy\,ds.$

While $\operatorname{div}\widetilde{b}$ has a disadvantageous sign, it simply
acts as a bounded potential. Simple manipulations give the Cacciopoli
inequality:

$\displaystyle\sup_{t\in(\tau,0)}\int_{B_{\varrho}}|\widetilde{\theta}|^{2}(y,t)\,dy+\int\kern-3.39996pt\int_{Q_{\varrho,\tau}}|\nabla\widetilde{\theta}|^{2}\varphi^{2}\,dy\,ds$
(3.9)
$\displaystyle\quad\lesssim\left(\frac{1}{\tau-T}+\frac{1}{(R-\varrho)^{2}}\right)\int\kern-3.39996pt\int_{Q_{R,T}\setminus
Q_{\varrho,\tau}}|\widetilde{\theta}|^{2}\,dy\,ds.$

The remainder of the proof proceeds as in Theorem 2.1 except in the $(y,t)$
variables. Namely, we have the interpolation inequality as in Corollary 2.8,
and $\widetilde{\theta}^{\beta}$ is a subsolution of (3.4) whenever $\beta\geq
1$. Therefore, we may perform Moser’s iteration verbatim. As in (2.42), the
$L^{2}$ norm on the RHS may be replaced by the $L^{\gamma}$ norm. Finally,
undoing the transformation yields the inequality (3.2) in the $(x,t)$
variables, since $(y,t)\in B_{R}\times\\{t\\}$ corresponds to $(x,t)\in
B_{\lambda(t)R}\times\\{t\\}$. ∎

The quantitative local boundedness property in Theorem 1.1 follows from
applying Proposition 3.1 and its rescalings on finitely many small time
intervals. In Remark 4.3, we justify that the constant depends on the
‘profile’ of $b$ and not just its norm.

##  4\. Counterexamples

###  4.1. Elliptic counterexamples

Let $n\geq 3$. Our counterexamples will be axisymmetric in ‘slab’ domains
$B_{R}\times(0,1)$, where $R>0$ is arbitrary and $B_{R}$ is a ball in
$\mathbb{R}^{n-1}$. We use the notation $x=(x^{\prime},z)$, where
$x^{\prime}\in\mathbb{R}^{n-1}$, $r=|x^{\prime}|$, and $z\in(0,1)$. Let

$\underaccent{\bar}{\theta}(x)=u(r)z,$ (4.1) $b(x)=V(r)e_{z}.$ (4.2)

Since $b$ is a shear flow in the $e_{z}$ direction, it is divergence free.
Then

$-\Delta\underaccent{\bar}{\theta}+b\cdot\nabla\underaccent{\bar}{\theta}=-z\Delta_{x^{\prime}}u+Vu.$
(4.3)

We will a construct a subsolution $\underaccent{\bar}{\theta}$ and
supersolution $\bar{\theta}$ using the _steady Schrödinger equation_

$-\Delta u+Vu=0$ (4.4)

in dimension $n-1$, where additionally $u\geq 0$ and $V\leq 0$. By the scale
invariance in $x^{\prime}$, it will suffice to construct a solution at a
single fixed length scale $R=R_{0}$. The way to proceed is well known. We
define

$u=\log\log\frac{1}{r},$ (4.5) $V=\frac{\Delta u}{u},$ (4.6)

for $r\leq R_{0}\ll 1$ so that $u$ is well defined. A simple calculation
verifies that $0\leq u\in H^{1}_{\rm loc}$, $\Delta u,V\leq 0$, and $\Delta
u,V\in L^{(n-1)/2}_{\rm loc}$.777Since Schrödinger solutions with critical
potentials $V$ belong to $L^{p}_{\rm loc}$ for all $p<\infty$ (see Han and Lin
[HL97], Theorem 4.4), it is natural to choose $u$ with a $\log$. The double
$\log$ ensures that $u$ has finite energy when $n=3$. Notice also that $\Delta
u=-(r\log r^{-1})^{-2}$ when $n=3$. Therefore, $Vu=\Delta u\in L^{1}_{\rm
loc}$, and the PDE (4.4) is satisfied in the sense of distributions. Using
(4.3), we verify that $\underaccent{\bar}{\theta}$ is a distributional
subsolution:

$-\Delta\underaccent{\bar}{\theta}+b\cdot\nabla\underaccent{\bar}{\theta}=-z\Delta_{x^{\prime}}u+Vu=(1-z)Vu\leq
0\text{ in }B_{R_{0}}\times(0,1),$ (4.7)

with equality at $\\{z=1\\}$. We also wish to control solutions from above.
Since $\Delta u\leq 0$, we define

$\bar{\theta}(x^{\prime},z)=u(r).$ (4.8)

Clearly, $\underaccent{\bar}{\theta}\leq\bar{\theta}$, and $\bar{\theta}$ is a
distributional supersolution:

$-\Delta\bar{\theta}+b\cdot\nabla\bar{\theta}=-\Delta_{x^{\prime}}u\geq
0\text{ in }B_{R_{0}}\times(0,1).$ (4.9)

We now construct smooth subsolutions and supersolutions approximating
$\underaccent{\bar}{\theta}$ and $\bar{\theta}$ according to the above
procedure. Let $\varphi$ be standard mollifier and

$\varphi_{\varepsilon}=\frac{1}{\varepsilon^{n-1}}\varphi\left(\frac{\cdot}{\varepsilon}\right).$
(4.10)

Define $u_{\varepsilon}=\varphi_{\varepsilon}\ast u$, $V_{\varepsilon}=\Delta
u_{\varepsilon}/u_{\varepsilon}$, $b_{\varepsilon}=V_{\varepsilon}(r)e_{z}$,
$\underaccent{\bar}{\theta}_{\varepsilon}=zu_{\varepsilon}(r)$, and
$\bar{\theta}_{\varepsilon}=u_{\varepsilon}(r)$. Then
$(\underaccent{\bar}{\theta}_{\varepsilon})$ and
$(\bar{\theta}_{\varepsilon})$ trap a family $(\theta_{\varepsilon})$ of
smooth solutions to the PDEs

$-\Delta\theta_{\varepsilon}+b_{\varepsilon}\cdot\nabla\theta_{\varepsilon}=0\text{
on }B_{R_{0}/2}\times(0,1)$ (4.11)

when $\varepsilon\in(0,R_{0}/2)$. Moreover, we have the desired estimates

$\sup_{\varepsilon\in(0,R_{0}/2)}\|\theta_{\varepsilon}\|_{L^{p}(B_{R_{0}/2}\times(0,1))}\leq\|\bar{\theta}\|_{L^{p}(B_{R_{0}}\times(0,1))}<+\infty,\quad
p\in[1,+\infty),$ (4.12)
$\sup_{\varepsilon\in(0,R_{0}/2)}\|V_{\varepsilon}\|_{L^{\frac{n-1}{2}}(B_{R_{0}/2})}\lesssim\sup_{\varepsilon}\|\Delta
u_{\varepsilon}\|_{L^{\frac{n-1}{2}}(B_{R_{0}/2})}\lesssim\|\Delta
u\|_{L^{\frac{n-1}{2}}(B_{R_{0}})}<+\infty,$ (4.13)

and the singularity as $\varepsilon\to 0^{+}$:

$\sup_{B_{\frac{R_{0}}{4}}\times(\frac{1}{4},\frac{3}{4})}\theta_{\varepsilon}\geq\sup_{B_{\frac{R_{0}}{4}}\times(\frac{1}{4},\frac{3}{4})}\underaccent{\bar}{\theta}_{\varepsilon}\to+\infty\text{
as }\varepsilon\to 0^{+}.$ (4.14)

######  Remark 4.1 (Line singularity).

The solutions constructed above are singular on the $z$-axis, as the maximum
principle demands.

######  Remark 4.2 (Time-dependent examples).

The above analysis of unbounded solutions for the steady Schrödinger equation
with critical potential is readily adapted to the parabolic PDE
$\partial_{t}u-\Delta u+Vu=f$ in
$B_{R_{0}}\times(-T_{0},0)\subset\mathbb{R}^{n+1}$, $n\geq 2$, (i) with
potential $V$ belonging to $L^{q}_{t}L^{p}_{x}$, $2/q+n/p=2$, $q>1$, and zero
force, or (ii) with force $f$ belonging to the same space and zero potential.
For example, one can define $u=\log\log(-t+r^{2})$, $V=-(\partial_{t}u-\Delta
u)/u$, and $f=0$. The case $q=1$ is an endpoint case in which solutions remain
bounded. These examples are presumably well known, although we do not know a
suitable reference.

###  4.2. Parabolic counterexamples

###### Proof of borderline cases: $L^{q}_{t}L^{p}_{x}$,
$\frac{2}{q}+\frac{n}{p}=2$, $q>1$.

1\. A heat subsolution. Let

$\Gamma(x,t)=(4\pi t)^{-n/2}e^{-\frac{|x|^{2}}{4t}}$ (4.15)

be the heat kernel. Let

$E(x,t)=(\Gamma-c_{n})_{+},$ (4.16)

where $c_{n}=(8\pi)^{-n/2}$. Then $E$ is globally Lipschitz away from $t=0$,
and $E(\cdot,t)$ is supported in the ball $B_{R(t)}$, where

$R(t)^{2}=2nt\log\frac{2}{t},\quad t<2,$ (4.17)

and $E$ vanishes in $t\geq 2$.

2\. A steady, compactly supported drift. There exists a divergence-free vector
field $U\in C^{\infty}_{0}(B_{4})$ satisfying

$U\equiv\vec{e}_{1}\text{ when }|x|\leq 2.$ (4.18)

Here is a construction: Let $\phi\in C^{\infty}_{0}(B_{4})$ be a radially
symmetric cut-off function such that $\phi\equiv 1$ on $B_{3}$. By applying
Bogovskii’s operator in the annulus $B_{4}\setminus\overline{B_{2}}$, see
[Gal11, Theorem III.3.3, p. 179], there exists $W\in
C^{\infty}_{0}(B_{4}\setminus\overline{B_{2}})$ solving

$\operatorname{div}W=-\operatorname{div}(\phi\vec{e}_{1})\in
C^{\infty}_{0}(B_{4}\setminus\overline{B_{2}}).$ (4.19)

Notably, the property of compact support is preserved. Finally, we define

$U=\phi\vec{e}_{1}+W.$ (4.20)

3\. Building blocks. Let $0\leq S\in C^{\infty}_{0}(0,1)$ and
$X\colon\mathbb{R}\to\mathbb{R}^{n}$ be the solution of the ODE

$\dot{X}(t)=S(t)\vec{e}_{1},\quad X(0)=-10n\vec{e}_{1}.$ (4.21)

Define

$b_{S}(x,t)=S(t)U\left(\frac{x-X(t)}{R(t)}\right),$ (4.22)

where $R(t)$ was defined in (4.17) above, and

$E_{S}(x,t)=E(x-X(t),t).$ (4.23)

Then $E_{S}$ is a subsolution:

$\left(\partial_{t}-\Delta+b_{S}\cdot\nabla\right)E_{S}\leq 0\quad\text{ on
}\mathbb{R}^{n}\times(0,1).$ (4.24)

If $[a,a^{\prime}]\subset(0,1)$, $S\in C^{\infty}_{0}(a,a^{\prime})$, and
$\int S\,dt\geq 20n$, then we have $E_{S}(\cdot,t)|_{B_{3}}\equiv 0$ when
$t\leq a$ or $t\geq a^{\prime}$. Additionally,
$E_{S}(\cdot,\widetilde{t})=E(\cdot,\widetilde{t})$ for some
$\widetilde{t}\in(a,a^{\prime})$.

We also consider the solution $\Phi_{S}$ to the PDE:

$\displaystyle\left(\partial_{t}-\Delta+b_{S}\cdot\nabla\right)\Phi_{S}$
$\displaystyle=0\quad\text{ on }\mathbb{R}^{n+1}\setminus\\{(X(0),0)\\}$
(4.25) $\displaystyle\Phi_{S}|_{t=0}$ $\displaystyle=\delta_{x=X(0)},$

which for short times $|t|\ll 1$ and negative times is equal to the heat
kernel $\Gamma(x-X(0),t)$. By the comparison principle,

$E_{S}\leq\Phi_{S}.$ (4.26)

We have the following measurements on the size of the drift:

$\|b_{S}\|_{L^{q}_{t}L^{p}_{x}(\mathbb{R}^{n+1})}^{q}\lesssim\int_{\mathbb{R}}S(t)^{q}R(t)^{\frac{nq}{p}}\,dt,$
(4.27)

where $1\leq p,q<+\infty$, and

$\|b_{S}\|_{L^{1}_{t}L^{\infty}_{x}(\mathbb{R}^{n+1})}=\|U\|_{L^{\infty}(\mathbb{R}^{n})}\int_{\mathbb{R}}S(t)\,dt.$
(4.28)

4\. Large displacement. For $S_{k}\in C^{\infty}_{0}(t_{k},t_{k}^{\prime})$,
$k\geq 1$, with $t_{k}=o_{k\to+\infty}(1)$ and
$[t_{k},t_{k}^{\prime}]\subset(0,1)$ disjoint-in-$k$, we consider the drifts
$b_{S_{k}}$. Let $M>0$. We claim that it is possible to choose $S_{k}$
satisfying

$\lVert b_{S_{k}}\rVert_{L^{1}_{t}L^{\infty}_{x}(\mathbb{R}^{n+1})}=M$ (4.29)

and

$\left\|\sum_{k\geq
1}b_{S_{k}}\right\|_{L^{q}_{t}L^{p}_{x}(\mathbb{R}^{n+1})}<+\infty$ (4.30)

for all $p,q\in[1,+\infty]$ satisfying $2/q+n/p=2$ and $q>1$. Indeed, consider

$\bar{S}(t)=\left(t\log t^{-1}\log\log t^{-1}\right)^{-1}$ (4.31)

when $t\leq c_{0}$ so that the above expression is well defined, and extended
smoothly on $[c_{0},1]$. We ask also that $t_{1}\leq c_{0}$. Since

$\int_{a}^{a^{\prime}}\bar{S}(t)dt=\log\log\log
t^{-1}\big{|}_{a}^{a^{\prime}}$ (4.32)

when $a^{\prime}\leq c_{0}$, we have

$\int_{t=0}^{1}\bar{S}(t)\,dt=+\infty,$ (4.33)

whereas

$\displaystyle\int_{t=0}^{1}\bar{S}(t)^{q}R(t)^{\frac{nq}{p}}\,dt$
$\displaystyle\leq O(1)+C_{n}\int_{t=0}^{c_{0}}(t\log
t^{-1})^{-q+\frac{nq}{2p}}(\log\log t^{-1})^{-q}\,dt$ (4.34)
$\displaystyle\leq O(1)+C_{n}\int_{t=0}^{c_{0}}(t\log t^{-1})^{-1}(\log\log
t^{-1})^{-q}\,dt<+\infty$

when $q\in(1,+\infty)$. The case $q=+\infty$ is similar. We choose
$S_{k}=\bar{S}(t)\varphi_{k}$ with suitable smooth cut-offs $\varphi_{k}$ to
complete the proof of the claim.

5\. Unbounded solution. We choose $M=20n$ and a suitable sequence of $S_{k}$
as above. We reorder the building blocks we defined above so that the $k$th
subsolution and $k$th drift are ‘activated’ on times
$(1-t_{k}^{\prime},1-t_{k})$. Define

$b_{k}(\cdot,t)=b_{S_{k}}(\cdot,t-(1-t_{k}^{\prime})+t_{k}),\quad
b=\sum_{k\geq 1}b_{k}$ (4.35)

and, for size parameters $A_{k}\geq 0$,

$E_{k}(\cdot,t)=E_{S_{k}}(\cdot,t-(1-t_{k}^{\prime})+t_{k})\mathbf{1}_{(1-t_{k}^{\prime},1-t_{k})},\quad
E=\sum_{k\geq 1}A_{k}E_{k}.$ (4.36)

Then $E$ is a subsolution of the PDE

$(\partial_{t}-\Delta+b\cdot\nabla)E\leq 0\quad\text{ on
}B_{3}\times(-\infty,1).$ (4.37)

We further define

$\Phi_{k}(\cdot,t)=\Phi_{S_{k}}(\cdot,t-(1-t_{k}^{\prime})+t_{k}),\quad\theta=\sum_{k\geq
1}A_{k}\Phi_{k},$ (4.38)

which is a solution of the PDE

$(\partial_{t}-\Delta+b\cdot\nabla)\theta=0\quad\text{ on
}B_{2}\times(-\infty,1).$ (4.39)

Since $E_{k}\leq\Phi_{k}$, we have that $E\leq\theta$ on
$B_{2}\times(-\infty,1)$. Additionally, we have

$\sup_{t}A_{k}\|E_{k}(\cdot,t)\|_{L^{\infty}(B_{1})}\geq
A_{k}\|E_{k}(\cdot,1-\widetilde{t}_{k})\|_{L^{\infty}(B_{1})}\gtrsim
A_{k}{t_{k}^{\prime}}^{-n/2},$ (4.40)

where $\widetilde{t}_{k}\in(t_{k},t_{k}^{\prime})$ satisfies
$X_{S_{k}}(\widetilde{t_{k}})=0$. Therefore, by the comparison principle and
(4.40), we have

$\limsup_{t\to
1_{-}}\|\theta(\cdot,t)\|_{L^{\infty}(B_{1})}\gtrsim\limsup_{k\to+\infty}A_{k}{t_{k}^{\prime}}^{-n/2}.$
(4.41)

To control the solution from above, we use

$\|\theta\|_{L^{\infty}_{t}L^{1}_{x}(\mathbb{R}^{n+1})}\leq\sum A_{k}.$ (4.42)

Therefore, it is possible to choose $A_{k}\to 0$ as $k\to+\infty$ while
keeping the $\limsup$ in (4.41) infinite. Hence, by ‘pruning’ the sequence of
$A_{k}$ (meaning we pass to a subsequence, without relabeling), we can always
ensure that $\|\theta\|_{L^{\infty}_{t}L^{1}_{x}(\mathbb{R}^{n+1})}<+\infty$.∎

######  Remark 4.3.

The sequence of solutions $\\{\theta_{k}\\}$ above demonstrates that the
constant in the quantitative local boundedness property in Theorem 1.1 for
drifts $b\in L^{1}_{t}L^{\infty}_{x}$ depends on the ‘profile’ of $b$ rather
than just its norm.

###### Proof of non-borderline cases: $L^{p}_{x}L^{q}_{t}$,
$\frac{3}{q}+\frac{n-1}{p}>2$, $p\leq q$.

This construction exploits rescaled copies of $E$ and is, in a certain sense,
self-similar.

_1\. Building blocks_. Let $(t_{k})\subset(0,1)$ be an increasing sequence,
with $t_{k}\to 1$ as $k\to+\infty$. Define $I_{k}=(t_{k},t_{k+1})$,
$R_{k}^{2}=|I_{k}|$.

Let $0\leq S\in C^{\infty}(0,1)$ satisfying $\int_{0}^{1}S(t)\,dt=M$ with
$M=20n$. Define $X_{k}\colon\mathbb{R}\to\mathbb{R}^{n}$ to be the solution of
the ODE

$\dot{X}_{k}(t-t_{k})=\underbrace{\frac{1}{|I_{k}|}S\left(\frac{t}{|I_{k}|}\right)}_{=:S_{k}(t-t_{k})}\vec{e}_{1},\quad
X_{k}(t_{k})=-10n\vec{e}_{1}.$ (4.43)

The ‘total speed’ has been normalized: $\int|\dot{X}_{k}|\,dt=\int S\,dt=M$.
Define also

$b_{k}(x,t)=S_{k}(t)U\left(\frac{x-X_{k}(t)}{R_{k}}\right)$ (4.44)

and

$E_{k}(x,t)=\frac{1}{R_{k}^{n}}E\left(\frac{x-X_{k}(t)}{R_{k}},\frac{t-t_{k}}{|I_{k}|}\right).$
(4.45)

Then $E_{k}$ is a subsolution:

$\left(\partial_{t}-\Delta+b_{k}\cdot\nabla\right)E_{k}\leq 0\text{ on
}\mathbb{R}^{n+1}\setminus\\{(X_{k}(t_{k}),t_{k})\\}$ (4.46)

and satisfies many of the same properties as $E_{k}$ in the previous
construction, among which is

$E_{k}(\cdot,\widetilde{t_{k}})=\frac{1}{R_{k}^{n}}E\left(\frac{\cdot}{R_{k}},\widetilde{t}\right)$
(4.47)

for some $\widetilde{t_{k}}\in(t_{k},t_{k+1})$ and $\widetilde{t}\in(0,1)$. We
define the solution $\theta_{k}$ to the PDE:

$\displaystyle\left(\partial_{t}-\Delta+b_{k}\cdot\nabla\right)\theta_{k}$
$\displaystyle=0\quad\text{ on
}\mathbb{R}^{n+1}\setminus\\{(X_{k}(t_{k}),t_{k})\\}$ (4.48)
$\displaystyle\theta_{k}|_{t=t_{k}}$ $\displaystyle=\delta_{x=X_{k}(t_{k})},$

which for short times $|t-t_{k}|\ll_{k}1$ and times $t<t_{k}$ is equal to the
heat kernel $\Gamma(x-X_{k}(t_{k}),t)$. The comparison principle implies

$E_{k}\leq\theta_{k}.$ (4.49)

2\. Estimating the drift. We now estimate the size of $b_{k}$. To begin, we
estimate the $L^{q}_{t}L^{p}_{x}$ norms, $\frac{2}{q}+\frac{n}{p}>2$. Using
the scalings from (4.44), we have

$\max|b_{k}|\leq\|U\|_{L^{\infty}}\|S\|_{L^{\infty}}|I_{k}|^{-1}$ (4.50)

and

$\|b_{k}\|_{L^{q}_{t}L^{p}_{x}(\mathbb{R}^{n+1})}\lesssim\|U\|_{L^{\infty}}\|S\|_{L^{\infty}}|I_{k}|^{\frac{1}{q}-1}R_{k}^{\frac{n}{p}}\lesssim
R_{k}^{\varepsilon(p,q)}=o_{k\to+\infty}(1),$ (4.51)

since $|I_{k}|=R_{k}^{2}$. Next, we estimate the
$L^{p}_{x^{\prime}}L^{\infty}_{x_{n}}L^{\widetilde{q}}_{t}$ norm, where
$\frac{2}{\widetilde{q}}+\frac{n-1}{p}>2$. We are most interested when
$\widetilde{q}=+\infty$ and $p=\frac{n-1}{2}-$, but it is not more effort to
estimate this. Importantly, we have

${\rm supp}\;b_{k}\subset
B^{\mathbb{R}^{n-1}}_{C_{n}R_{k}}\times(-C_{n},C_{n})\times I_{k}.$ (4.52)

Using this and (4.50), we have

$\|b_{k}\|_{L^{p}_{x^{\prime}}L^{\infty}_{x_{n}}L^{\widetilde{q}}_{t}(\mathbb{R}^{n+1})}\lesssim\|U\|_{L^{\infty}}\|S\|_{L^{\infty}}R_{k}^{\frac{n-1}{p}}|I_{k}|^{\frac{1}{\widetilde{q}}-1}\lesssim
R_{k}^{\varepsilon(p,\widetilde{q})}=o_{k\to+\infty}(1).$ (4.53)

Interpolating between (4.51) and (4.53) with
$(p,\widetilde{q})=(\frac{n-1}{2}-,\infty)$, we thus obtain

$\|b_{k}\|_{L^{p}_{x}L^{q}_{t}(\mathbb{R}^{n+1})}=o_{k\to+\infty}(1)$ (4.54)

when $\frac{3}{q}+\frac{n-1}{p}>2$ and $p\leq q$. After ‘pruning’ the sequence
in $k$ (meaning we pass to a subsequence, without relabeling), we have

$\|b\|_{L^{q}_{t}L^{p}_{x}(\mathbb{R}^{n+1})}\leq\sum\|b_{k}\|_{L^{q}_{t}L^{p}_{x}(\mathbb{R}^{n+1})}<+\infty,\quad\frac{2}{q}+\frac{n}{p}>2$
(4.55)

and

$\|b\|_{L^{p}_{x}L^{q}_{t}(\mathbb{R}^{n+1})}\leq\sum\|b_{k}\|_{L^{p}_{x}L^{q}_{t}(\mathbb{R}^{n+1})}<+\infty,\quad\frac{3}{q}+\frac{n-1}{p}>2,\quad
p\leq q.$ (4.56)

_3\. Concluding_. The remainder of the proof proceeds as before, with the
notable difference that we do not need to reorder the blocks in time. To
summarize, we have

$\sup_{t}A_{k}\|E_{k}(\cdot,t)\|_{L^{\infty}(B_{1})}\geq
A_{k}\|E_{k}(\cdot,\widetilde{t_{k}})\|_{L^{\infty}(B_{1})}\gtrsim
A_{k}{R_{k}}^{-n},$ (4.57)

where $\widetilde{t}_{k}\in(t_{k},t_{k+1})$ satisfies
$X_{k}(\widetilde{t_{k}})=0$, and hence,

$\limsup_{t\to
1_{-}}\|\theta(\cdot,t)\|_{L^{\infty}(B_{1})}\gtrsim\limsup_{k\to+\infty}A_{k}{R_{k}}^{-n}$
(4.58)

To control the solution from above, we again use (4.42) and choose $A_{k}\to
0$ as $k\to+\infty$ while maintaining that the RHS (4.58) is infinite. By
again ‘pruning’ the sequence in $k$, we have
$\|\theta\|_{L^{\infty}_{t}L^{1}_{x}(\mathbb{R}^{n+1})}<+\infty$. This
completes the proof. ∎

######  Remark 4.4 (An open question).

As mentioned in the introduction, we do not construct counterexamples in the
endpoint cases $L^{p}_{x}L^{q}_{t}$, $\frac{3}{q}+\frac{n-1}{p}=2$, except
when $p=q=\frac{n+2}{2}$ or $(p,q)=(\frac{n-1}{2},+\infty)$ (steady example
constructed above). This seems to suggest, perhaps, that local boundedness
should also fail on the line between these two points, but that the
counterexamples may be more subtle. It would be interesting to construct these
examples. Since each ‘block’ above is uniformly bounded in the desired spaces,
we can say that, if local boundedness were to hold there, it must be depend on
the ‘profile’ of $b$ and not just its norm, as in Remark 4.3.

##  5\. Upper bounds on fundamental solutions

For the Gaussian-like upper bounds on fundamental solutions, we consider a
divergence-free vector field $b\in
C^{\infty}_{0}(\mathbb{R}^{n}\times[0,+\infty))$ and $\alpha_{0}\geq 0$,
$C_{0},M_{0}>0$ satisfying the following two properties:

I. _Local boundedness_ : For each $x_{0}\in\mathbb{R}^{n}$, $R>0$,
$t_{0}\in\mathbb{R}_{+}$, parabolic cylinder
$Q_{R}(x_{0},t_{0})=B_{R}(x_{0})\times(t_{0}-R^{2},t_{0})\subset\mathbb{R}^{n}\times\mathbb{R}_{+}$,
and Lipschitz solution $u\in W^{1,\infty}(Q_{R}(x_{0},t_{0}))$, we have

$\|u\|_{L^{\infty}(Q_{R/2}(x_{0},t_{0}))}\leq
C_{0}\left(\frac{M_{0}}{R^{\alpha_{0}+2}}+\frac{1}{R^{2}}\right)^{\frac{n+2}{4}}\|u\|_{L^{2}(Q_{R}(x_{0},t_{0}))}$
(5.1)

II. _Global, revised form boundedness condition_ , which we denote
($\widetilde{\text{FBC}}$):

> For each $x_{0}\in\mathbb{R}^{n}$, $R>0$, open interval
> $I\subset\mathbb{R}_{+}$, and Lipschitz function $u\in
> W^{1,\infty}(B_{R}\times I)$, there exists a measurable set
> $A=A(x_{0},R,b)\subset(R/2,R)$ with $|A|\geq R/4$ and satisfying
>
> $\displaystyle\frac{1}{|A|}\int\kern-3.39996pt\int_{B_{A}\times
> I}\frac{|u|^{2}}{2}(b\cdot
> n)\,dx\,dt\leq\frac{M_{0}}{\varepsilon^{\alpha_{0}+1}R^{\alpha_{0}+2}}\int\kern-3.39996pt\int_{B_{R}\times
> I}|u|^{2}\,dx\,dt$ ($\widetilde{\text{FBC}}$)
> $\displaystyle\quad+\varepsilon\left(\frac{1}{R^{2}}\int\kern-3.39996pt\int_{B_{R}\times
> I}|u|^{2}\,dx\,dt+\int\kern-3.39996pt\int_{B_{R}\times I}|\nabla
> u|^{2}\,dx\,dt+\sup_{t\in I}\int_{B_{R}}|u(x,t)|^{2}\,dx\right),$
>
> where $n$ is the unit outer normal direction and $B_{A}=\cup_{r\in
> A}\partial B_{r}$.

Under the above conditions, we have

######  Theorem 5.1 (Upper bounds on fundamental solutions).

If $b\in C^{\infty}_{0}(\mathbb{R}^{n}\times[0,+\infty))$ is a divergence-free
vector field satisfying I and II, then the fundamental solution
$\Gamma=\Gamma(x,t;y,s)$ to the parabolic operator
$L=\partial_{t}-\Delta+b\cdot\nabla$ satisfies the following upper bounds:

$\displaystyle\Gamma(x,t;y,s)\leq
C(t-s)\left(\frac{1}{t-s}+\frac{M_{0}}{(t-s)^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}$
(5.2)
$\displaystyle\quad\times\left[\exp\left(-\frac{c|x-y|^{2}}{t-s}\right)+\exp\left(-\frac{c|x-y|^{1+\frac{1}{1+\alpha_{0}}}}{((t-s)M_{0})^{\frac{1}{1+\alpha_{0}}}}\right)\right]$
(5.3)

for all $x,y\in\mathbb{R}^{n}$ and $0\leq s<t<+\infty$, where $C,c>0$ may
depend on $C_{0}$, $\alpha_{0}$, and $n$.

######  Remark 5.2 (Comments on I and II).

Notice that ($\widetilde{\text{FBC}}$) in II bounds the _outflux_ of scalar
through the domain. The _influx_ , which was bounded explicitly in (FBC), is
handled implicitly by I.

The parameter $\alpha_{0}$ in ($\widetilde{\text{FBC}}$) does not track the
same information as the parameter $\alpha$ in (FBC). The constant $M_{0}$ in
($\widetilde{\text{FBC}}$) has dimensions $L^{\alpha_{0}}$, whereas the
constant $M$ in (FBC) is dimensionless. The upper bound in (5.2) is
dimensionally correct.

Upon optimizing in $\varepsilon$, one discovers that
($\widetilde{\text{FBC}}$) is equivalent to an interpolation inequality, see
the derivation in Remark 2.6. The above form is reflective of how the estimate
is utilized below.

It would be possible to incorporate a parameter $\delta_{0}\in(0,1/2)$ such
that $|A|\geq\delta_{0}R$, as in (FBC).

We verify that the above conditions are satisfied in the setting of the main
theorems:

######  Lemma 5.3 (Verifying $\widetilde{\text{FBC}}$).

Let $p,q\in[1,+\infty]$ and $b\in
C^{\infty}_{0}(\mathbb{R}^{n}\times[0,+\infty))$ be a divergence-free vector
field.

1\. If

$b\in L^{q}_{t}L^{p}_{x}(\mathbb{R}^{n}\times\mathbb{R}_{+}),\quad
1\leq\zeta_{0}:=\frac{2}{q}+\frac{n}{p}<2,$ (5.4)

then $b$ satisfies I and II with
$\alpha_{0}=\frac{\zeta_{0}-1}{\theta_{2}}(=\frac{1}{\theta_{2}}-2)$,
$M_{0}=C\|b\|_{L^{q}_{t}L^{p}_{x}(\mathbb{R}^{n}\times\mathbb{R}_{+})}^{\frac{1}{\theta_{2}}}+\frac{1}{4}$,
$C_{0}=C$, and $\theta_{2}=1-\frac{\zeta_{0}}{2}$.

2\. If

$b\in L^{p}_{x}L^{q}_{t}(\mathbb{R}^{n}\times\mathbb{R}_{+}),\quad p\leq
q,\quad 1\leq\zeta_{0}:=\frac{3}{q}+\frac{n-1}{p}<2,$ (5.5)

then $b$ satisfies I and II with
$\alpha_{0}=\frac{\zeta_{0}-1}{\theta_{2}}(=\frac{1}{\theta_{2}}-2)$,
$M_{0}=C\|b\|_{L^{p}_{x}L^{q}_{t}(\mathbb{R}^{n}\times\mathbb{R}_{+})}^{\frac{1}{\theta_{2}}}+\frac{1}{4}$,
$C_{0}=C$, and $\theta_{2}=1-\frac{\zeta_{0}}{2}$.

###### Proof.

I follows from Proposition 2.3, Theorem 2.1, and the translation and scaling
symmetries. II follows from rescaling $R_{0}$ in Remark 2.6. ∎

###### Proof of Theorem 5.1.

By a translation of the coordinates, we need only consider
$\Gamma(x_{0},t_{0};0,0)$ for $t_{0}>0$ and $x_{0}\in\mathbb{R}^{n}$. We adapt
a method due to E. B. Davies [Dav87]. Let $\psi=\psi(|x|)=\psi(r)$ be a
bounded radial Lipschitz function to be specified such that $\psi=0$ when
$r<|x_{0}|/2$ and $\psi$ is a constant when $r\geq|x_{0}|$. For now, we record
the property $|\nabla\psi|\leq\gamma$, where $\gamma>0$ will be specified
later.

_1\. Weighted energy estimates_. Let $f\in C^{\infty}_{0}(\mathbb{R}^{n})$ and
$u$ be the solution to the equation $Lu=0$ in
$\mathbb{R}^{n}\times\mathbb{R}^{+}$ with initial condition
$u(0,\cdot)=e^{-\psi}f$. For $t\geq 0$, denote

$J(t)=\frac{1}{2}\int_{\mathbb{R}^{n}}e^{2\psi(y)}u^{2}(y,t)\,dy.$ (5.6)

Then, by integration by parts, for $t>0$, we have

$\displaystyle\dot{J}(t)$
$\displaystyle=\int_{\mathbb{R}^{n}}e^{2\psi}u\partial_{t}u\,dy=\int_{\mathbb{R}^{n}}e^{2\psi}u(\Delta
u-b\cdot\nabla u)\,dy$ $\displaystyle=-\int_{\mathbb{R}^{n}}e^{2\psi}|\nabla
u|^{2}\,dy-2\int_{\mathbb{R}^{n}}e^{2\psi}u\nabla
u\cdot\nabla\psi\,dy-\int_{\mathbb{R}^{n}}e^{2\psi}b\cdot\nabla\frac{u^{2}}{2}\,dy$
$\displaystyle\leq-\frac{1}{2}\int_{\mathbb{R}^{n}}e^{2\psi}|\nabla
u|^{2}\,dy+C\gamma^{2}J(t)+\int_{\mathbb{R}^{n}}e^{2\psi}b\cdot\nabla\psi
u^{2}\,dy,$

where we applied Young’s inequality. Hence,

$\displaystyle
J(t)+\frac{1}{2}\int_{0}^{t}\int_{\mathbb{R}^{n}}e^{2\psi}|\nabla
u|^{2}\,dy\,ds$ (5.7) $\displaystyle\quad\leq
J(0)+C\gamma^{2}\int_{0}^{t}J(s)\,ds+\int_{0}^{t}\int_{\mathbb{R}^{n}}e^{2\psi}b\cdot\nabla\psi
u^{2}\,dy\,ds.$

Now for each $t>0$, we choose $\psi^{\prime}(r)=\gamma\mathbf{1}_{A}$, where
$A=A(x_{0},|x_{0}|,b)$ is the set of ‘good slices’ from
($\widetilde{\text{FBC}}$) in II. By ($\widetilde{\text{FBC}}$) applied to
$e^{\psi}u$, we have

$\displaystyle
J(t)+\frac{1}{2}\int_{0}^{t}\int_{\mathbb{R}^{n}}e^{2\psi}|\nabla
u|^{2}\,dy\,ds$ $\displaystyle\leq
J(0)+C\gamma^{2}\int_{0}^{t}J(s)\,ds+\gamma\int_{0}^{t}\int_{B_{A}}e^{2\psi}b\cdot\frac{y}{|y|}u^{2}\,dy\,ds$
$\displaystyle\leq
J(0)+C\gamma^{2}\int_{0}^{t}J(s)\,ds+\gamma\varepsilon|A|\left(\sup_{s\in[0,t]}J(s)+\int_{0}^{t}\int_{B_{|x_{0}|}}e^{2\psi}|\nabla
u|^{2}\,dy\,ds\right)$
$\displaystyle\quad+\gamma\varepsilon|A|\int_{0}^{t}\int_{B_{|x_{0}|}}e^{2\psi}|\nabla\psi|^{2}u^{2}\,dy\,ds+\left(\frac{\gamma
M_{0}|A|}{\varepsilon^{\alpha_{0}+1}(|x_{0}|/2)^{\alpha_{0}+2}}+\frac{\gamma\varepsilon|A|}{|x_{0}|^{2}}\right)\int_{0}^{t}\int_{B_{|x_{0}|}}e^{2\psi}u^{2}\,dy\,ds.$

Recall that $|A|\leq|x_{0}|/2$. We set $\varepsilon=(100|x_{0}|\gamma)^{-1}$
and take the supremum in $t$ in order to absorb the 3rd term on the right-hand
side, which has coefficient $\gamma\varepsilon|A|$:

$\displaystyle\frac{1}{2}\sup_{s\in[0,t]}J(s)+\frac{1}{4}\int_{0}^{t}\int_{\mathbb{R}^{n}}e^{2\psi}|\nabla
u|^{2}\,dy\,ds$ $\displaystyle\quad\leq
J(0)+\left(C\gamma^{2}+CM_{0}\gamma^{2+\alpha_{0}}+|x_{0}|^{-2}\right)\int_{0}^{t}J(s)\,ds.$

Now, by the Gronwall inequality,

$J(t)\leq CJ(0)e^{(C\gamma^{2}+CM_{0}\gamma^{2+\alpha_{0}}+|x_{0}|^{-2})t}.$
(5.8)

_2\. Duality argument_. For $s<t$, we define an operator

$P^{\psi}_{s\to
t}f(x)=e^{\psi(x)}\int_{\mathbb{R}^{n}}\Gamma(x,t;y,s)e^{-\psi(y)}f(y)\,dy.$
(5.9)

By the local boundedness estimate (5.1) in I, for any $x\in\mathbb{R}^{n}$ and
$R=\sqrt{t}$,

$\displaystyle e^{-2\psi(x)}|P_{0\to t}^{\psi}f(x)|^{2}=|u(x,t))|^{2}\leq
C_{0}\left(\frac{1}{t}+\frac{M_{0}}{t^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}\int_{0}^{t}\int_{B_{\sqrt{t}}(x)}u^{2}(y,s)\,dy\,ds.$

Thus, allowing $C$ to depend on $C_{0}$, from (5.8) we have

$\displaystyle|P_{0\to t}^{\psi}f(x)|^{2}\leq
C\left(\frac{1}{t}+\frac{M_{0}}{t^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}\int_{0}^{t}\int_{B_{\sqrt{t}}(x)}e^{2\psi(x)}u^{2}(y,s)\,dy\,ds$
$\displaystyle\leq
C\left(\frac{1}{t}+\frac{M_{0}}{t^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}\sup_{y\in
B_{\sqrt{t}}(x)}e^{2(\psi(x)-\psi(y))}\int_{0}^{t}J(s)\,ds$
$\displaystyle\overset{\eqref{eq1.23}}{\leq}Ct\left(\frac{1}{t}+\frac{M_{0}}{t^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}e^{2\gamma\min\\{\sqrt{t},|x_{0}|/2\\}}\times
e^{(C\gamma^{2}+CM_{0}\gamma^{2+\alpha_{0}}+|x_{0}|^{-2})t}J(0).$

Since $J(0)=\|f\|_{L^{2}}^{2}/2$, the above inequality together with a
translation in time implies that

$\displaystyle\|P_{s\to t}^{\psi}\|_{L^{2}\to L^{\infty}}^{2}\leq
C(t-s)\left(\frac{1}{t-s}+\frac{M_{0}}{(t-s)^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}$
$\displaystyle\quad\cdot
e^{2\gamma\min\\{\sqrt{t-s},|x_{0}|/2\\}}e^{(C\gamma^{2}+CM_{0}\gamma^{2+\alpha_{0}}+|x_{0}|^{-2})(t-s)}J(0).$

By duality, we also have

$\displaystyle\|P_{s\to t}^{\psi}\|_{L^{1}\to L^{2}}^{2}\leq
C(t-s)\left(\frac{1}{t-s}+\frac{M_{0}}{(t-s)^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}$
$\displaystyle\quad\times
e^{2\gamma\min\\{\sqrt{t-s},|x_{0}|/2\\}}e^{(C\gamma^{2}+CM_{0}\gamma^{2+\alpha_{0}}+|x_{0}|^{-2})(t-s)}J(0).$

Therefore,

$\displaystyle\|P_{0\to t}^{\psi}\|_{L^{1}\to L^{\infty}}\leq\|P_{0\to
t/2}^{\psi}\|_{L^{1}\to L^{2}}\|P_{t/2\to t}^{\psi}\|_{L^{2}\to L^{\infty}}$
$\displaystyle\leq
Ct\left(\frac{1}{t}+\frac{M_{0}}{t^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}e^{2\gamma\min\\{\sqrt{t},|x_{0}|/2\\}}e^{(C\gamma^{2}+CM_{0}\gamma^{2+\alpha_{0}}+|x_{0}|^{-2})t}.$

From (5.9), the above inequality implies that for any $x,y\in\mathbb{R}^{n}$,

$\displaystyle e^{\psi(x)-\psi(y)}\Gamma(x,t;y,0)$ $\displaystyle\leq
Ct\left(\frac{1}{t}+\frac{M_{0}}{t^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}e^{2\gamma\min\\{\sqrt{t},|x_{0}|/2\\}}e^{(C\gamma^{2}+CM_{0}\gamma^{2+\alpha_{0}}+|x_{0}|^{-2})t}.$

In particular, we have

$\displaystyle\Gamma(x_{0},t;0,0)$ $\displaystyle\leq
Ce^{-\psi(x_{0})}t\left(\frac{1}{t}+\frac{M_{0}}{t^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}e^{2\gamma\min\\{\sqrt{t},|x_{0}|/2\\}}e^{(C\gamma^{2}+CM_{0}\gamma^{2+\alpha_{0}}+|x_{0}|^{-2})t}$
$\displaystyle\leq
Ct\left(\frac{1}{t}+\frac{M_{0}}{t^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}\exp\left\\{2\gamma\min\Big{(}\sqrt{t},\frac{|x_{0}|}{2}\Big{)}-\frac{\gamma|x_{0}|}{4}+C\gamma^{2}t+CM_{0}\gamma^{2+\alpha_{0}}t+|x_{0}|^{-2}t\right\\}$

since $\psi(0)=0$ and $\psi(x_{0})\geq\gamma|x_{0}|/4$.

_3\. Optimizing $\gamma$_. In order to have a negative term in the
exponential, we first consider

$|x_{0}|\geq 16\sqrt{t}.$ (5.10)

Then

$2\min\Big{(}\sqrt{t},\frac{|x_{0}|}{2}\Big{)}-\frac{|x_{0}|}{4}\leq-\frac{|x_{0}|}{8}$
(5.11)

Notice that final term in the exponential is already controlled:

$t|x_{0}|^{-2}\lesssim 1.$ (5.12)

Hence,

$\Gamma(x_{0},t;0,0)\lesssim
t\left(\frac{1}{t}+\frac{M_{0}}{t^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}\exp\left(-\frac{1}{8}\gamma|x_{0}|+C\gamma^{2}t+CM_{0}\gamma^{2+\alpha_{0}}t\right).$
(5.13)

The new expression inside the exponential is

$-\frac{1}{8}\gamma|x_{0}|+\underbrace{C\gamma^{2}t}_{A}+\underbrace{CM_{0}\gamma^{2+\alpha_{0}}t}_{B}.$
(5.14)

We divide space into an ‘inner region’ and ‘outer region’, and we anticipate
that the rate of decay may be modified in the outer region, where $B$
dominates.

_3a. Outer region_. We consider scalings of $\gamma$ in which
$-\frac{1}{8}\gamma|x_{0}|$ overtakes $B$. Consider

$M_{0}\gamma^{2+\alpha_{0}}t=\varepsilon\gamma|x_{0}|.$ (5.15)

Then $\frac{1}{8}\gamma|x_{0}|\gg B$ when

$\gamma^{1+\alpha_{0}}=\varepsilon M_{0}^{-1}\frac{|x_{0}|}{t}$ (5.16)

and $\varepsilon\ll 1$. In this scaling, we have that $B\geq A$
($CM_{0}\gamma^{2+\alpha_{0}}\geq C\gamma^{2}$) when

$M_{0}\gamma^{\alpha_{0}}\geq 1,$ (5.17)

or

$\frac{M_{0}^{\frac{1}{\alpha_{0}}}|x_{0}|}{t}\geq\varepsilon^{-1}.$ (5.18)

In this region, under the additional assumption (5.10), we have the
exponential bound

$\Gamma(x_{0},t;0,0)\lesssim
t\left(\frac{1}{t}+\frac{M_{0}}{t^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}\exp\left(-\frac{c|x_{0}|^{1+\frac{1}{1+\alpha_{0}}}}{(tM_{0})^{\frac{1}{1+\alpha_{0}}}}\right)$
(5.19)

_3b. Inner region_. We now consider scalings of $\gamma$ in which
$-\frac{1}{8}\gamma|x_{0}|$ overtakes the term $A$. Consider

$\gamma^{2}t=\varepsilon\gamma|x_{0}|.$ (5.20)

Then $\frac{1}{8}\gamma|x_{0}|\gg A$ when

$\gamma=\varepsilon|x_{0}|/t$ (5.21)

and $\varepsilon\ll 1$. In this scaling, we have that $A\geq B$
($C\gamma^{2}\geq CM_{0}\gamma^{2+\alpha_{0}}$) when

$M_{0}\gamma^{\alpha_{0}}\leq 1,$ (5.22)

or

$M_{0}^{\frac{1}{\alpha_{0}}}\frac{|x_{0}|}{t}\leq\varepsilon^{-1}.$ (5.23)

In this region, under the additional assumption (5.10), we have the
exponential bound

$\Gamma(x_{0},t;0,0)\lesssim
t\left(\frac{1}{t}+\frac{M_{0}}{t^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}\exp\left(-\frac{c|x_{0}|^{2}}{t}\right).$
(5.24)

_3c. Patching_. Combining (5.19) and (5.24), we have

$\Gamma(x_{0},t;0,0)\lesssim
t\left(\frac{1}{t}+\frac{M_{0}}{t^{1+\frac{\alpha_{0}}{2}}}\right)^{\frac{n+2}{2}}\left[\exp\left(-\frac{c|x_{0}|^{2}}{t}\right)+\exp\left(-\frac{c|x_{0}|^{1+\frac{1}{1+\alpha_{0}}}}{(tM_{0})^{\frac{1}{1+\alpha_{0}}}}\right)\right].$
(5.25)

under the assumption (5.10). On the other hand, when $|x_{0}|\leq 16\sqrt{t}$,
the fundamental solution is controlled by the Nash estimate [Nas58]:

$\|\Gamma(\cdot,t;0,0)\|_{L^{\infty}(\mathbb{R}^{n})}\lesssim
t^{-\frac{n}{2}},$ (5.26)

which is independent of the divergence-free drift. This contributes to the
prefactor in (5.25). The proof is complete. ∎

#### Acknowledgments

DA thanks Vladimír Šverák for encouraging him to answer this question and
helpful discussions. DA also thanks Tobias Barker and Simon Bortz for helpful
discussions, especially concerning Remark 4.2, and their patience. DA was
supported by the NDSEG Fellowship and NSF Postdoctoral Fellowship Grant No.
2002023. HD was partially supported by the Simons Foundation, Grant No.
709545, a Simons Fellowship, and the NSF under agreement DMS-2055244.

## References

  * [ABES19] Pascal Auscher, Simon Bortz, Moritz Egert, and Olli Saari.  On regularity of weak solutions to linear parabolic systems with measurable coefficients.  J. Math. Pures Appl. (9), 121:216–243, 2019.
  * [Aro67] D. G. Aronson.  Bounds for the fundamental solution of a parabolic equation.  Bull. Amer. Math. Soc., 73:890–896, 1967.
  * [BS20a] Peter Bella and Mathias Schäffner.  Non-uniformly parabolic equations and applications to the random conductance model.  arXiv preprint arXiv:2009.11535, 2020.
  * [BS20b] Peter Bella and Mathias Schäffner.  On the regularity of minimizers for scalar integral functionals with $(p,q)$-growth.  Anal. PDE, 13(7):2241–2257, 2020.
  * [BS21] Peter Bella and Mathias Schäffner.  Local boundedness and Harnack inequality for solutions of linear nonuniformly elliptic equations.  Comm. Pure Appl. Math., 74(3):453–477, 2021.
  * [Dav87] E. B. Davies.  Explicit constants for Gaussian upper bounds on heat kernels.  Amer. J. Math., 109(2):319–333, 1987.
  * [DG57] Ennio De Giorgi.  Sulla differenziabilità e l’analiticità delle estremali degli integrali multipli regolari.  Mem. Accad. Sci. Torino. Cl. Sci. Fis. Mat. Nat. (3), 3:25–43, 1957\.
  * [DK18] Hongjie Dong and Seick Kim.  Fundamental solutions for second-order parabolic systems with drift terms.  Proc. Amer. Math. Soc., 146(7):3019–3029, 2018.
  * [FG85] Eugene B. Fabes and Nicola Garofalo.  Parabolic B.M.O. and Harnack’s inequality.  Proc. Amer. Math. Soc., 95(1):63–69, 1985.
  * [Fil13] N. Filonov.  On the regularity of solutions to the equation $-\Delta u+b\cdot\nabla u=0$.  Zap. Nauchn. Sem. S.-Peterburg. Otdel. Mat. Inst. Steklov. (POMI), 410(Kraevye Zadachi Matematicheskoĭ Fiziki i Smezhnye Voprosy Teorii Funktsiĭ. 43):168–186, 189, 2013.
  * [FR96] Jens Frehse and Michael Ružička.  Existence of regular solutions to the steady Navier-Stokes equations in bounded six-dimensional domains.  Annali della Scuola Normale Superiore di Pisa-Classe di Scienze, 23(4):701–719, 1996.
  * [FS18] Nikolay Filonov and Timofey Shilkin.  On some properties of weak solutions to elliptic equations with divergence-free drifts.  In Mathematical analysis in fluid mechanics—selected recent results, volume 710 of Contemp. Math., pages 105–120. Amer. Math. Soc., [Providence], RI, [2018] ©2018.
  * [FSSC98] Bruno Franchi, Raul Serapioni, and Francesco Serra Cassano.  Irregular solutions of linear degenerate elliptic equations.  Potential Anal., 9(3):201–216, 1998.
  * [Gal11] G. P. Galdi.  An introduction to the mathematical theory of the Navier-Stokes equations.  Springer Monographs in Mathematics. Springer, New York, second edition, 2011.  Steady-state problems.
  * [Giu03] Enrico Giusti.  Direct methods in the calculus of variations.  World Scientific Publishing Co., Inc., River Edge, NJ, 2003.
  * [HL97] Qing Han and Fanghua Lin.  Elliptic partial differential equations, volume 1 of Courant Lecture Notes in Mathematics.  New York University, Courant Institute of Mathematical Sciences, New York; American Mathematical Society, Providence, RI, 1997.
  * [HLMP21] Steve Hofmann, Linhan Li, Svitlana Mayboroda, and Jill Pipher.  The dirichlet problem for elliptic operators having a BMO anti-symmetric part.  Mathematische Annalen, June 2021.
  * [Ign14] Mihaela Ignatova.  On the continuity of solutions to advection-diffusion equations with slightly super-critical divergence-free drifts.  Adv. Nonlinear Anal., 3(2):81–86, 2014.
  * [IKR14] Mihaela Ignatova, Igor Kukavica, and Lenya Ryzhik.  The Harnack inequality for second-order elliptic equations with divergence-free drifts.  Commun. Math. Sci., 12(4):681–694, 2014.
  * [IKR16] Mihaela Ignatova, Igor Kukavica, and Lenya Ryzhik.  The Harnack inequality for second-order parabolic equations with divergence-free drifts of low regularity.  Comm. Partial Differential Equations, 41(2):208–226, 2016.
  * [Kon07] Michalis Kontovourkis.  On elliptic equations with low-regularity divergence-free drift terms and the steady-state Navier-Stokes equations in higher dimensions.  University of Minnesota, 2007.
  * [KS19] Seick Kim and Georgios Sakellaris.  Green’s function for second order elliptic equations with singular lower order coefficients.  Comm. Partial Differential Equations, 44(3):228–270, 2019.
  * [Lie96] Gary M. Lieberman.  Second order parabolic differential equations.  World Scientific Publishing Co., Inc., River Edge, NJ, 1996.
  * [LP19] Linhan Li and Jill Pipher.  Boundary behavior of solutions of elliptic operators in divergence form with a BMO anti-symmetric part.  Comm. Partial Differential Equations, 44(2):156–204, 2019.
  * [LSUt68] O. A. Ladyženskaja, V. A. Solonnikov, and N. N. Ural’ tseva.  Linear and quasilinear equations of parabolic type.  Translations of Mathematical Monographs, Vol. 23. American Mathematical Society, Providence, R.I., 1968.  Translated from the Russian by S. Smith.
  * [LZ04] Vitali Liskevich and Qi S. Zhang.  Extra regularity for parabolic equations with drift terms.  Manuscripta Math., 113(2):191–209, 2004.
  * [Moo17] Connor Mooney.  Finite time blowup for parabolic systems in two dimensions.  Arch. Ration. Mech. Anal., 223(3):1039–1055, 2017.
  * [Mos61] Jürgen Moser.  On Harnack’s theorem for elliptic differential equations.  Comm. Pure Appl. Math., 14:577–591, 1961.
  * [Mos64] Jürgen Moser.  A Harnack inequality for parabolic differential equations.  Comm. Pure Appl. Math., 17:101–134, 1964.
  * [Mos67] Jürgen Moser.  Correction to: “A Harnack inequality for parabolic differential equations”.  Comm. Pure Appl. Math., 20:231–236, 1967.
  * [Mos71] J. Moser.  On a pointwise estimate for parabolic differential equations.  Comm. Pure Appl. Math., 24:727–740, 1971.
  * [Mou19] Mihalis Mourgoglou.  Regularity theory and Green’s function for elliptic equations with lower order terms in unbounded domains.  arXiv preprint arXiv:1904.04722, 2019.
  * [MS04] Pierre D. Milman and Yu. A. Semenov.  Global heat kernel bounds via desingularizing weights.  J. Funct. Anal., 212(2):373–398, 2004.
  * [Nas58] John Nash.  Continuity of solutions of parabolic and elliptic equations.  Amer. J. math, 80(4):931–954, 1958.
  * [NUt11] A. I. Nazarov and N. N. Ural’ tseva.  The Harnack inequality and related properties of solutions of elliptic and parabolic equations with divergence-free lower-order coefficients.  Algebra i Analiz, 23(1):136–168, 2011.
  * [Osa87] Hirofumi Osada.  Diffusion processes with generators of generalized divergence form.  J. Math. Kyoto Univ., 27(4):597–619, 1987.
  * [QX19a] Zhongmin Qian and Guangyu Xi.  Parabolic equations with divergence-free drift in space $L^{\ell}_{t}L^{q}_{x}$.  Indiana Univ. Math. J., 68(3):761–797, 2019.
  * [QX19b] Zhongmin Qian and Guangyu Xi.  Parabolic equations with singular divergence-free drift vector fields.  J. Lond. Math. Soc. (2), 100(1):17–40, 2019.
  * [Sak20] Georgios Sakellaris.  Scale invariant regularity estimates for second order elliptic equations with lower order coefficients in optimal spaces.  arXiv preprint arXiv:2005.14086, 2020.
  * [Sem06] Yu. A. Semenov.  Regularity theorems for parabolic equations.  J. Funct. Anal., 231(2):375–417, 2006.
  * [SSvZ12] Gregory Seregin, Luis Silvestre, Vladimír Šverák, and Andrej Zlatoš.  On divergence-free drifts.  J. Differential Equations, 252(1):505–540, 2012.
  * [Str88] Michael Struwe.  On partial regularity results for the Navier-Stokes equations.  Comm. Pure Appl. Math., 41(4):437–458, 1988.
  * [SVZ13] Luis Silvestre, Vlad Vicol, and Andrej Zlatoš.  On the loss of continuity for super-critical drift-diffusion equations.  Archive for Rational Mechanics and Analysis, 207(3):845–877, 2013\.
  * [Tao13] Terence Tao.  Localisation and compactness properties of the Navier-Stokes global regularity problem.  Anal. PDE, 6(1):25–107, 2013.
  * [Tru71] Neil S. Trudinger.  On the regularity of generalized solutions of linear, non-uniformly elliptic equations.  Arch. Rational Mech. Anal., 42:50–62, 1971.
  * [Wu21] Bian Wu.  On supercritical divergence-free drifts.  arXiv preprint arXiv:2106.02408, 2021.
  * [Zha04] Qi S. Zhang.  A strong regularity result for parabolic equations.  Communications in Mathematical Physics, 244(2):245–260, 2004.
  * [Zha20] Xicheng Zhang.  Maximum principle for non-uniformly parabolic equations and applications.  arXiv preprint arXiv:2012.05026, 2020.
  * [Zhi04] V. V. Zhikov.  Remarks on the uniqueness of the solution of the Dirichlet problem for a second-order elliptic equation with lower order terms.  Funktsional. Anal. i Prilozhen., 38(3):15–28, 2004.