#  Gradient estimate for solutions of second-order elliptic equations

Vladimir Maz’ya  Department of Mathematics, Linköping University, SE-581 83
Linköping, Sweden and RUDN, 6 Miklukho-Maklay St, Moscow, 117198, Russia
<EMAIL_ADDRESS>and  Robert McOwen  Department of Mathematics,
Northeastern University, Boston, MA 02115<EMAIL_ADDRESS>
(Date: November 22, 2021)

###### Abstract.

We obtain a local estimate for the gradient of solutions to a second-order
elliptic equation in divergence form with bounded measurable coefficients that
are square-Dini continuous at the single point $x=0$. In particular, we treat
the case of solutions that are not Lipschitz continuous at $x=0$. We show that
our estimate is sharp.

###### Key words and phrases:

Gradient estimate, weak solution, square-Dini condition

###### 1991 Mathematics Subject Classification:

35B45, 35B40, 35J15

##  1\. Introduction

We consider a uniformly elliptic equation in divergence form:

(1) $\partial_{j}(a_{ij}(x)\,\partial_{i}u)=0\quad\hbox{for $x\in B_{1}$},$

where we have used the summation convention and the coefficients
$a_{ij}=a_{ji}$ are real-valued bounded measurable functions in the ball
$B_{1}:=\\{x\in{\mathbb{R}}^{n}:|x|<1\\}$, $n\geq 2$. A weak solution $u\in
H^{1,2}(B_{1})$, i.e. $u,\nabla u\in L^{2}(B_{1})$, is known to be Hölder
continuous by De Giorgi [4] and Nash [12], but may have unbounded gradient.
Under additional regularity of the $a_{ij}$ such as Dini continuity or Dini
mean oscillation in a domain, weak solutions are known to be $C^{1}$: cf. [7]
and [5]. Gradient bounds in a domain have also been obtained for nonlinear
problems, including minimal regularity on the data: cf. [11], [8], [1], [2],
and [3]. In [9] we obtained conditions on the coefficients $a_{ij}$ just at
$x=0$ which imply that all weak solutions of (1) are Lipschitz continuous at
$x=0$. In this paper we allow the possibility of weak solutions that are not
Lipschitz continuous at $x=0$, but obtain bounds on the growth of their
gradient in mean value. We are not aware of any previous results in the
literature of this type.

Let us describe our assumptions on the coefficients $a_{ij}(x)$. These are
most conveniently stated when $a_{ij}(0)=\delta_{ij}$, which can be achieved
by a change of coordinates. We want the $a_{ij}(x)$ to be square-Dini
continuous at $x=0$, which means

(2) $\sup_{|x|=r}|a_{ij}(x)-\delta_{ij}|\leq\omega(r)\quad\hbox{as}\ r=|x|\to
0,$

where the modulus of continuity $\omega(r)$ is a nondecreasing, positive,
continuous function for $r>0$ satisfying

(3) $\int_{0}^{1}\frac{\omega^{2}(r)}{r}\,dr<\infty.$

Since $\omega(r)$ does not vanish rapidly as $r\to 0$, it is natural and often
convenient to assume

(4) $\hbox{for some $\kappa>0$},\ \omega(r)\,r^{-1+\kappa}\ \hbox{is
nonincreasing for $0<r<1$.}$

As seen in [9], the conditions (2)-(3) allow the regularity of solutions of
(1) at $x=0$ to be determined by the asymptotic behavior of solutions of the
$n$-dimensional dynamical system

(5) $\frac{d\phi}{dt}+R(e^{-t})\phi=0\quad\hbox{for}\ 0<t<\infty,$

where $R(r)$ is the matrix function

(6)
$R(r):={\diagup\hskip-11.9501pt\int}_{S^{n-1}}\left(A(r\theta)-nA(r\theta)\theta\otimes\theta\right)\,d\theta\quad\hbox{for
$0<r<1$.}$

Here $A=(a_{ij})$, $A\theta\otimes\theta$ is the outer product of the vectors
$A\theta$ and $\theta$, and $d\theta$ denotes standard surface measure on
$S^{n-1}$; the slashed integral in (6) and throughout this paper denotes mean
value. We note that111Throughout this paper we use $c$ to denote a constant
that does not depend upon the solution $u$ but whose value may change with
each occurence. $|R(r)|\leq c\,\omega(r)$, but the matrix $R$ need not be
symmetric. In [9] we found that, if the dynamical system (5) is uniformly
stable as $t\to\infty$, then a weak solution of (1) must be Lipschitz
continuous at $x=0$. In this paper, we do not assume (5) is uniformly stable,
but seek an estimate on the growth of the gradient of a weak solution.

To formulate our main result, for a square matrix $M$ let $\mu[M]$ denote the
largest eigenvalue of the symmetric matrix $(M+M^{*})/2$, where $M^{*}$
denotes the adjoint (i.e. transpose) of $M$.

######  Theorem 1.

Under the above conditions on $a_{ij}$, if $u\in H^{1,2}(B_{1})$ is a weak
solution of (1), then

(7) $\left({\diagup\hskip-11.9501pt\int}_{r<|x|<2r}|\nabla
u(x)|^{2}\,dx\right)^{1/2}\leq
c\,\|u\|_{L^{2}(B_{1})}\,\exp\left(\int_{r}^{1}\mu[-R(\rho)]\,\frac{d\rho}{\rho}\right)\quad\hbox{for}\
0<r<1/2.$

We will give an example in Section 6 that shows the estimate (7) is sharp. The
proof of (7) will only concern $0<r<r_{0}$ for $r_{0}$ sufficiently small,
since (7) is clearly true for $r_{0}<r<1$. Under additional assumptions on the
$a_{ij}(x)$ for $x\not=0$, it is possible to replace the left hand side of (7)
by a local pointwise estimate; see, for example, Theorem 2 in [10]. While our
motivation for this study was to estimate the growth of the gradient of
solutions, Theorem 1 also applies to decay. In particular, if

(8) $\int_{r}^{1}\mu[-R(\rho)]\,\frac{d\rho}{\rho}\to-\infty\quad\hbox{as}\
r\to 0,$

then the $L^{2}$-mean of the gradient of solutions must vanish as $r\to 0$. In
fact, in [9] it was shown that (8) implies the stronger conclusion that $u$ is
differentiable at $x=0$ and $\nabla u(0)=0$.

Let us make some additional comments. First of all, let us denote our
estimator by

(9) $E(r):=\exp\left[\int_{r}^{1}\mu[-R(\rho)]\,\frac{d\rho}{\rho}\right].$

If $\mu[-R(r)]\geq 0$, then $E(r)$ is nonincreasing in $r$, and we may have
$E(r)\to\infty$ as $r\to 0$, but not quickly since $\mu[-R(r)]\to 0$ as $r\to
0$. Similarly, if $\mu[-R(r)]\leq 0$, then $E(r)$ is nondecreasing in $r$, and
when (8) occurs we have $E(r)\to 0$ as $r\to 0$, but not quickly. In fact, for
any $0<\lambda<1$, we can check that $E(r)\,r^{-\lambda}$ is decreasing and
$E(r)\,r^{-\lambda}$ is increasing on $0<r<r_{0}$ for $r_{0}$ sufficiently
small. Since we could replace $B_{1}$ in (1) by $B_{r_{0}}$, we will simply
assume

(10) $E(r)\,r^{-\lambda}\ \hbox{is decreasing and}\ E(r)\,r^{\lambda}\
\hbox{is increasing for $0<r<1$.}$

Secondly, regarding the proof, we will write the weak solution $u$ in the
following form:

(11a) $u(x)=u_{0}(|x|)+\vec{v}(|x|)\cdot x+w(x)\quad\hbox{for}\ 0<|x|<1,$
where (11b)
$u_{0}(r):={\diagup\hskip-11.9501pt\int}_{S^{n-1}}u(r\theta)\,d\theta,\quad
v_{k}(r)=\frac{n}{r}{\diagup\hskip-11.9501pt\int}_{S^{n-1}}u(r\theta)\,\theta_{k}\,d\theta,\quad\theta_{k}=x_{k}/|x|,$
and (11c)
${\diagup\hskip-11.9501pt\int}_{S^{n-1}}w(r\theta)\,d\theta=0={\diagup\hskip-11.9501pt\int}_{S^{n-1}}w(r\theta)\,\theta_{k}\,d\theta\quad\hbox{for}\
k=1,\dots,n.$

We shall show that the functions $\vec{v}$ and $\vec{v}\,^{\prime}$ satisfy a
$2n$-dimensional dynamical system that depends on $w$, and that $w$ satisfies
a Poisson equation that depends on $\vec{v},\vec{v}\,^{\prime}$; this enables
us to obtain estimates on $\vec{v}$, $\vec{v}\,^{\prime}$, and $\nabla w$, and
then we are able to estimate $u_{0}^{\prime}$. In the estimates for $w$, we
will use $L^{p}$-means over annuli: for $f\in
H^{1,p}_{\ell\mathrm{o}c}({\mathbb{R}}^{n}\backslash\\{0\\})$ where $1\leq
p<\infty$, let us define

(12)
$M_{p}(f,r):=\left({\diagup\hskip-11.9501pt\int}_{A_{r}}|f(x)|^{p}\,dx\right)^{1/p},\quad\hbox{where}\
A_{r}=\\{x\in{\mathbb{R}}^{n}:r<|x|<2r\\}.$

The $f$ in (12) may be scalar or vector-valued. Using this notation, (7) may
be written

(13) $M_{2}(\nabla u,r)\leq c\,E(r)\,\|u\|_{L^{2}(B_{1})}\quad\hbox{for}\
0<r<1/2.$

##  2\. Estimates for an $n$-dimensional Dynamical System

If we let $R_{1}(t)=R(e^{-t})$ for $R(t)$ as in (6) and
$\varpi(t)=\omega(e^{-t})$, then we can write the dynamical system (5) as

(14) $\frac{d\phi}{dt}=-R_{1}(t)\phi\quad\hbox{for}\ t>0,$

where $|R_{1}(t)|\leq c\,\varpi(t)$ for $t>0$. If we take scalar product with
$\phi$ on both sides of (14), we obtain

$\frac{1}{2}\frac{d}{dt}\,|\phi|^{2}=-(R_{1}\phi,\phi)=-\frac{1}{2}([R_{1}(t)+R_{1}^{*}(t)]\phi,\phi).$

Using the variational characterization of the largest eigenvalue of a
symmetric matrix, we have

$\frac{1}{2}\frac{d}{dt}\,|\phi|^{2}\leq\mu[-R_{1}(t)]\,|\phi|^{2},$

and we can integrate this differential inequality to obtain

$|\phi(t)|\leq|\phi(0)|\exp\left(\int_{0}^{t}\mu[-R_{1}(s)]\,ds\right).$

If we define

(15)
${\mathcal{E}}(t):=E(e^{-t})=\exp\left(\int_{0}^{t}\mu[-R_{1}(s)]\,ds\right),$

then we have proved the following:

######  Proposition 1.

The fundamental matrix $\Phi$ for (14) with $\Phi(0)=I$ satisfies the
following:

(16) $|\Phi(t)|\leq{\mathcal{E}}(t)\quad\hbox{for}\ t>0$

and

(17)
$|\Phi(t)\Phi^{-1}(s)|\leq\exp\left(\int_{s}^{t}\mu[-R_{1}(\upsilon)]\,d\upsilon\right)=\frac{{\mathcal{E}}(t)}{{\mathcal{E}}(s)},\quad\hbox{for}\
0<s<t.$

Of course, the fundamental matrix $\Phi$ is useful in solving the
nonhomogeneous equation

(18) $\frac{d\phi}{dt}+R_{1}(t)\phi=f(t).$

######  Corollary 1.

If ${\mathcal{E}}^{-1}(t)\,f(t)\in L^{1}(0,\infty)$ then the solution $\phi$
of (18) satisfies

$|\phi(t)|\leq{\mathcal{E}}(t)\left(|\phi(0)|+\|{\mathcal{E}}^{-1}\,f\|_{L^{1}}\right).$

###### Proof.

This follows immediately from (17) and the variation of constants formula:

$\phi(t)=\Phi(t)\left(\phi(0)+\int_{0}^{t}\Phi^{-1}(s)\,f(s)\,ds\right).$

∎

##  3\. Estimates for a $2n$-dimensional Dynamical System

In this section we consider the $2n\times 2n$ system on $(0,\infty)$

(19a) $\frac{d}{dt}\begin{pmatrix}\phi\\\
\psi\end{pmatrix}+\begin{pmatrix}0&0\\\
0&-nI\end{pmatrix}\begin{pmatrix}\phi\\\
\psi\end{pmatrix}+{\mathcal{R}}(t)\begin{pmatrix}\phi\\\
\psi\end{pmatrix}=F(t),$ where i) ${\mathcal{R}}$ is a $2n\times 2n$ matrix of
the form (19b) ${\mathcal{R}}(t)=\begin{pmatrix}R_{1}(t)&R_{2}(t)\\\
R_{3}(t)&R_{4}(t)\end{pmatrix}\quad\hbox{with}\
\boldsymbol{\bigr{|}}R_{j}(t)\boldsymbol{\bigr{|}}\leq\varpi(t)\ \hbox{on}\
0<t<\infty,$ and ii) $F=(F_{1},F_{2})$ satisfies (19c)
${\mathcal{E}}^{-1}F_{1}\in L^{1}(0,\infty),$ and there exists $\delta>0$ so
that for any choice of $\alpha\in[n-\delta,n)$ there is a constant
$c_{\alpha}$ so that (19d) $e^{\alpha t}\int_{t}^{\infty}|F_{2}(s)|e^{-\alpha
s}\,ds\leq c_{\alpha}\,\varpi(t)\,{\mathcal{E}}(t)\quad\hbox{for}\
0<t<\infty.$

We want $\psi$ to satisfy the “finite energy condition”

(20) $\int_{0}^{\infty}(|\psi|^{2}+|\psi_{t}|^{2})e^{-nt}dt<\infty.$

######  Proposition 2.

Under the above hypotheses, all solutions $(\phi,\psi)$ of (19) that satisfy
(20) must satisfy the following estimates:

(21a) $|\phi(t)|\leq
c\,{\mathcal{E}}(t)(c_{\alpha}+|\phi(0)|+\|{\mathcal{E}}^{-1}F_{1}\|_{L^{1}})$
(21b) $|\psi(t)|\leq
c\,\varpi(t)\,{\mathcal{E}}(t)(c_{\alpha}+|\phi(0)|+\|{\mathcal{E}}^{-1}F_{1}\|_{L^{1}}).$

###### Proof.

As in the previous section, we let $\Phi$ denote the fundamental solution for
$\dot{\phi}+R_{1}\phi=0$ on $[0,\infty)$ with $\Phi(0)=I$. We need to estimate
the fundamental solution $\Psi(t)$ for $\dot{\psi}+R_{4}\psi=0$ with
$\Psi(0)=I$. Since $\varpi(t)\to 0$ as $t\to\infty$ and
$\int_{0}^{\infty}\varpi^{2}(t)\,dt<\infty$, for any $\delta\in(0,1)$ we have
$\varpi(t)<\delta$ and $\int_{t_{0}}^{\infty}<\delta$ for $t_{0}$ sufficiently
large. Without loss of generality, we will assume that

(22)
$\varpi(t)<\delta\quad\hbox{and}\quad\int_{0}^{\infty}\varpi^{2}(t)\,dt<\delta.$

Note that this $\delta$ can be made small independent of other constants that
may appear. We can now use Gronwall’s inequality to conclude

(23) $|\Psi(t)\Psi^{-1}(s)|\leq e^{\delta|t-s|}\quad\hbox{for}\ t,s>0.$

Applying variation of constants in the first equation in (19a), we can write

(24a)
$\phi(t)=\Phi(t)\phi(0)+\Phi(t)\int_{0}^{t}\Phi^{-1}(s)(F_{1}(s)-R_{2}(s)\psi(s))\,ds.$
Variation of constants in the second equation in (19a) combined with (20)
yields (24b)
$\psi(t)=e^{nt}\Psi(t)\int_{t}^{\infty}\Psi^{-1}(s)(R_{3}(s)\phi(s)-F_{2}(s))\,e^{-ns}\,ds.$

If we plug (24b) into (24a) we obtain the following equation for $\phi$:

(25a) $\phi(t)+S\phi(t)=\eta_{0}(t)+\eta_{1}(t)+\eta_{2}(t),$ where (25b)
$\displaystyle S\phi(t)$
$\displaystyle=-\Phi(t)\int_{0}^{t}\Phi^{-1}(s)\,R_{2}(s)\,e^{ns}\,\Psi(s)\int_{s}^{\infty}\Psi^{-1}(\upsilon)\,R_{3}(\upsilon)\,\phi(\upsilon)\,e^{-nu}\,d\upsilon\,ds,$
$\displaystyle\eta_{0}(t)=\Phi(t)\phi(0),\qquad\quad\eta_{1}(t)=\Phi(t)\int_{0}^{t}\Phi^{-1}(s)\,F_{1}(s)\,ds,$
$\displaystyle\eta_{2}(t)$
$\displaystyle=\Phi(t)\int_{0}^{t}\Phi^{-1}(s)\,R_{2}(s)\,e^{ns}\,\Psi(s)\int_{s}^{\infty}\Psi^{-1}(\upsilon)\,F_{2}(\upsilon)\,e^{-n\upsilon}\,d\upsilon\,ds.$

We want to show (25) has a solution $\phi$ and estimate it in terms of
$\varpi$, ${\mathcal{E}}$, and $F$.

Let $X=C([0,\infty),{\mathbb{R}}^{n})$ with the norm

(26)
$\|\phi\|_{X}:=\sup\left\\{\frac{|\phi(t)|}{{\mathcal{E}}(t)}:0<t<\infty\right\\}.$

Assume $\|\phi\|_{X}=1$, so $|\phi(t)|\leq{\mathcal{E}}(t)$. Using this, (17),
and (23), we obtain the estimate

(27) $\displaystyle|S\phi(t)|$
$\displaystyle\leq\int_{0}^{t}|\Phi(t)\Phi^{-1}(s)\,|R_{2}(s)|\,e^{(n-\delta)s}\int_{s}^{\infty}e^{(\delta-n)\upsilon}|R_{3}(\upsilon)|{\mathcal{E}}(\upsilon)\,d\upsilon\,ds$
$\displaystyle\leq
c\,{\mathcal{E}}(t)\int_{0}^{t}({\mathcal{E}}(s))^{-1}\,\varpi^{2}(s)\,e^{(n-\delta)s}\int_{s}^{\infty}e^{(\delta-n)\upsilon}{\mathcal{E}}(\upsilon)\,d\upsilon\,ds.$

Regarding the last integral, we can use integration by parts to obtain

$\int_{s}^{\infty}e^{(\delta-n)\upsilon}{\mathcal{E}}(\upsilon)\,d\upsilon=\frac{e^{(\delta-n)s}{\mathcal{E}}(s)}{n-\delta}+\frac{1}{n-\delta}\int_{s}^{\infty}e^{(\delta-n)\upsilon}{\mathcal{E}}(\upsilon)\,\mu[-R_{1}(\upsilon)]\,d\upsilon.$

But, using $|\mu[-R_{1}(t)]|\leq c\,\varpi(t)$ and (22), we can take $\delta$
small enough that

$\int_{s}^{\infty}e^{(\delta-n)\upsilon}{\mathcal{E}}(\upsilon)\,d\upsilon\leq\frac{e^{(\delta-n)s}{\mathcal{E}}(s)}{n-\delta}+\frac{1}{n-\delta}\int_{s}^{\infty}e^{(\delta-n)\upsilon}{\mathcal{E}}(\upsilon)\,d\upsilon.$

By iterating this, we obtain

(28)
$\int_{s}^{\infty}e^{(\delta-n)\upsilon}{\mathcal{E}}(\upsilon)\,d\upsilon\leq
c\,e^{(\delta-n)s}{\mathcal{E}}(s).$

Using this in (27) we obtain

$|S\phi(t)|\leq c\,{\mathcal{E}}(t)\int_{0}^{t}\varpi^{2}(s)\,ds\leq
c\,{\mathcal{E}}(t)\,\delta.$

Taking $\delta<c^{-1}$, we indeed have $\|S\|_{X\to X}<1$. To verify that
$\eta_{0},\eta_{1}\in X$ we use Proposition 1 and Corollary 1 to conclude

$|\eta_{0}(t)+\eta_{1}(t)|\leq{\mathcal{E}}(t)\left(|\phi(0)|+\|{\mathcal{E}}^{-1}F_{1}\|_{L^{1}}\right).$

To show $\eta_{2}\in X$, we use (17), (28), and (19d) with $\alpha=n-\delta$
to estimate

$\displaystyle|\eta_{2}(t)|$ $\displaystyle\leq
c\,{\mathcal{E}}(t)\int_{0}^{t}({\mathcal{E}}(s))^{-1}\,\varpi(s)\,e^{(n-\delta)s}\int_{s}^{\infty}e^{(\delta-n)\upsilon}\,F_{2}(\upsilon)\,d\upsilon\,ds$
$\displaystyle\leq
c\,{\mathcal{E}}(t)\int_{0}^{t}({\mathcal{E}}(s))^{-1}\,c_{\alpha}\,\varpi^{2}(s)\,ds<c\,c_{\alpha}\,{\mathcal{E}}(t).$

We conclude that there is a solution $\phi\in X$ of (25) satisfying (21a).
Using (24b) we find $\psi$ and hence the solution $(\phi,\psi)$ of (19). To
estimate $\psi$ we use (19c), (24b), and (21a) to obtain

$\displaystyle|\psi(t)|$ $\displaystyle\leq
e^{(n-\delta)t}\int_{t}^{\infty}e^{(\delta-n)s}(\varpi(s)|\phi(s)|+|F_{2}(s)|)\,ds$
$\displaystyle\leq
c\,\varpi(t)\,{\mathcal{E}}(t)\left(c_{\alpha}+|\phi(0)|+\|{\mathcal{E}}^{-1}F_{1}\|_{L^{1}}\right),$

which is (21b). ∎

##  4\. Potential Theory Estimates

In the proof of Theorem 1, we will encounter the equation $-\Delta w=g$ in
${\mathbb{R}}^{n}\backslash\\{0\\}$, where $g$ is a distribution with certain
orthogonality properties that we will now describe. For $f\in
L^{1}_{{\ell\mathrm{o}c}}({\mathbb{R}}^{n}\backslash\\{0\\})$, let us define

(29)
$Pf(r\theta):={\diagup\hskip-11.9501pt\int}_{S^{n-1}}f(r\phi)\,d\phi+n\theta_{k}{\diagup\hskip-11.9501pt\int}_{S^{n-1}}\phi_{k}\,f(r\phi)\,d\phi.$

Using

${\diagup\hskip-11.9501pt\int}_{S^{n-1}}\theta_{k}\theta_{\ell}\,d\theta=\frac{1}{n}\delta_{k\ell}\quad\hbox{for}\
k,\ell=1,\dots,n,$

we see that $P^{2}f=Pf$, and $Pf$ is the projection of $f$ onto the functions
on $S^{n-1}$ spanned by $1,\theta_{1},\dots,\theta_{n}$. Let us define

(30) $f(r\theta)^{\perp}:=f(r\theta)-Pf(r\theta).$

To control the growth of solutions near $x=0$ we will not only use the
$L^{p}$-mean over annuli as defined in (12), but the following to control
first and second derivatives:

(31) $\displaystyle M_{1,p}(f,r):$ $\displaystyle=rM_{p}(\nabla
w,r)+M_{p}(w,r)$ $\displaystyle M_{2,p}(f,r):$
$\displaystyle=r^{2}M_{p}(D^{2}w,r)+M_{1,p}(w,r).$

The two equations that we want to solve are

(32a) $-\Delta w=f^{\perp}\quad\hbox{in}\ {\mathbb{R}}^{n}\backslash\\{0\\}$
and (32b) $-\Delta w=[\,{\rm div}\vec{f}\,]^{\perp}\quad\hbox{in}\
{\mathbb{R}}^{n}\backslash\\{0\\}.$

In both cases, we will obtain $w$ through convolution with the fundamental
solution $\Gamma(|x|)$ for $-\Delta$. The following was proved in [9] (see
Propositions 1 and 2).

######  Proposition 3.

Let $n\geq 2$ and $1<p<\infty$.

  * a)

If $f\in L^{p}_{{\ell\mathrm{o}c}}({\mathbb{R}}^{n}\backslash\\{0\\})$
satisfies

$\int_{|x|<1}|f(x)|\,|x|^{2}\,dx<\infty\quad\hbox{and}\quad\int_{|x|>1}|f(x)|\,|x|^{-n}\,dx<\infty,$

then convolution by $\Gamma$ defines a solution $w\in
H^{2,p}_{{\ell\mathrm{o}c}}({\mathbb{R}}^{n}\backslash\\{0\\})$ of (32a)
satisfying

$M_{2,p}(w,r)\leq
c\left(r^{-n}\int_{0}^{r}M_{p}(f,\rho)\,\rho^{n+1}\,d\rho+r^{2}\int_{r}^{\infty}M_{p}(f,\rho)\,\rho^{-1}\,d\rho\right).$

  * b)

If $\vec{f}=(f_{i})$ with $f_{i}\in
L^{p}_{{\ell\mathrm{o}c}}({\mathbb{R}}^{n}\backslash\\{0\\})$ satisfies

$\int_{|x|<1}|\vec{f}(x)|\,|x|\,dx<\infty\quad\hbox{and}\quad\int_{|x|>1}|\vec{f}(x)|\,|x|^{-n-1}\,dx<\infty,$

then convolution by $\Gamma$ defines a solution $w\in
H^{1,p}_{{\ell\mathrm{o}c}}({\mathbb{R}}^{n}\backslash\\{0\\})$ of (32b)
satisfying

$M_{1,p}(w,r)\leq
c\left(r^{-n}\int_{0}^{r}M_{p}(\vec{f},\rho)\,\rho^{n}\,d\rho+r^{2}\int_{r}^{\infty}M_{p}(\vec{f},\rho)\,\rho^{-2}\,d\rho\right).$

The proof of this proposition in [9] uses an expansion of $\Gamma$ in
spherical harmonics. For example, for $|x|<|y|$, $\Gamma(|x-y|)$ can be
written as a convergent series222We note that the coefficients $a_{k,m}$ were
unfortunately omitted in the corresponding formula (29) in [9].

(33)
$\Gamma(|x-y|)=\sum_{k=0}^{\infty}\,\frac{|x|^{k}}{|y|^{n-2+k}}\sum_{m=1}^{N(k)}a_{k,m}\,\varphi_{k,m}\left(\hat{x}\right)\,\varphi_{k,m}\left(\hat{y}\right),$

where $\varphi_{k,m}$ for $m=1,\dots,N(k)$ is an orthonormal basis for the
space ${\mathcal{H}}(k)$ of spherical harmonics of degree $k$, $N(k)={\rm
dim}({\mathcal{H}}(k)$, $\hat{x}=x/|x|$, and $a_{k,m}$ are certain
coefficients.

##  5\. Proof of Theorem 1

We assume that $u\in H^{1,2}(B_{1})$ is a weak solution of (1). If we
introduce a smooth function $\chi(r)$ satisfying $\chi(r)=1$ for $0\leq r\leq
1/4$ and $\chi(r)=0$ for $r\geq 1/2$, then

$\partial_{j}(a_{ij}\partial_{i}(\chi
u))=\partial_{j}(a_{ij}\chi^{\prime}\theta_{i}u)+\chi^{\prime}\theta_{j}a_{ij}\partial_{i}u,$

where (using the definition of weak solution)

$\int_{B_{1}}\chi^{\prime}\theta_{j}\,a_{ij}\partial_{i}u\,dx=\int_{B_{1}}a_{ij}\partial_{i}u\,\partial_{j}\chi\,dx=0.$

To use the results of the previous section, we need a problem defined on
${\mathbb{R}}^{n}$, so let us extend $a_{ij}$ and $\omega(r)$ to
${\mathbb{R}}^{n}$ by defining

(34) $a_{ij}(x)=\delta_{ij}\ \hbox{for}\
|x|>1\quad\hbox{and}\quad\omega(r)=\omega(1)\ \hbox{for}\ r>1.$

Now let us consider the problem

(35a) $\partial_{j}(a_{ij}(x)\,\partial_{i}\widetilde{u})={\rm
div}\vec{f}+f_{0},\quad\hbox{for $x\in{\mathbb{R}}^{n}$},$ where (35b)
$f_{j}:=a_{ij}\chi^{\prime}\theta_{i}u\quad\hbox{and}\quad
f_{0}:=\chi^{\prime}\theta_{j}a_{ij}\partial_{i}u$ are supported in
$1/4<|x|<1/2$ and satisfy (35c) $\|\vec{f}\|_{L^{2}}+\|f_{0}\|_{L^{2}}\leq
c\,\|u\|_{H^{1,2}(B_{1/2})}\leq c\,\|u\|_{L^{2}(B_{1})},$ where we have also
used the Cacciopolli inequality, Moreover, as observed above, we have (35d)
$\int_{{\mathbb{R}}^{n}}f_{0}(x)\,dx=0.$

We want $\widetilde{u}\in H^{1,2}({\mathbb{R}}^{n})$ to be a weak solution of
(35a) in the form

(36) $\widetilde{u}(x)=u_{0}(|x|)+\vec{v}(|x|)\cdot x+w(x)\quad\hbox{for}\
x\in{\mathbb{R}}^{n}\backslash\\{0\\},$

where the $u_{0}$, $\vec{v}$ and $w$ still satisfy the conditions (11b) and
(11c), In particular, if we let $C^{\infty}_{0}({\mathbb{R}}^{n}$ denote the
smooth functions with compact support, we have

(37)
$\int_{{\mathbb{R}}^{n}}a_{ij}(x)\,\partial_{i}\widetilde{u}\,\partial_{j}\eta\,dx=\int_{{\mathbb{R}}^{n}}f_{i}\,\partial_{i}\eta\,dx-\int_{{\mathbb{R}}^{n}}f_{0}\,\eta\,dx\quad\hbox{for
all}\ \eta\in C^{\infty}_{0}({\mathbb{R}}^{n}).$

After we find $\widetilde{u}$ and estimates on $u_{0}$, $\vec{v}$, and $w$, we
will show $\widetilde{u}=\chi\,u$.

If we plug (11a) into (37) and take $\eta=\eta(r)$, we obtain (cf. (55a) in
[9]):

(38a)
$\alpha(r)\,u_{0}^{\prime}(r)+r\,\vec{\beta}(r)\cdot\vec{v}\,^{\prime}(r)+\vec{\gamma}(r)\cdot\vec{v}(r)+p[\nabla
w](r)=\widetilde{f}(r)+r^{1-n}\int_{0}^{r}\overline{f_{0}}(\rho)\rho^{n-1}\,d\rho,$
where (38b)
$\displaystyle\alpha(r)={\diagup\hskip-11.9501pt\int}_{S^{n-1}}a_{ij}(r\theta)\theta_{i}\theta_{j}\,ds_{\theta},$
$\displaystyle\qquad\beta_{k}(r)={\diagup\hskip-11.9501pt\int}_{S^{n-1}}a_{ij}(r\theta)\theta_{i}\theta_{j}\theta_{k}\,ds_{\theta},$
$\displaystyle\gamma_{j}(r)={\diagup\hskip-11.9501pt\int}_{S^{n-1}}a_{ij}(r\theta)\theta_{i}\,ds_{\theta},$
$\displaystyle\qquad p[\nabla
w](r)={\diagup\hskip-11.9501pt\int}_{S^{n-1}}a_{ij}(r\theta)\,\partial_{j}w(r\theta)\,\theta_{i}\,ds_{\theta},$
$\displaystyle\widetilde{f}(r)={\diagup\hskip-11.9501pt\int}_{S^{n-1}}f_{i}(r\theta)\theta_{i}\,d\theta,$
$\displaystyle\qquad\overline{f_{0}}(r)={\diagup\hskip-11.9501pt\int}_{S^{n-1}}f_{0}(r\theta)\,d\theta.$
Using (2), for $0<r<1$ we see that these coefficients satisfy (38c)
$\displaystyle|\alpha(r)-1|,$ $\displaystyle|\vec{\beta}(r)|,\
|\vec{\gamma}(r)|\,\leq\omega(r)\,$ $\displaystyle|p[\nabla w](r)|$
$\displaystyle\leq\omega(r)\,{\diagup\hskip-11.9501pt\int}_{S^{n-1}}|\nabla
w|\,ds,$

while for $r>1$ we have $\alpha(r)=1$ and $\beta_{k}(r)=\gamma_{k}(r)=p[\nabla
w](r)=0$. Without loss of generality, we may assume $\alpha(r)>0$ for all
$r>0$.

Similarly, we can take $\eta=\eta(r)x_{\ell}$ in (37) and obtain a second-
order liner system of ODEs. But we can use (38a) to eliminate $u_{0}$ and then
reduce the second-order system for $\vec{v}$ to a first-order system for
$(\vec{v},\vec{v}_{r})$; of course these systems also depend on $w$. The
first-order system is simplified if we change the independent variable to
$t=-\log r\in(-\infty,\infty)$. This brings the dynamical system (19) into
play. If we follow the calculations in [9], we find that the first-order
system for $(\vec{v},\vec{v}_{r})$ that depends on $w$ may be converted to the
form

(39) $\frac{d}{dt}\begin{pmatrix}\phi\\\
\psi\end{pmatrix}+\begin{pmatrix}0&0\\\
0&-nI\end{pmatrix}\begin{pmatrix}\phi\\\
\psi\end{pmatrix}+{\mathcal{R}}(t)\begin{pmatrix}\phi\\\
\psi\end{pmatrix}=F(t,\nabla w)+G(t),$

where ${\mathcal{R}}(t)=0$ for $t<0$, but for $t>0$ it is of the form (19b)
with

$R_{1}(t)={\diagup\hskip-11.9501pt\int}_{S^{n-1}}\left(A(r\theta)-nA(r\theta)\theta\otimes\theta\right)\,d\theta\,(1+O(\varpi^{2}(t))\quad\hbox{as}\
t\to\infty.$

The term $F(t,\nabla w)$ vanishes for $t<0$, but

(40) $|F(t,\nabla w)|\leq
c\,\varpi(t)\,{\diagup\hskip-11.9501pt\int}_{S^{n-1}}|\nabla
w|\,d\theta\quad\hbox{for}\ t>0,$

and the term $G(t)$ comes from $\vec{f}$ and $f_{0}$, so is supported in $\log
2\leq t\leq 2\log 2$ and we can estimate

(41) $\|{\mathcal{E}}^{-1}G\|_{L^{1}}\leq
c\,(\|\vec{f}\|_{L^{2}}+\|f_{0}\|_{L^{2}}).$

Moreover, the relationship between $(\phi,\psi)$ and $(\vec{v},\vec{v}_{t})$
satisfies:

(42) $\begin{pmatrix}\phi\\\
\psi\end{pmatrix}-\frac{1}{n^{2}}\begin{pmatrix}n\vec{v}-\vec{v}_{t}\\\
\vec{v}_{t}\end{pmatrix}\leq
c\,\varpi(t)\left(|\vec{v}(t)|+|\vec{v}_{t}(t)|+{\diagup\hskip-11.9501pt\int}_{S^{n-1}}|\nabla
w|d\theta\right).$

Now, given $w$ with certain properties, we solve (39) with initial conditions
$\phi(0)=0=\psi(0)$ to find $\phi,\psi$ and hence $\vec{v}$, $\vec{v}_{r}$.
However, we need to separately control the dependence of $\vec{v}$ on $w$, so
let us write $\vec{v}=\vec{v}\,^{w}+\vec{v}\,^{\circ}$ where $\vec{v}\,^{w}$
is the solution with $G=0$ and $\vec{v}\,^{\circ}$ is the solution with $F=0$.
Similarly, when we use (38a) to solve for $u_{0}$ in terms of $\vec{v}$ and
$w$, we want to write $u_{0}=u_{0}^{w}+u_{0}^{\circ}$ where $u_{0}^{w}$
depends on $w$ (including $\vec{v}\,^{w}$), while $u_{0}^{\circ}$ is
independent of $w$:

$\displaystyle u^{\prime}_{0}$
$\displaystyle=(u_{0}^{w})^{\prime}+(u_{0}^{\circ})^{\prime}=\alpha^{-1}\left(-r\vec{\beta}\cdot(\vec{v}\,^{w})^{\prime}-\vec{\gamma}\cdot\vec{v}\,^{w}-p[\nabla
w]\right)$
$\displaystyle+\alpha^{-1}\left(\widetilde{f}(r)+r^{1-n}\int_{0}^{r}\overline{f_{0}}(\rho)\rho^{n-1}\,d\rho-r\vec{\beta}\cdot(\vec{v}\,^{\circ})^{\prime}-\vec{\gamma}\cdot\vec{v}\,^{\circ}\right)$

Finally, let us find the PDE that $w$ must satisfy. As in [9], let us
introduce the matrix

(43) $\Omega_{ij}(x)=a_{ij}(x)-\delta_{ij},$

which satisfies $|\Omega_{ij}(x)|\leq\omega(|x|)$ for $|x|<1$ and
$\Omega(x)=0$ for $|x|>1$. Recalling that $Pw=0$ and the notation (30), we
find that $w$ satisfies

$\Delta w+[{\rm div}(\Omega\nabla w)]^{\perp}+[{\rm
div}(\Omega\nabla(\vec{v}\cdot x))]^{\perp}+[{\rm
div}(\Omega\nabla(u_{0})]^{\perp}=[\partial_{i}f_{i}+f_{0}]^{\perp}.$

But we want to write this as

(44) $\displaystyle\Delta w$ $\displaystyle+[{\rm div}(\Omega\nabla
w)]^{\perp}+[{\rm div}(\Omega\nabla(\vec{v}\,^{w}\cdot x))]^{\perp}+[{\rm
div}(\Omega\nabla(u^{w}_{0})]^{\perp}=$
$\displaystyle[\partial_{i}f_{i}+f_{0}-{\rm
div}(\Omega\nabla(\vec{v}\,^{\circ}\cdot x))-{\rm
div}(\Omega\nabla(u_{0}^{\circ})]^{\perp}.$

If we apply $\Delta^{-1}$ (i.e. convolution by the fundamental solution) to
both sides of (44), we obtain

(45a) $w+T_{1}(w)+T_{2}(w)+T_{3}(w)=\xi,$ where (45b)
$T_{1}(w)=\Delta^{-1}\left[{\rm div}(\Omega\nabla w)\right]^{\perp}$ (45c)
$T_{2}(w)=\Delta^{-1}\left[{\rm div}(\Omega\nabla(\vec{v}^{\,w}\cdot
x))\right]^{\perp}$ (45d) $T_{3}(w)=\Delta^{-1}\left[{\rm
div}(\Omega\nabla(u_{0}^{w}))\right]^{\perp}$ (45e)
$\xi=\Delta^{-1}\left[\partial_{i}f_{i}+f_{0}-{\rm
div}(\Omega\nabla(\vec{v}\,^{\circ}\cdot x))-{\rm
div}(\Omega\nabla(u^{\circ}_{0})\right]^{\perp}.$

We want to show that (45) admits a solution $w\in Y$, where $Y$ denotes those
functions in $y\in H^{1,2}({\mathbb{R}}^{n}\backslash\\{0\\})$ for which the
following norm is finite:

(46)
$\|y\|_{Y}=\sup_{0<r<1}\frac{M_{1,2}(y,r)}{\omega(r)\,r\,E(r)}+\sup_{r>1}\frac{M_{1,2}(y,r)}{\delta\,r^{-n}}.$

Here $\delta>0$ will be small, but we may still assume $\varpi(t)$ satisfies
(22). For $\omega(r)$ this means $\omega(r)<\delta$ and
$\int_{0}^{1}\omega^{2}(r)\,r^{-1}\,dr<\delta$ for $0<r<1$.

We need to show $\xi\in Y$ and that each of the mappings $T_{j}:Y\to Y$ is
small. For future reference, we note that $y\in Y$ implies

(47) $M_{2}(\nabla y,r)\leq\begin{cases}\omega(r)\,E(r)\,\|y\|_{Y}&0<r<1\\\
\delta\,r^{-n-1}\,\|y\|_{Y}&r>1.\end{cases}$

First let us confirm that $T_{1}:Y\to Y$ with small norm. Let $y\in Y$ with
$\|y\|_{Y}\leq 1$. Use Proposition 3 with $\vec{f}=\Omega\nabla y$ to conclude

$M_{1,2}(T_{1}(y),r)\leq c\left(r^{-n}\int_{0}^{r}M_{2}(\Omega\nabla
y,\rho)\,\rho^{n}\,d\rho+r^{2}\int_{r}^{\infty}M_{2}(\Omega\nabla
y,\rho)\,\rho^{-2}\,d\rho\right).$

For any function $y\in
H^{1,2}_{{\ell\mathrm{o}c}}({\mathbb{R}}^{n}\backslash\\{0\\})$, note that
$M_{2}(\Omega\nabla y,r)=0$ for $r>1$ and $M_{2}(\Omega\nabla
y,r)\leq\omega(r)\,M_{2}(\nabla y,r)$ for $0<r<1$; but for $y\in Y$ with
$\|y\|_{Y}\leq 1$ this last estimate becomes $M_{2}(\Omega\nabla
y,r)\leq\omega^{2}(r)\,E(r)$. Consequently, for $0<r<1$ we have

(48) $M_{1,2}(T_{1}(y),r)\leq
c\left(r^{-n}\int_{0}^{r}\omega^{2}(\rho)\,E(\rho)\rho^{n}\,d\rho+r^{2}\int_{r}^{1}\omega^{2}(\rho)\,E(\rho)\rho^{-2}\,d\rho\right).$

Similar to (28), we can integrate by parts to show

(49) $r^{-n}\int_{0}^{r}E(\rho)\rho^{n}\,d\rho\leq c\,r\,E(r).$

We can use this, the monotonicity of $\omega$, and (22) to estimate the first
integral in (48):

$r^{-n}\int_{0}^{r}\omega^{2}(\rho)\,E(\rho)\rho^{n}\,d\rho\leq
c\,\delta\,\omega(r)\,r\,E(r).$

We can also use (4) and (10) to estimate the second integral in (48):

$r^{2}\int_{r}^{1}\omega^{2}(\rho)\,E(\rho)\rho^{-2}\,d\rho\leq\delta\,r^{2}\omega(r)r^{-1+\kappa}E(r)\,r^{-\lambda}\int_{r}^{1}\rho^{-1-\kappa-\lambda}\,d\rho,$

where we have chosen $0<\lambda<\kappa$. We can evaluate this last integral to
obtain

$r^{2}\int_{r}^{1}\omega^{2}(\rho)\,E(\rho)\rho^{-2}\,d\rho\leq\frac{\delta}{\kappa-\lambda}\,\omega(r)\,r\,E(r)(1+r^{\kappa-\lambda})\leq
c\,\delta\,\omega(r)\,r\,E(r).$

Meanwhile for $r>1$, we have

$M_{1,2}(T_{1}(y),r)\leq c\,r^{-n}\int_{0}^{1}M_{2}(\Omega\nabla
y,\rho)\rho^{n}\,d\rho\leq
c\,\delta^{2}\,r^{-n}\int_{0}^{1}E(\rho)\rho^{n}\,d\rho\leq
c\,\delta^{2}\,r^{-n}.$

These estimates together show $\|T_{1}(y)\|_{Y}\leq c\,\delta$. Taking
$\delta$ small enough, we have $\|T_{1}\|_{Y\to Y}<1/3$.

For $T_{2}$ we need to use Proposition 2 to obtain estimates on
$\vec{v}\,^{y}$ for $y\in Y$, so let us first confirm that $F(t,\nabla y)$ in
(39) satisfies the conditions (19c) and (19d). For (19c) we use (40), to
compute

$\displaystyle\int_{0}^{\infty}{\mathcal{E}}^{-1}(t)\,|F_{1}(t,\nabla y)|dt$
$\displaystyle\leq\int_{0}^{\infty}{\mathcal{E}}^{-1}(t)\,\varpi(t)\,{\diagup\hskip-11.9501pt\int}|\nabla
y|\,d\theta\,dt$
$\displaystyle=c\int_{0}^{1}E^{-1}(r)\,\omega(r)\,{\diagup\hskip-11.9501pt\int}|\nabla
y|\,d\theta\,\frac{dr}{r}$
$\displaystyle=c\sum_{j=0}^{\infty}\int_{2^{-j-1}}^{2^{-j}}E^{-1}(r)\,\omega(r)\,{\diagup\hskip-11.9501pt\int}|\nabla
y|\,d\theta\,\frac{dr}{r}$

But (4) and (10) imply that for $2^{-j-1}<r<2^{-j}$ we have
$E^{-1}(r)<2^{\lambda}\,E^{-1}(2^{-j-1})$ and
$\omega(r)<2^{1-\kappa}\omega(2^{-j-1})$, so

$\displaystyle\int_{0}^{\infty}{\mathcal{E}}^{-1}(t)\,|F_{1}(t,\nabla y)|dt$
$\displaystyle\leq
c\,\sum_{j=0}^{\infty}E^{-1}(2^{-j-1})\,\omega(2^{-j-1})\,M_{1}(\nabla
y,2^{-j-1})$ $\displaystyle\leq
c\,\sum_{j=0}^{\infty}E^{-1}(2^{-j-1})\,\omega(2^{-j-1})\,M_{2}(\nabla
y,2^{-j-1})$ $\displaystyle\leq
c\,\sum_{j=0}^{\infty}\omega^{2}(2^{-j-1})\,\|y\|_{Y}\leq
c\,\left(\int_{0}^{1}\frac{\omega^{2}(r)}{r}dr\right)\,\|y\|_{Y}.$

Hence we have confirmed (19c) with

(50) $\|{\mathcal{E}}^{-1}(t)\,F_{1}(t,\nabla y)\|_{L^{1}(0,\infty)}\leq
c\,\delta\,\|y\|_{Y}.$

To confirm (19d) and estimate $c_{\alpha}$ for $\alpha\in[n-\delta,n)$, we use
the same tools to compute

$\displaystyle e^{\alpha t}\int_{t}^{\infty}|F(t,\nabla y)|e^{-\alpha s}\,ds$
$\displaystyle=r^{-\alpha}\int_{0}^{r}\omega(\rho){\diagup\hskip-11.9501pt\int}|\nabla
y|d\theta\,\rho^{\alpha}\,\frac{d\rho}{\rho}$ $\displaystyle\leq
c\,r^{-\alpha}\,\omega(r)\sum_{j=0}^{\infty}\int_{2^{-j-1}r}^{2^{-j}r}{\diagup\hskip-11.9501pt\int}|\nabla
y|\,d\theta\,\rho^{\alpha}\,\frac{d\rho}{\rho}$ $\displaystyle\leq
c\,r^{-\alpha}\,\omega(r)\sum_{j=0}^{\infty}M_{1}(\nabla
y,2^{-j-1}r)\,(2^{-j}r)^{\alpha}$ $\displaystyle\leq
c\,r^{-\alpha}\,\omega(r)\sum_{j=0}^{\infty}M_{2}(\nabla
y,2^{-j-1}r)\,(2^{-j}r)^{\alpha}$ $\displaystyle\leq
c\,r^{-\alpha}\,\omega(r)\,\|y\|_{Y}\,\sum_{j=0}^{\infty}\omega(2^{-j-1}r)\,E(2^{-j-1}r)\,(2^{-j-1}r)^{\alpha}$
$\displaystyle\leq
c\,r^{-\alpha}\,\omega(r)\,\|y\|_{Y}\int_{0}^{r}\omega(\rho)\,E(\rho)\,\rho^{\alpha}\,\frac{d\rho}{\rho}$
$\displaystyle\leq
c\,r^{-\alpha}\,\omega^{2}(r)\,E(r)\,r^{\alpha/2}\,\|y\|_{Y}\int_{0}^{r}\rho^{\alpha/2-1}\,{d\rho}$
$\displaystyle=c\,\omega^{2}(r)\,{E}(r)\,\|y\|_{Y}=c\,\varpi^{2}(t)\,{\mathcal{E}}(t)\,\|y\|_{Y}.$

We conclude that (19d) holds with

(51) $c_{\alpha}=c\,\delta\,\|y\|_{Y}.$

Now let us confirm $T_{2}:Y\to Y$ with small norm. We want to use Proposition
3 with $\vec{f}=\Omega\nabla(\vec{v}\,^{y}\cdot x)$, so we need to estimate
$M_{2}(\nabla(\vec{v}\,^{y}\cdot x))$ for $|x|<1$. Using (42) and then
Proposition 2 with (50) and (51) we have

$\displaystyle|\nabla(\vec{v}\,^{y}\cdot x))|$ $\displaystyle\leq
c(r|(\vec{v}\,^{y})^{\prime}|+|\vec{v}\,^{y}|)\leq
c\,(1+\varpi(t))\,(|\phi^{y}(t)|+|\psi^{y}(t)|)+c\,\varpi(t){\diagup\hskip-11.9501pt\int}|\nabla
y|\,d\theta$ $\displaystyle\leq
c\,\delta\,{\mathcal{E}}(t)\,\|y\|_{Y}=c\,\delta\,E(r)\,\|y\|_{Y}.$

Assuming $\|y\|_{Y}\leq 1$ and applying Proposition 3, we have for $0<r<1$

$\displaystyle M_{1,2}(T_{2}[y],r)$ $\displaystyle\leq
c\,\delta\left(r^{-n}\int_{0}^{r}\omega(\rho)E(\rho)\rho^{n}\,d\rho+r^{2}\int_{r}^{1}\omega(\rho)E(\rho)\,\rho^{-2}\,d\rho\right)$
$\displaystyle\leq c\,\delta\,\omega(r)\,r\,E(r)$

where we have used (49), (4), and (10) as we did for $T_{1}$. Meanwhile, for
$r>1$ we have

$M_{1,2}(T_{2}[y],r)\leq
c\,\delta^{2}\,r^{-n}\int_{0}^{1}E(\rho)\,\rho^{n}\,d\rho\leq
c\,\delta^{2}\,r^{-n}.$

So $\|T_{2}[y]\|_{Y}<c\,\delta$ and by taking $\delta$ small enough we can
arrange $\|T_{2}\|_{Y\to Y}<1/3$.

The last map to consider is $T_{3}$, so we need to estimate $M_{2}(\nabla
u_{0}^{y},r)$ for $0<r<1$. Recall that $u_{0}^{y}$ comes from (38a) without
$\vec{f}$ and $f_{0}$, so

$|\nabla u_{0}^{y}(r)|\leq
c\left(\omega(r)\,r\,|(\vec{v}\,^{y})^{\prime}(r)|+\omega(r)|\vec{v}\,^{y}(r)|+p[\nabla
w](r)\right)\leq c\,\delta\,\omega(r)\,E(r)\,\|y\|_{Y}.$

The rest of the proof that $\|T_{3}\|_{Y\to Y}<1/3$ is the same as for
$T_{2}$.

Finally we need to confirm that $\xi\in Y$. We take the terms one at a time.
For $\xi_{1}=\Delta^{-1}[{\rm div}\vec{f}]^{\perp}$ where $\vec{f}\in
L^{2}({\mathbb{R}}^{n})$ has support in $1/4<|x|<1/2$, we apply Proposition 3b
to conclude

$M_{1,2}(\xi_{1},r)\leq\begin{cases}c\,\|\vec{f}\|_{L^{2}}\,r^{2}&0<r<1\\\
c\,\|\vec{f}\|_{L^{2}}\,r^{-n}&r>1.\end{cases}$

Similarly, for $\xi_{2}=\Delta^{-1}[f_{0}]^{\perp}$ where $f_{0}\in
L^{2}({\mathbb{R}}^{n})$ has support in $1/4<|x|<1/2$, we apply Proposition 3a
to conclude

$M_{1,2}(\xi_{2},r)\leq
M_{2,2}(\xi_{2},r)\leq\begin{cases}c\,\|f_{0}\|_{L^{2}}\,r^{2}&0<r<1\\\
c\,\|f_{0}\|_{L^{2}}\,r^{-n}&r>1.\end{cases}$

For $\xi_{3}=\Delta^{-1}[{\rm div}(\Omega\nabla(\vec{v}\,^{\circ}\cdot
x)]^{\perp}$, we need to estimate $M_{2}(\nabla(\vec{v}\,^{\circ}\cdot x),r)$
for $0<r<1$. But by Proposition 2 applied to (39) with $F=0$ and $G$
satisfying (41), we know that

$|\nabla(\vec{v}\,^{\circ}\cdot x)|\leq c\,(|\phi(t)|+|\psi(t)|)\leq
c\,(\|\vec{f}\|_{L^{2}}+\|f_{0}\|_{L^{2}}).$

This enables us to estimate (using tricks as above)

$\displaystyle M_{1,2}(\xi_{3},r)$ $\displaystyle\leq
c\,(\|\vec{f}\|_{L^{2}}+\|f_{0}\|_{L^{2}})\left(r^{-n}\int_{0}^{r}\omega(\rho)\,E(\rho)\,\rho^{n}\,d\rho+r^{2}\int_{r}^{1}\omega(\rho)\,E(\rho)\,\rho^{-2}\,d\rho\right)$
$\displaystyle\leq
c\,\begin{cases}(\|\vec{f}\|_{L^{2}}+\|f_{0}\|_{L^{2}})\,\omega(r)\,r\,E(r)&0<r<1,\\\
(\|\vec{f}\|_{L^{2}}+\|f_{0}\|_{L^{2}})\,r^{-n}&r>1.\end{cases}$

Finally, for $\xi_{4}=\Delta^{-1}[{\rm
div}(\Omega\nabla(u^{\circ}_{0})]^{\perp}$ we need to estimate
$M_{2}(\nabla(u^{\circ}_{0}),r)$ for $0<r<1$. But from (38a) we know

$(u_{0}^{\circ})^{\prime}(r)=\alpha^{-1}(r)\left[\widetilde{f}(r)+r^{1-n}\int_{0}^{r}\overline{f_{0}}(\rho)\,d\rho\right],$

which has support in $1/4<r<1/2$, so $M_{1,2}(\xi_{4},r)$ satisfies the same
estimates as for $M_{1,2}(\xi_{1},r)$ and $M_{1,2}(\xi_{1},r)$. We summarize
this by

(52) $\|\xi\|_{Y}\leq c\,(\|\vec{f}\|_{L^{2}}+\|f_{0}\|_{L^{2}})\leq
c\,\|u\|_{L^{2}(B_{1})}.$

We have shown that there exists $w\in Y$ with $\|w\|_{Y}\leq
c\,\|u\|_{L^{2}(B_{1})}$ such that the function $\widetilde{u}(x)$ as in (36)
is a weak solution of (35a). Moreover, $\widetilde{u}$ satisfies the desired
estimate (7) since $M_{2}(u^{\prime}_{0},r)$, $M_{2}(r\,\vec{v},^{\prime}r)$,
$M_{2}(\nabla w,r)\leq c\,\omega(r)\,E(r)\,\|u\|_{L^{2}(B_{1})}$ and
$M_{2}(\vec{v},r)\leq c\,E(r)\,\|u\|_{L^{2}(B_{1})}$. It only remains to show
that $\widetilde{u}=\chi u$. But $z=\widetilde{u}-\chi u$ is a solution of
$\partial_{j}(a_{ij}\partial_{i}u)=0$ in ${\mathbb{R}}^{n}$ that vanishes at
infinity, so $z=0$ by the maximum principle. $\Box$

##  6\. The Gilbarg-Serrin Example

In [6], Gilbarg and Serrin considered coefficients of the form

(53) $a_{ij}(r,\theta)=\delta_{ij}+g(r)\,\theta_{i}\theta_{j},\quad\hbox{where
$r=|x|$ and $\theta_{i}=x_{i}/r$.}$

They assumed $g(r)$ is a bounded function satisfying $g(r)>-1+\varepsilon$,
which guarantees uniform ellipticity. In [9] it was shown that the system (14)
reduces to the scalar equation

(54) $\frac{d\phi}{dt}=\frac{n-1}{n}\,{\rm g}(t)\,\phi.$

There is no matrix, so the unique “eigenvalue” of $-R(r)$ is
$\frac{n-1}{n}\,g(r)$ and

$E(r)=\exp\left[\frac{n-1}{n}\int_{r}^{1}\frac{g(\rho)}{\rho}\,d\rho\right].$

In [10] it was shown that, provided $g(r)$ has finite variation on $(0,1)$,
then there is a solution of (1) of the form (11a) with $u_{0}=0$,

(55)
$\vec{v}(r)=c\,\vec{e}_{j}\exp\left[\frac{n-1}{n}\int_{r}^{1}\frac{g(\rho)}{\rho}\,d\rho\right]\left(1+o(1))\right)\quad\hbox{as}\
r\to 0,$

and $w$ satisfying $M_{2}(\nabla w,r)\leq c\,\omega(r)\,{E}(r)$. Here
$\vec{e}_{j}$ is any one of the $n$ coordinate unit vectors, so we have
$|\vec{v}|=c\,E(r)(1+o(1))$ and we see that the estimate (7) in Theorem 1 is
sharp in this case.

Acknowledgement: This paper has been supported by the RUDN University
Strategic Academic Leadership Program.

## References

  * [1] A. Cianchi, V. Maz’ya, Global boundedness of the gradient for a class of nonlinear elliptic systems, Arch. Ration. Mech. Anal., 212 (2014), p.129-177.
  * [2] A. Cianchi, V. Maz’ya, Global gradient estimates in elliptic problems under minimal data and domain regularity, Commun. Pure Appl. Anal. 14 (2015), no. 1, 285–311.
  * [3] C. De Filippis, G. Mingione, Lipschitz bounds and nonautonomous integrals, Arch. Rat. Mech. Anal, doi: 10.1007/s00205-021-01698-5
  * [4] E. De Giorgi, Sulla differenziabilità e l’analiticità delle estremali degil integrali multipli regolari, Mem. Accad. Sci. Torino Cl. Sci. Fis. Mat. Natur. (3) 3 (1957), 25-43.
  * [5] H. Dong, S. Kim, On $C^{1}$, $C^{2}$, and weak type $(1,1)$ estimates for linear elliptic operators, Comm. Partial Differential Equations, 42 (2017), 417-435.
  * [6] D. Gilbarg, J. Serrin, On isolated singularities of solutions of second-order elliptic equations, J. Analyse Math.  4 (1955/56), 309-340.
  * [7] P. Hartman, A. Wintner, On uniform Dini conditions in the theory of linear partial differential equations of elliptic type, American J. Math. 77 (1955), 329-354.
  * [8] T. Kuusi, G. Mingione, Universal potential estimates, J. Funct. Anal. 262 (2012), 4205-4269.
  * [9] V.G. Maz’ya, R.McOwen, Differentiablilty of solutions to second-order elliptic equations via dynamical systems, J. Differential Equations, 250 (2010), 1137-1168.
  * [10] V.G. Maz’ya, R. McOwen, Gilbarg-Serrin equation and Lipschitz regularity, arXiv:2108.05764.
  * [11] G. Mingione, Gradient potential estimates, J. Eur. Math. Soc. 13 (2011), 459-486.
  * [12] J. Nash, Continuity of solutions of parabolic and elliptic equations, Amer. J. Math. 80 (1958), 931-954.


$\mathbf{j}\in\mathbf{Z}^{2}_{\neq 0}$, and using redundancies created by the
reality condition $\overline{q}=q$, we can express (5.97) as the family of
equations

(5.98)
$\displaystyle\begin{cases}\dot{a}_{\mathbf{j}}=\displaystyle{\sum_{\mathbf{j}+\mathbf{k}-\mathbf{l}=0}\theta_{\mathbf{k}\mathbf{l}}(a_{\mathbf{k}}a_{\mathbf{l}}+b_{\mathbf{k}}b_{\mathbf{l}})+\sum_{\mathbf{j}-\mathbf{k}-\mathbf{l}=0}\theta_{\mathbf{k}\mathbf{l}}(b_{\mathbf{k}}b_{\mathbf{l}}-a_{\mathbf{k}}a_{\mathbf{l}})-\mathbf{1}_{\mathscr{D}}(\mathbf{j})\lambda_{\mathbf{j}}a_{\mathbf{j}}+\beta_{a_{\mathbf{j}}}}\\\
\dot{b}_{\mathbf{j}}=\displaystyle{\sum_{\mathbf{j}+\mathbf{k}-\mathbf{l}=0}\theta_{\mathbf{k}\mathbf{l}}(a_{\mathbf{k}}b_{\mathbf{l}}-b_{\mathbf{k}}a_{\mathbf{l}})-\sum_{\mathbf{j}-\mathbf{k}-\mathbf{l}=0}\theta_{\mathbf{k}\mathbf{l}}(a_{\mathbf{k}}b_{\mathbf{l}}+b_{\mathbf{k}}a_{\mathbf{l}})-\mathbf{1}_{\mathscr{D}}(j)\lambda_{\mathbf{j}}b_{\mathbf{j}}+\beta_{b_{\mathbf{j}}}}\end{cases}$

where all indices $\mathbf{j}$, $\mathbf{k}$, $\mathbf{l}$ belong to the first
quadrant

(5.99)
$\displaystyle\mathbf{Z}_{+}^{2}=\\{\mathbf{j}=(j_{1},j_{2})\in\mathbf{Z}^{2}_{\neq
0}\,:\,j_{2}>0\\}\cup\\{\mathbf{j}=(j_{1},j_{2})\in\mathbf{Z}_{\neq
0}^{2}\,:\,j_{2}=0\text{ and }j_{1}>0\\}.$

Ignoring the dissipative and forced terms in equation (5.98), we note that for
any triple $\mathbf{j},\mathbf{k},\mathbf{l}\in\mathbf{Z}^{2}_{+}$ with
$\mathbf{j}+\mathbf{k}-\mathbf{l}=0$ (which is the same as
$\mathbf{l}-\mathbf{j}-\mathbf{k}=0$), we can break apart the nonlinear term
in (5.98) that has only the these indices into four groups of three equations
given by:

(5.100)
$\displaystyle\begin{cases}\dot{a}_{\mathbf{j}}=\theta_{\mathbf{k}\mathbf{l}}a_{\mathbf{k}}a_{\mathbf{l}}\\\
\dot{a}_{\mathbf{k}}=\theta_{\mathbf{j}\mathbf{l}}a_{\mathbf{j}}a_{\mathbf{l}}\\\
\dot{a}_{\mathbf{l}}=-\theta_{\mathbf{j}\mathbf{k}}a_{\mathbf{j}}a_{\mathbf{k}}\end{cases},\,\,\,\begin{cases}\dot{a}_{\mathbf{j}}=\theta_{\mathbf{k}\mathbf{l}}b_{\mathbf{k}}b_{\mathbf{l}}\\\
\dot{b}_{\mathbf{k}}=\theta_{\mathbf{j}\mathbf{l}}a_{\mathbf{j}}b_{\mathbf{l}}\\\
\dot{b}_{\mathbf{l}}=-\theta_{\mathbf{j}\mathbf{k}}a_{\mathbf{j}}b_{\mathbf{k}}\end{cases},\,\,\,\begin{cases}\dot{b}_{\mathbf{j}}=\theta_{\mathbf{k}\mathbf{l}}a_{\mathbf{k}}b_{\mathbf{l}}\\\
\dot{a}_{\mathbf{k}}=\theta_{\mathbf{j}\mathbf{l}}b_{\mathbf{j}}b_{\mathbf{l}}\\\
\dot{b}_{\mathbf{l}}=-\theta_{\mathbf{j}\mathbf{k}}b_{\mathbf{j}}a_{\mathbf{k}}\end{cases},\,\,\,\begin{cases}\dot{b}_{\mathbf{j}}=-\theta_{\mathbf{k}\mathbf{l}}b_{\mathbf{k}}a_{\mathbf{l}}\\\
\dot{b}_{\mathbf{k}}=-\theta_{\mathbf{j}\mathbf{l}}b_{\mathbf{j}}a_{\mathbf{l}}\\\
\dot{a}_{\mathbf{l}}=\theta_{\mathbf{j}\mathbf{k}}b_{\mathbf{j}}b_{\mathbf{k}}.\end{cases}$

Importantly, each of the systems above conserves the _relative enstrophy_ and
_relative energy_ ; that is, for example, for the first equation in (5.100),
the _relative enstrophy_ is

(5.101)
$\displaystyle\mathcal{E}(a_{\mathbf{j}},a_{\mathbf{k}},a_{\mathbf{l}}):=a_{\mathbf{j}}^{2}+a_{\mathbf{k}}^{2}+a_{\mathbf{l}}^{2},$

which is conserved, and the _relative energy_ is

(5.102) $\displaystyle
E_{\mathbf{j}\mathbf{k}\mathbf{l}}(a_{\mathbf{j}},a_{\mathbf{k}},a_{\mathbf{l}}):=\frac{a_{\mathbf{j}}^{2}}{|\mathbf{j}|^{2}}+\frac{a_{\mathbf{k}}^{2}}{|\mathbf{k}|^{2}}+\frac{a_{\mathbf{l}}^{2}}{|\mathbf{l}|^{2}},$

which is also conserved. We let
$V_{a_{\mathbf{j}}a_{\mathbf{k}}a_{\mathbf{l}}}$,
$V_{a_{\mathbf{j}}b_{\mathbf{k}}b_{\mathbf{l}}}$,
$V_{b_{\mathbf{j}}a_{\mathbf{k}}b_{\mathbf{l}}}$ and
$V_{b_{\mathbf{j}}b_{\mathbf{k}}a_{\mathbf{l}}}$ denote the vector fields
corresponding to the systems in (5.100) from left to right, respectively. That
is, for example, $V_{a_{\mathbf{j}}b_{\mathbf{k}}b_{\mathbf{l}}}$ has
$a_{\mathbf{j}}$th coordinate given by
$\theta_{\mathbf{k}\mathbf{l}}b_{\mathbf{k}}b_{\mathbf{l}}$,
$b_{\mathbf{k}}$th coordinate given by
$\theta_{\mathbf{j}\mathbf{l}}a_{\mathbf{j}}b_{\mathbf{l}}$,
$b_{\mathbf{l}}$th coordinate given by
$-\theta_{\mathbf{j}\mathbf{k}}a_{\mathbf{j}}b_{\mathbf{k}}$ with all other
coordinates (in $\mathbf{R}^{\infty}$) equal to $0$.

Let $V_{\text{damp}}$ denote the vector field (on $\mathbf{R}^{\infty}$) such
that for every $\mathbf{j}\in\mathbf{Z}_{+}^{2}$, the $a_{\mathbf{j}}$th entry
is given by
$-\textbf{1}_{\mathscr{D}}(\mathbf{j})\lambda_{\mathbf{j}}a_{\mathbf{j}}$ and
$b_{\mathbf{j}}$th entry is given by
$-\mathbf{1}_{\mathscr{D}}(\mathbf{j})\lambda_{\mathbf{j}}b_{\mathbf{j}}$.
Fixing $\mathfrak{m}\in\mathbf{N}$, for any $\mathbf{j}\in\mathbf{Z}_{+}^{2}$
we let

$\displaystyle\beta_{a_{\mathbf{j}}}$
$\displaystyle=\sum_{\ell=1}^{\mathfrak{m}}\beta_{a_{\mathbf{j}}}^{\ell}\quad\text{
and
}\quad\beta_{b_{\mathbf{j}}}=\sum_{\ell=1}^{\mathfrak{m}}\beta_{b_{\mathbf{j}}}^{\ell}$

where
$\beta_{a_{\mathbf{j}}}^{\ell},\beta_{b_{\mathbf{j}}}^{\ell}\in\mathbf{R}$.
For any $\ell=1,2,\ldots,\mathfrak{m}$, we let $V_{\ell}$ denote the vector
field (on $\mathbf{R}^{\infty}$) such that for every
$\mathbf{j}\in\mathbf{Z}_{+}^{2}$, the $a_{\mathbf{j}}$the entry is
$\beta_{a_{\mathbf{j}}}^{\ell}$ and the $b_{\mathbf{j}}$th entry is
$\beta_{b_{\mathbf{j}}}^{\ell}$.

Consider now the collection of vector fields

(5.103)
$\displaystyle\mathscr{S}:=\\{V_{\text{damp}},V_{1},\ldots,V_{\mathfrak{m}}\\}\cup\bigcup_{\begin{subarray}{c}\mathbf{j}+\mathbf{k}-\mathbf{l}=0\\\
\mathbf{j},\mathbf{k},\mathbf{l}\in\mathbf{Z}_{+}^{2}\\\
\langle\mathbf{j},\mathbf{k}^{\perp}\rangle\neq
0\end{subarray}}\\{V_{a_{\mathbf{j}}a_{\mathbf{k}}a_{\mathbf{l}}},V_{a_{\mathbf{j}}b_{\mathbf{k}}b_{\mathbf{l}}},V_{b_{\mathbf{j}}a_{\mathbf{k}}b_{\mathbf{l}}},V_{b_{\mathbf{j}}b_{\mathbf{k}}a_{\mathbf{l}}}\\}$

Note that, if $\langle\mathbf{j},\mathbf{k}^{\perp}\rangle=0$, then the vector
fields
$V_{a_{\mathbf{j}}a_{\mathbf{k}}a_{\mathbf{l}}},V_{a_{\mathbf{j}},b_{\mathbf{k}}b_{\mathbf{l}}},V_{b_{\mathbf{j}},a_{\mathbf{k}}b_{\mathbf{l}}},V_{b_{\mathbf{j}}b_{\mathbf{k}}a_{\mathbf{l}}}$
with $\mathbf{j}+\mathbf{k}=\mathbf{l}$ are identically zero, hence the
restriction above in (5.103).

The focus in Section 5 is on a finite-dimensional version of the above set
$\mathscr{S}$ where we project the above vector fields onto indices belonging
to the finite lattice

(5.104)
$\displaystyle\mathscr{Z}^{2}:=\\{j=(j_{1},j_{2})\in\mathbf{Z}_{+}^{2}\,:\,|j_{1}|\leq
N,\,|j_{2}|\leq N\\}$

as well as project the vector fields onto the relevant finite-dimensional
Euclidean space. That is, we still consider the splitting in (5.103), but the
vector fields are now vector fields on
$\mathbf{R}^{2|\mathscr{Z}^{2}|}=\mathbf{R}^{2N(N+2)}=:\mathbf{R}^{d}$ rather
than $\mathbf{R}^{\infty}$ and all indices $\mathbf{j},\mathbf{k},\mathbf{l}$
belong to $\mathscr{Z}^{2}$. In Section 5, we make the slight abuse of
notation and give these vector fields the same names as above. That is, our
sought after splitting is described by the collection of vector fields on
$\mathbf{R}^{d}$ as in (5.8).

## References

  * [1] M. Abramowitz and I.A. Stegun.  Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.  Dover, 9th edition, 1964.
  * [2] A. Agazzi, J. C. Mattingly, and O. Melikechi.  Random splitting of fluid models: Positive Lyapunov exponents.  arXiv preprint arXiv:2210.02958, 2022.
  * [3] A. Agazzi, J.C. Mattingly, and O. Melikechi.  Random splitting of fluid models: Unique ergodicity and convergence.  Comm. Math. Phys., 401(1):497–549, 2023.
  * [4] Ludwig Arnold.  Random Dynamical Systems.  Springer Monographs in Mathematics. Springer-Verlag, Berlin, 1998.
  * [5] Y. Bakhtin and T. Hurth.  Invariant densities for dynamical systems with random switching.  Nonlinearity, 25(10):2937–2952, 2012.
  * [6] J. Bedrossian and K. Liss.  Stationary measures for stochastic differential equations with degenerate damping.  Probability Theory and Related Fields, pages 1–78, 2024.
  * [7] M. Benaïm, S. Le Borgne, F. Malrieu, and P.-A. Zitt.  Qualitative properties of certain piecewise deterministic Markov processes.  Ann. Inst. Henri Poincaré Probab. Stat., 51(3):1040–1075, 2015\.
  * [8] F. Bonetto, J. L. Lebowitz, and L. Rey-Bellet.  Fourier’s law: a challenge to theorists.  In Mathematical physics 2000, pages 128–150. Imp. Coll. Press, London, 2000.
  * [9] Evan Camrud.  Geometric Ergodicity of the Degenerate Stochastic Lorenz 96 Model in Four Dimensions.  PhD thesis, Iowa State University, 2022.
  * [10] N. Cuneo, P. Eckmann J, M. Hairer, and L. Rey-Bellet.  Non-equilibrium steady states for networks of oscillators.  Electron. J. of Probabi., 23:1 – 28, 2018.
  * [11] Randal D., Gersende F., E. Moulines, and P. Soulier.  Practical drift conditions for subgeometric rates of convergence.  Ann. of Appl. Probab., 14(3):1353 – 1377, 2004.
  * [12] Giuseppe Da Prato.  An Introduction to Infinite-Dimensional Analysis.  Springer Science & Business Media, 2006.
  * [13] M. H. A. Davis.  Piecewise-deterministic Markov processes: a general class of nondiffusion stochastic models.  J. Roy. Statist. Soc. Ser. B, 46(3):353–388, 1984.  With discussion.
  * [14] R. Douc, G. Fort, and A. Guillin.  Subgeometric rates of convergence of $f$-ergodic strong Markov processes.  Stoch. Proc. Appl., 119(3):897–923, 2009.
  * [15] J.-P. Eckmann, C.-A. Pillet, and L. Rey-Bellet.  Non-equilibrium statistical mechanics of anharmonic chains coupled to two heat baths at different temperatures.  Comm. Math. Phys., 201(3):657–697, 1999.
  * [16] M. Hairer.  Convergence of Markov processes.  Lecture notes, 18:26, 2010.
  * [17] M. Hairer and J. C. Mattingly.  Slow energy dissipation in anharmonic oscillator chains.  Comm. Pure Appl. Math., 62(8):999–1032, 2009.
  * [18] Rafail Khasminskii.  Stochastic Stability of Differential Equations.  Springer, 2012.
  * [19] Y. Kifer.  Ergodic theory of random transformations, volume 10 of Progress in Probability and Statistics.  Birkhäuser Boston, Inc., Boston, MA, 1986.
  * [20] O. Kiselev.  Uniform asymptotic behaviour of Jacobi-sn near a singular point. The lost formula from handbooks for elliptic functions.  arXiv preprint arXiv:1510.06602, 2015.
  * [21] Sean P. Meyn and Richard L. Tweedie.  Markov Chains and Stochastic Stability.  Springer Science & Business Media, 2012.
  * [22] P. Tuominen and R. L. Tweedie.  Subgeometric rates of convergence of f-ergodic Markov chains.  Adv. in Appl. Probab., 26(3):775–798, 1994.
  * [23] Brendan Williamson.  On SDEs with Partial Damping Inspired by the Navier-Stokes Equations.  PhD thesis, Duke University, 2019.