# The Fundamental Solution of an Elliptic Equation with Singular Drift

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: August 31, 2022)

###### Abstract.

For $n\geq 3$, we study the existence and asymptotic properties of the
fundamental solution for elliptic operators in nondivergence form,
${\mathcal{L}}(x,\partial_{x})=a_{ij}(x)\partial_{i}\partial_{j}+b_{k}(x)\partial_{k}$,
where the $a_{ij}$ have modulus of continuity $\omega(r)$ satisfying the
square-Dini condition and the $b_{k}$ are allowed mild singularities of order
$r^{-1}\omega(r)$. A singular integral is introduced that controls the
existence of the fundamental solution. We give examples that show the singular
drift $b_{k}\partial_{k}$ may act as a perturbation that does not dramatically
change the fundamental solution of
${\mathcal{L}}^{o}=a_{ij}\partial_{i}\partial_{j}$, or it could change an
operator ${\mathcal{L}}^{o}$ that does not have a fundamental solution to one
that does.

###### Key words and phrases:

Elliptic equations, nondivergence form, square-Dini condition, singular drift,
fundamental solution, asymptotics

###### 1991 Mathematics Subject Classification:

35A08 (Primary), 35B40, 35J15 (Secondary)

##  0\. Introduction

Consider a 2nd-order operator in nondivergence form with a 1st-order term

(1)
${\mathcal{L}}(x,\partial_{x})u=a_{ij}(x)\,\partial_{i}\partial_{j}u+b_{k}(x)\,\partial_{k}u,\quad\hbox{for
$x\in U$},$

where the coefficients $a_{ij}=a_{ji}$ and $b_{k}$ are measurable, real-valued
functions in an open set $U\subset{\mathbb{R}}^{n}$, $n\geq 3$; here and
throughout the paper we use the summation convention for repeated indices. We
assume ${\mathcal{L}}(x,\partial_{x})$ is pointwise elliptic, i.e. the matrix
$a_{ij}(x)$ is positive definite for each $x\in U$. Our objective is to study
the existence and asymptotics of its fundamental solution, i.e. a function
$F(x,y)$ satisfying $F(x,\cdot)\in L^{1}_{{\ell\mathrm{o}c}}(U)$ for each
$x\in U$ and

(2) $-{\mathcal{L}}(x,\partial_{x})F(x,y)=\delta(x-y)\quad\hbox{for}\ x,y\in
U,$

in a distributional sense that needs to be made clear. In the classical case
that the coefficients are sufficiently smooth (e.g. Hölder continuous), the
existence of a fundamental solution and its asymptotic properties are well-
known: cf. Miranda [17].

In [11] we studied this problem for
${\mathcal{L}}^{o}(x,\partial_{x})=a_{ij}(x)\,\partial_{i}\partial_{j}$ when
the modulus of continuity $\omega(r)$ for the $a_{ij}$ satisfies the square-
Dini condition

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

Notice that this assumption is weaker than $\lambda$-Hölder continuity since
(3) is satisfied by $\omega(r)=r^{\lambda}$ for $0<\lambda<1$; (3) is also
weaker than the Dini-condition, $\int_{0}^{1}r^{-1}\omega(r)\,dr<\infty$. The
hypotheses and conclusions of [11] are most easily stated if we fix $y=0$ and
let $U=B_{\varepsilon}=\\{x\in{\mathbb{R}}^{n}:|x|<\varepsilon\\}$ for
$\varepsilon$ sufficiently small. By an affine change of coordinates, we may
arrange $a_{ij}(0)=\delta_{ij}$, so we may assume

(4) $\sup_{|x|=r}|a_{ij}(x)-\delta_{ij}|\leq\omega(r)\quad\hbox{for}\
0<r<\varepsilon.$

It was found in [11] that the existence of a solution of

(5) $-{\mathcal{L}}^{o}(x,\partial_{x})F(x)=\delta(x)\quad\hbox{in}\
B_{\varepsilon},$

for $\varepsilon$ sufficiently small depends on the behavior of the integral

(6) $I^{o}(r)=\frac{1}{|S^{n-1}|}\int_{r<|z|<\varepsilon}\left(\hbox{tr}({\bf
A}_{z})-n\frac{\left\langle{\bf
A}_{z}z,z\right\rangle}{|z|^{2}}\right)\frac{dz}{|z|^{n}}$

as $r\to 0$; here ${\bf A}_{x}$ denotes the matrix $(a_{ij})$ and $\langle{\bf
A}_{z}z,z\rangle=a_{ij}(z)z_{i}z_{j}$. If $I^{o}(0)=\lim_{r\to 0}I^{o}(r)$
exists and is finite, then we showed there is a solution of (5). However, if
$I^{o}(r)\to-\infty$ as $r\to 0$, then we found a solution $F(x)$ of

(7) ${\mathcal{L}}^{o}(x,\partial_{x})F(x)=0\quad\hbox{in}\ B_{\varepsilon},$

which has a singularity at $x=0$, but $F(x)=o(|x|^{2-n})$ as $|x|\to 0$ and
$F(x)$ is not a solution of (5); this violates the “extended maximum
principle” of [5].

In the present paper, we generalize the results of [11] by including the 1st-
order term $b_{k}\partial_{k}$, which is called a drift term: cf. [2], [1]. If
the $b_{k}$ are bounded in $U$, they will not effect the existence of the
fundamental solution, so we will allow the drift to be singular at $x=0$, but
satisfy the condition:

(8) $\sup_{|x|=r}|x|\,|b_{k}(x)|\leq c\,\omega(r)\quad\hbox{for}\
0<r<\varepsilon.$

(Here and throughout this paper, $c$ denotes a generic constant.) We again
construct an unbounded solution $Z(x)$ of

(9) ${\mathcal{L}}(x,\partial_{x})Z(x)=0\quad\hbox{for}\ x\in
B_{\varepsilon}\backslash\\{0\\},$

and then check to see whether $Z$ has the appropriate singular behavior as
$|x|\to 0$ so that, in a distributional sense,

(10) $-{\mathcal{L}}(x,\partial_{x})Z(x)=C_{0}\,\delta(x)\quad\hbox{for some
constant $C_{0}$.}$

If so, then setting $F(x):=C_{0}^{-1}\,Z(x)$ defines a solution of

(11) $-{\mathcal{L}}(x,\partial_{x})F(x)=\delta(x)\quad\hbox{for}\ x\in
B_{\varepsilon}.$

To determine whether this can be done, the quantity (6) is generalized to

(12) $I(r)=\frac{1}{|S^{n-1}|}\int_{r<|z|<\varepsilon}\left(\hbox{tr}({\bf
A}_{z})-n\frac{\left\langle{\bf A}_{z}z,z\right\rangle}{|z|^{2}}+\langle{\bf
B}_{z},z\rangle\right)\frac{dz}{|z|^{n}},$

where ${\bf B}_{z}$ denotes the vector with components $b_{k}(z)$ and
$\langle{\bf B}_{z},z\rangle=b_{k}(z)\,z_{k}$. In general, $I(r)$ may have a
singularity at $r=0$, but it is weaker than logarithmic: for any $\lambda>0$
there exists $c_{\lambda}>0$ such that

(13) $|I(r)|\leq\lambda\,|\log r|+c_{\lambda}\quad\hbox{for}\
0<r<\varepsilon.$

Before stating the main results of this paper, we need some notation and an
additional assumption. For an open set $U$, $1<p<\infty$, and integer
$m=0,1,2$, let $W^{m,p}_{\ell\mathrm{o}c}(U)$ denote the Sobolev space of
functions whose derivatives up to order $m$ are locally $L^{p}$-integrable in
$U$. The $L^{p}$-mean of $w\in
L^{p}_{\ell\mathrm{o}c}({\mathbb{R}}^{n}\backslash\\{0\\})$ on the annulus
$A_{r}:=\\{x:r<|x|<2r\\}$ is defined by

(14a)
$M_{p}(w,r):=\left({\diagup\hskip-11.9501pt\int}_{A_{r}}|w(x)|^{p}\,dx\right)^{1/p}.$
Similarly, we define (14b)
$M_{1,p}(w,r)=rM_{p}(Dw,r)+M_{p}(w,r)\quad\hbox{for}\ w\in
W^{1,p}_{\ell\mathrm{o}c}({\mathbb{R}}^{n}\backslash\\{0\\}),$ (14c)
$M_{2,p}(w,r)=r^{2}M_{p}(D^{2}w,r)+M_{1,p}(w,r)\quad\hbox{for}\ w\in
W^{2,p}_{\ell\mathrm{o}c}({\mathbb{R}}^{n}\backslash\\{0\\}).$

For $p=\infty$ we can analogously define $M_{\infty}(w,r)$,
$M_{1,\infty}(w,r)$ and $M_{2,\infty}(w,r)$. The $L^{p}$-mean over annuli
centered at $y\not=0$ will be denoted by $M_{p}(w,r;y)$ and similarly for
$M_{1,p}(w,r;y)$ and $M_{2,p}(w,r;y)$. We also define

(15) $\sigma(r)=\int_{0}^{r}\omega^{2}(\rho)\frac{d\rho}{\rho},$

which satisfies $\sigma(r)\to 0$ as $r\to 0$ because of (3). Finally, as a
modulus of continuity we want $\omega(r)$ to be nondecreasing for $0<r<1$, but
we also assume for some $\kappa\in(0,1)$ that

(16) $\omega(r)\,r^{-1+\kappa}\quad\hbox{is nonincreasing on $0<r<1$.}$

This is natural since we are interested in moduli of continuity that vanish
slower than $r$ as $r\to 0$.

Our first result generalizes Theorem 1 in [11], which only applied to
${\mathcal{L}}^{o}(x,\partial_{x})$.

######  Theorem 1.

For $n\geq 3$ and $p\in(1,\infty)$, assume $a_{ij}=a_{ji}$ satisfy (4) while
the $b_{k}$ satisfy (8). Then for $\varepsilon>0$ sufficiently small, there is
a solution of (9) in the form

(17a)
$Z(x)=\int_{|x|}^{\varepsilon}s^{1-n}\,e^{I(s)}\,ds\,(1+\zeta(|x|))+v(x)\quad\hbox{for}\
|x|<\varepsilon,$ where $I(r)$ is given by (12), $\zeta(r)$ satisfies (17b)
$M_{2,p}(\zeta,r)\leq c\,\max(\omega(r),\sigma(r)),$ and $v(x)$ satisfies
(17c) $M_{2,p}(v,r)\leq c\,r^{2-n}\,e^{I(r)}\,\omega(r).$

Moreover, for any $u\in
W^{2,p}_{\ell\mathrm{o}c}(\overline{B_{\varepsilon}}\backslash\\{0\\})$ that
is a strong solution of ${\mathcal{L}}(x,\partial_{x})u=0$ in
$\overline{B_{\varepsilon}}\backslash\\{0\\}$ subject to the growth condition

(18) $M_{2,p}(u,r)\leq c\,r^{1-n+\varepsilon_{0}}\ \hbox{where}\
\varepsilon_{0}>0,$

there exist constants $c_{0}$, $c_{1}$ (depending on $u$) such that

(19) $u(x)=c_{0}\,Z(x)+c_{1}+w(x),$

where $M_{2,p}(w,r)\leq c\,r^{1-\varepsilon_{1}}$ for any $\varepsilon_{1}>0$.

We will prove this result in Section 1, but let us here observe that $I(r)$
satisfies $|I^{\prime}(r)|\leq c\,r^{-1}\,\omega(r)$, so we may integrate by
parts and take $p>n$ to conclude

(20) $Z(x)=\frac{|x|^{2-n}\,e^{I(|x|)}}{n-2}\left(1+\xi(x)\right),$

where $M_{1,\infty}(\xi,r)\leq c\,\max(\omega(r),\sigma(r))$. This shows that
$I(r)$ controls how closely $Z$ adheres to the fundamental solution of the
Laplacian.

Our second result generalizes Theorem 2 in [11] and shows that the existence
and finiteness of the limit $I(0)=\lim_{r\to 0}I(r)$ determines whether $Z$
solves (10) for some constant $C_{0}$.

######  Theorem 2.

Under the assumptions of Theorem 1 and $\varepsilon>0$ sufficiently small:

  * (i)

If $I(0)=\lim_{r\to 0}I(r)$ exists and is finite, then we can solve (10) in
$B_{\varepsilon}$ with $C_{0}=|S^{n-1}|\,e^{I(0)}$.

  * (ii)

If $I(r)\to-\infty$ as $r\to 0$, then solving (10) in $B_{\varepsilon}$ yields
$C_{0}=0$, and so $Z$ solves $-{\mathcal{L}}(x,\partial_{x})Z(x)=0$ in
$B_{\varepsilon}$, despite its singularity at $x=0$.

We will prove this result in Section 2, but let us observe that in case (i) we
have found a solution of (11) of the form

(21) $F(x)=\frac{|x|^{2-n}}{(n-2)|S^{n-1}|}\,(1+o(|x|))\quad\hbox{as}\ |x|\to
0,$

so the comparison with the fundamental solution of the Laplacian is even more
explicit. On the other hand, if $I(r)\to+\infty$ as $r\to 0$, then the
singular solution $Z(x)$ grows more rapidly as $|x|\to 0$ than the fundamental
solution for the Laplacian, and we are not able to solve (10).

Let us consider a simple example to illustrate the effect of the drift term on
the fundamental solution.

Example 1. Let $a_{ij}=\delta_{ij}$ so that our operator (1) becomes

(22) ${\mathcal{L}}(x,\partial_{x})=\Delta+b_{k}(x)\partial_{x},$

and the quantity $I(r)$ in (12) reduces to

(23) $I(r)=\frac{1}{|S^{n-1}|}\int_{r<|x|<\varepsilon}\frac{\langle{\bf
B}_{z},z\rangle}{|z|^{n}}\,dz.$

If $\omega(r)$ satisfies the Dini condition, then from condition (8) we easily
conclude that $I(0)$ exists and is finite; but this finite limit may exist
without the Dini condition: e.g. we could take
$b_{k}(x)=x_{k}\,\sin(|x|^{-1})\,\omega(|x|)/|x|^{2}$. In any case, provided
$I(0)$ is finite, we have a solution $F(x)$ of (10) that is comparable to the
fundamental solution of the Laplacian.

In [2], Cranston and Zhao consider operators of the form
${\mathcal{L}}=\frac{1}{2}\Delta+b\cdot\nabla$ with vector field $b$. Assuming
$U$ is a bounded Lipschitz domain and $b(x)$ satisfies the conditions

(24a) $\lim_{r\to 0}\,\sup_{x\in
U}\int_{|x-y|<r}\frac{|b(y)|^{2}}{|x-y|^{n-2}}dy=0=\lim_{r\to 0}\,\sup_{x\in
U}\int_{|x-y|<r}\frac{|b(y)|}{|x-y|^{n-1}}dy,$ they conclude that the Green’s
function $G(x,y)$ for ${\mathcal{L}}$ in $U$ exists and is comparable to the
Green’s function $G_{0}(x,y)$ for ${\mathcal{L}}_{0}=\frac{1}{2}\Delta$, i.e.
(24b) $c^{-1}\,G_{0}(x,y)\leq G(x,y)\leq c\,G_{0}(x,y)\quad\hbox{for}\ x,y\in
U,\ x\not=y.$

If $U$ contains the origin, then setting $F(x):=G(x,0)$ defines a solution of
(11), so let us compare our hypotheses and conclusions with those of [2]. If
$b$ has a singularity at $x=0$ of the form $|b(x)|=|x|^{-1}\,\omega(|x|),$ the
conditions in (24a) require

$\int_{0}^{r}\frac{\omega^{2}(\rho)}{\rho}\,d\rho<\infty\quad\hbox{and}\quad\int_{0}^{r}\frac{\omega(\rho)}{\rho}\,d\rho<\infty.$

The first of these is the square-Dini condition (3) that we have required,
while the second is the Dini condition that we have not required: to conclude
that the fundamental solution exists at $x=0$, we only require the function
$I(r)$ given in (23) to have a finite limit $I(0)=\lim_{r\to 0}I(r)$.
Moreover, (20) is a sharper estimate than (24b).

In Example 1, the drift term plays the role of a perturbation which, if not
too large, does not affect the existence of the fundamental solution. We now
consider an example where singular drift can convert an operator
${\mathcal{L}}^{o}$ that does not have a fundamental solution to one for which
a fundamental solution exists!

Example 2. Consider ${\mathcal{L}}^{o}=a_{ij}\partial_{i}\partial_{j}$ with
coefficients

(25) $a_{ij}(x)=\delta_{ij}+g(r)\frac{x_{i}x_{j}}{|x|^{2}},$

where $|g(r)|\leq\omega(r)$ with $\omega$ satisfying (3). Coefficients of the
form (25) were first considered by Gilbarg & Serrin [5], and have proven
useful in both nondivergence and divergence form equations (cf. [11],
[12],[13]). We can use (6) to calculate

(26) $I^{o}(r)=(1-n)\int_{r}^{1}g(\rho)\frac{d\rho}{\rho}.$

Note that $I^{o}(0)\\!=\\!\lim_{r\to 0}I^{o}(r)$ exists and is finite when
$\omega$ satisfies the Dini condition, but this finite limit may exist without
Dini continuity: e.g. $g(r)=\sin(r^{-1})\omega(r)$.

Now let us assume $g(r)>0$ and $\omega$ does not satisfy the Dini condition.
Then $I^{o}(r)\to-\infty$ as $r\to 0$, so the fundamental solution for
${\mathcal{L}}^{o}$ does not exist at $y=0$. However, if we add the first-
order coefficients

(27) $b_{k}(x)=(n-1)\,\frac{x_{k}}{|x|^{2}}\,[g(|x|)+g^{2}(|x|)]$

to obtain ${\mathcal{L}}$ as in (1), then we can use (12) to calculate

(28) $I(r)=(n-1)\int_{r}^{1}\frac{g^{2}(\rho)}{\rho}\,d\rho.$

Since $g^{2}(r)\leq\omega^{2}(r)$ and $\omega$ satisfies (3), we see that
$I(r)$ has a finite limit as $r\to 0$, and so we can solve (10) to conclude
the fundamental solution exists at $y=0$.

As previously stated, we will prove Theorem 1 in Section 1 and Theorem 2 in
Section 2; but in those sections we will also state and prove Corollaries 1
and 2 respectively, which show how the formulas for the singular solution $Z$
change when $y\not=0$ and $a_{ij}(y)\not=\delta_{ij}$. To formulate these
results, we need to generalize (4): for a given $y\in U$, choose $\varepsilon$
so that $0<\varepsilon<{\rm dist}(y,\partial U)$ and assume

(29) $\sup_{|x-y|=r}\|{\bf A}_{x}-{\bf A}_{y}\|\leq\omega(r)\quad\hbox{for}\
0<r<\varepsilon.$

We continue to assume the $b_{k}$ are only singular at $x=0$ and satisfy (8).
We also generalize the function $I(r)$ in (12) as

(30a)
$I_{y}(r)=\frac{1}{|S^{n-1}|}\int_{r<|z-y|<\varepsilon}H(z,y)\,\frac{dz}{|z-y|^{n}}\,,$
where the integrand $H(z,y)$ is (30b) ${\rm tr}({\bf A}_{z}{\bf
A}_{y}^{-1})-n\frac{\langle{\bf A}_{z}{\bf A}_{y}^{-1/2}(z-y),{\bf
A}_{y}^{-1/2}(z-y)\rangle}{|z-y|^{2}}+\langle{\bf B}_{z}{\bf
A}_{y}^{-1/2},(z-y)\rangle,$

with ${\bf A}_{y}^{-1}$ denoting the inverse matrix. If $y=0$ and ${\bf
A}_{0}=I$, then (30) coincides with (12); if $y\not=0$ and we also stipulate
$0<\varepsilon<|y|$, then ${\bf B}_{z}$ is bounded on $|z-y|<\varepsilon$, so
the last term in (30b) will not play a role in whether the limit $I_{y}(0)$
exists and is finite. Now we state the main result of this paper.

######  Theorem 3.

Suppose ${\mathcal{L}}(x,\partial_{x})$ as in (1) is an elliptic operator in a
bounded open set $U\subset{\mathbb{R}}^{n}$, $n\geq 3$, where the coefficients
$a_{ij}=a_{ji}$ are continuous functions with modulus of continuity
$\omega(r)$ satisfying (3). Suppose $U$ contains the origin and the $b_{k}$
satisfy (8) but otherwise are bounded in $U$. For each $y\in U$ assume that
the limit $I_{y}(0)=\lim_{r\to 0}I_{y}(r)$ exists and is finite. Then
${\mathcal{L}}(x,\partial_{x})$ has a fundamental solution $F(x,y)$ in $U$ and
it has the asymptotic behavior

(31) $F(x,y)=\frac{\langle{\bf
A}_{y}^{-1}(x-y),(x-y)\rangle^{\frac{2-n}{2}}}{(n-2)|S^{n-1}|\sqrt{\det{\bf
A}_{y}}}(1+o(1)))\quad\hbox{as}\ x\to y.$

This is proved in Section 3. Note that the leading asymptotic in (31) is
familiar as the “Levi function” that occurs in the classical case (cf. [17]).
As a consequence of Theorem 3, we see that the singular drift may affect the
existence of the fundamental solution but does not play a role in its
asymptotic behavior as $x\to y$.

Many of the arguments in this paper also appear in [11], but we have repeated
them here for the convenience of the reader. From a more general perspective,
the asymptotic analysis used here is related to that developed in [9].

Let us compare our results (in [11] and this paper) with recent work
estimating the singularity of the Green’s function for
${\mathcal{L}}^{o}(x,\partial_{x})=a_{ij}(x)\partial_{i}\partial_{j}$ when the
coefficients $a_{ij}$ satisfy the Dini mean oscillation (DMO) condition: cf.
[7], [8], [4]. If $U$ is a bounded $C^{1,1}$-domain and the $a_{ij}=a_{ji}$
are continuous functions on $\overline{U}$, then it is well-known that the
Green’s function $G(x,y)$ exists. If the $a_{ij}$ satisfy the DMO condition in
$U$, then [4] shows that, for any $x_{0}\in U$, the following limit holds:

(32) $\lim_{x\to x_{0}}|x-x_{0}|^{2-n}|G(x,x_{0})-G_{x_{0}}(x,x_{0})|=0,$

where $G_{x_{0}}$ denotes that Green’s function for the constant coefficient
operator ${\mathcal{L}}^{o}(x_{0},\partial_{x})$ in $U$. Since the Green’s
function is a particular fundamental solution and since (32) is equivalent to
(31), it is natural to compare the hypotheses of the two results. As shown in
[3], there are coefficients that are DMO but do not satisfy our square-Dini
condition (3). On the other hand, there are coefficients of the form (25)
which satisfy (3) but are not DMO. In fact, as shown in [13], with

(33) $g(r)=\sin(|\log r|)\,|\log r|^{-\gamma},$

the coefficients in (25) are DMO only for $\gamma>1$, but they satisfy (3) for
$\gamma>1/2$ and the limit $I(0)$ exists and is finite for all $\gamma>0$.
This example with $1/2<\delta\leq 1$ shows that the results of [4] do not
cover our results for ${\mathcal{L}}^{o}(x,\partial_{x})$, let alone the
operator (1) with singular drift.

##  1\. Construction of the singular solution $Z$

In this section, we will not only prove Theorem 1, but we will state and prove
Corollary 1, which shows how the formulas change when we no longer assume
$y=0$ and $a_{ij}(y)=\delta_{ij}$. Instead of constructing $Z(x)$ in a small
ball, we replace the condition that $\omega(r)$ satisfies (3) with

(34) $\sigma(1)=\int_{0}^{1}\omega^{2}(\rho)\frac{d\rho}{\rho}<\mu^{2}$

where $\mu>0$ is sufficiently small, and then show existence in the unit ball
$B_{1}$. In fact, with $\kappa\in(0,1)$ as in (16), this also implies

(35) $\omega(r)<c_{\kappa}\,\mu\quad\hbox{for}\ 0<r\leq 1,$

since

$\mu^{2}>\int_{0}^{r}\omega^{2}(\rho)\frac{d\rho}{\rho}\geq\omega^{2}(r)r^{-2+2\kappa}\int_{0}^{r}\rho^{1-2\kappa}d\rho=\frac{\omega^{2}(r)}{2(1-\kappa)}.$

Moreover, it will be useful to consider ${\mathcal{L}}$ on all of
${\mathbb{R}}^{n}$, so we assume

(36) $a_{ij}(x)=\delta_{ij}\ \hbox{and}\ b_{k}(x)=0\ \hbox{for}\ |x|>1,$

and construct a solution of ${\mathcal{L}}Z=0$ on
${\mathbb{R}}^{n}\backslash\\{0\\}$.

Proof of Theorem 1. As in [11], we use spherical means: for a function $f(x)$
we denote its mean value over the sphere $|x|=r$ by $\overline{f}(r)$:

$\overline{f}(r)={\diagup\hskip-11.9501pt\int}_{S^{n-1}}f(r\theta)\,d\theta,$

where $S^{n-1}$ is the unit sphere, the slashed integral denotes mean value,
$r=|x|$, $\theta=x/|x|\in S^{n-1}$, and $d\theta$ denotes standard surface
measure on $S^{n-1}$. Let us write

(37a) $Z(x)=h(|x|)+v(x),\quad\hbox{where}\ h(r):=\overline{Z}(r),$ so that
(37b) $\overline{v}(r)=0.$

If we take the spherical mean of the equation ${\mathcal{L}}(h+v)=0$, we
obtain

(38a)
$\alpha(r)\,h^{\prime\prime}+\left[\frac{\alpha_{n}(r)-\alpha(r)+\beta(r)}{r}\right]\,h^{\prime}+\overline{{\mathcal{L}}\,v}(r)=0,$
where (38b)
$\alpha(r):={\diagup\hskip-11.9501pt\int}_{S^{n-1}}a_{ij}(r\theta)\theta_{i}\theta_{j}\,d\theta,\qquad\alpha_{n}(r):={\diagup\hskip-11.9501pt\int}_{S^{n-1}}a_{ii}(r\theta)\,d\theta,\\\
$ and (38c)
$\beta(r):=r\,{\diagup\hskip-11.9501pt\int}_{S^{n-1}}b_{k}(r\theta)\theta_{k}\,d\theta.$

In terms of these, instead of (12) we can write

(39) $I(r)=\int_{r}^{1}(\alpha_{n}(s)-n\,\alpha(s)+\beta(s))\frac{ds}{s}.$

Using (4), we see that

(40a)
$\alpha(r)=1+O(\omega(r))\quad\hbox{and}\quad\alpha_{n}(r)=n+O(\omega(r))\quad\hbox{as}\
r\to 0,\\\ $ and using (8) we see that (40b) $|\beta(r)|\leq\omega(r).$

Hence the integrand in (39) is bounded by $c\,\omega(s)\,s^{-1}$. Since
$\omega(r)$ need not satisfy the Dini condition, we do not know whether $I(r)$
has a finite limit as $r\to 0$; this is, of course, the significance of
Theorem 2.

Turning to the $\overline{{\mathcal{L}}\,v}(r)$ term, if we write
$a_{ij}\partial_{i}\partial_{j}v=\widetilde{a}_{ij}\partial_{i}\partial_{j}v+\Delta
v$ where $\widetilde{a}_{ij}:=(a_{ij}-\delta_{ij})$ and use $\overline{\Delta
v}=\Delta\overline{v}=0$, we see that

$\overline{{\mathcal{L}}v}(r)=\overline{\widetilde{a}_{ij}\partial_{i}\partial_{j}v}(r)+\overline{b_{k}\partial_{k}v}(r).$

Notice that $|\overline{{\mathcal{L}}v}(r)|\leq
c\,\omega(r)\,(|\overline{D^{2}v}|+r^{-1}|\overline{Dv}|)$ for $0<r<1$ and
$\overline{{\mathcal{L}}v}(r)=0$ for $r>1$. Since $\alpha(r)\to 1$ as $r\to 0$
and $\alpha(r)=1$ for $r>1$, we may assume $\alpha(r)\geq\delta>0$ for
$0<r<\infty$. Hence we may divide (38a) by $\alpha(r)$ and replace
$h^{\prime}$ by $g$ to obtain

(41a) $g^{\prime}+\left[\frac{n-1+R(r))}{r}\right]\,g=B[v](r),$ where (41b)
$R(r)=\frac{\alpha_{n}(r)+\beta(r)}{\alpha(r)}-n$ satisfies $|R(r)|\leq
c\,\omega(r)$ as $r\to 0$ and $R(r)=0$ for $r>1$, and $B[v](r)$ satisfies
(41c) $|B[v](r)|\leq
c\,\omega(r)\left(|\overline{D^{2}v}(r)|+r^{-1}|\overline{Dv}(r)|\right)\quad\hbox{for}\
0<r<1$ and $B[v](r)=0$ for $r>1$. Moreover, the monotonicity of $\omega(r)$
and (16) imply (41d) $\sup_{r<\rho<2r}\omega(\rho)\leq c\,\omega(r),$ so we
consequently obtain (41e) $r^{2}M_{p}(B[v],r)\leq c\,\omega(r)\,M_{2,p}(v,r).$

Solving (41a) involves the integrating factor $r^{n-1}E_{-}(r)$, where we
introduce

(42)
$E_{\pm}(r)=\exp\left[\pm\int_{r}^{\infty}R(s)\,\frac{ds}{s}\right]=\exp\left[\pm\int_{r}^{1}R(s)\,\frac{ds}{s}\right]=\frac{1}{E_{\mp}(r)}.$
Notice that $E_{-}(r)E_{+}(\rho)=\exp(\int_{\rho}^{r}R(s)\,s^{-1}\,ds)$ so by
(35) we have (43a)
$\left(\frac{\rho}{r}\right)^{c_{\kappa}\,\mu}\leq\exp\left(\pm\int_{\rho}^{r}R(s)\frac{ds}{s}\right)\leq\left(\frac{r}{\rho}\right)^{c_{\kappa}\,\mu}\quad\hbox{for}\
0<\rho<r\leq 1.$ In particular, we have (43b) $c_{1}\,E_{\pm}(r)\leq
E_{\pm}(\rho)\leq c_{2}\,E_{\pm}(r)\quad\hbox{for}\ r<\rho<2r,$ and for any
$f\in L^{p}_{{\ell\mathrm{o}c}}({\mathbb{R}}^{n}\backslash\\{0\\})$ and
$\nu\in{\mathbb{R}}$ we have (43c) $M_{p}(|x|^{\nu}E_{\pm}(|x|)\,f(x),r)\leq
c\,r^{\nu}E_{\pm}(r)\,M_{p}(f,r).$

While $E_{+}(r)$ is used to solve (41a), we observe that it is equivalent to
$e^{I(r)}$. In fact,

(44a) $E_{+}(r)=A\,e^{I(r)}(1+\tau(r))$ where $|R(s)(1-\alpha(s))|\leq
c\,\omega^{2}(s)$ implies (44b)
$A=\exp\left[\int_{0}^{1}R(s)(1-\alpha(s))s^{-1}\,ds\right]\quad\hbox{is
finite and positive,}$ and (44c)
$\tau(r)=\exp\left[-\int_{0}^{r}R(s)[1-\alpha(s)]\,\frac{ds}{s}\right]-1\quad\hbox{satisfies}\
|\tau(r)|\leq c\,\sigma(r).$

Hence, for some constants $c_{1},c_{2}$ we have

(45) $c_{1}\,E_{+}(r)\leq e^{I(r)}\leq c_{2}\,E_{+}(r).$

We consider (41a) as an ODE that depends on $v\in Y$, where $Y$ is the Banach
space of functions $v\in
W^{2,p}_{{\ell\mathrm{o}c}}({\mathbb{R}}^{n}\backslash\\{0\\})$ for which

(46)
$\|v\|_{Y}:=\sup_{0<r<1}\frac{M_{2,p}(v,r)\,r^{n-2}}{\omega(r)\,e^{I(r)}}+\sup_{r>1}\frac{M_{2,p}(v,r)\,r^{n-1}}{\mu}<\infty.$

If we let $\phi(r)=r^{n-1}E_{-}(r)\,g(r)$, then (41a) implies
$\phi^{\prime}(r)=r^{n-1}E_{-}(r)\,B[v](r).$ This can be integrated to find

(47) $\phi(r)=\phi(0)+\int_{0}^{r}\rho^{n-1}E_{-}(\rho)\,B[v](\rho)\,d\rho,$

where $\phi(0)$ is an arbitrary constant. Of course, for (47) to be valid we
need to know that $r^{n-1}E_{-}(r)\,B[v](r)$ is integrable at $r=0$. In fact,
we will show below that for $v\in Y$ we have

(48) $\int_{0}^{r}\rho^{n-1}\,E_{-}(\rho)\,|B[v](\rho)|\,d\rho\leq
c\,\mu^{2}\quad\hbox{for all}\ r>0,$

where $c$ may be taken uniformly for all $v\in Y$ with $\|v|_{Y}\leq 1$. Hence
(47) is valid and we conclude

(49a)
$g(r)=h^{\prime}(r)=r^{1-n}\,E_{+}(r)\left[\phi(0)+\int_{0}^{r}\rho^{n-1}\,E_{-}(\rho)\,B[v](\rho)\,d\rho\right]$
and (49b)
$h^{\prime\prime}(r)=\frac{1-n-R(r)}{r^{n}}\,E_{+}(r)\left[\phi(0)+\int_{0}^{r}\rho^{n-1}E_{-}(\rho)B[v](\rho)\,d\rho\right]+B[v](r).$

To verify (48), we observe that $v\in Y$ implies $M_{2,p}(v,r)\leq
c\,\omega(r)\,r^{2-n}\,E_{+}(r)$, so we can use (41c), (41d), (43b), and
Hölder’s inequality to conclude

$\displaystyle\int_{r}^{2r}\rho^{n-1}\,E_{-}(\rho)\,|B[v](\rho)|\,d\rho$
$\displaystyle\leq
c\,E_{-}(r)\,\omega(r)\int_{r}^{2r}\rho^{n-1}\left(|\overline{D^{2}v}(\rho)|+\rho^{-1}|\overline{Dv}(\rho)|\right)d\rho$
$\displaystyle\leq
c\,E_{-}(r)\omega(r)\int_{r<|x|<2r}\left(|D^{2}v(x)|+|x|^{-1}|Dv(x)|\right)dx$
$\displaystyle\leq c\,E_{-}(r)\,\omega(r)\,r^{n-2}\,M_{2,p}(v,r)\leq
c\,\omega^{2}(r).$

Now if we write

$\int_{0}^{r}\rho^{n-1}\,E_{-}(\rho)\,|B[v](\rho)|\,d\rho=\sum_{j=0}^{\infty}\int_{r/2^{j+1}}^{r/2^{j}}\rho^{n-1}\,E_{-}(\rho)\,|B[v](\rho)|\,d\rho,$

then we obtain

(50) $\displaystyle\int_{0}^{r}\rho^{n-1}\,E_{-}(\rho)\,|B[v](\rho)|\,d\rho$
$\displaystyle\leq
c\sum_{j=0}^{\infty}\omega^{2}\left(\frac{r}{2^{j+1}}\right)$
$\displaystyle\leq
c\int_{0}^{r}\omega^{2}(\rho)\frac{d\rho}{\rho}=c\,\sigma(r)<c\,\mu^{2}.$

This confirms (48) with $c$ uniform for $\|v\|_{Y}\leq 1$.

We also have a PDE for $v$ that depends upon $h$. From
${\mathcal{L}}Z-\overline{{\mathcal{L}}Z}=0$ we find:

(51) $\displaystyle-\Delta v=$
$\displaystyle\widetilde{a}_{ij}\partial_{i}\partial_{j}h-\overline{\widetilde{a}_{ij}\partial_{i}\partial_{j}h}+\widetilde{a}_{ij}\partial_{i}\partial_{j}v-\overline{\widetilde{a}_{ij}\partial_{i}\partial_{j}v}$
$\displaystyle+b_{k}\partial_{k}h-\overline{b_{k}\partial_{k}h}++b_{k}\partial_{k}v-\overline{b_{k}\partial_{k}v}.$

For a given $v\in Y$, we solve (41a) for $h$ and use (49a) and (49b) to write

$b_{k}\partial_{k}h=r^{-n}E_{+}(r)\left[\phi(0)+\int_{0}^{r}\rho^{n-1}E_{-}(\rho)B[v](\rho)d\rho\right]\psi_{1}$

and

$\widetilde{a}_{ij}\partial_{i}\partial_{j}h=r^{-n}E_{+}(r)\left[\phi(0)+\int_{0}^{r}\rho^{n-1}E_{-}(\rho)B[v](\rho)\,d\rho\right]\,\psi_{2}+B[v]\,\widetilde{a}_{ij}\theta_{i}\theta_{j},$

where

(52)
$\psi_{1}(r\theta)=r\,b_{k}(r\theta)\theta_{k}\quad\hbox{and}\quad\psi_{2}(r\theta)=\widetilde{a}_{ii}(r\theta)-(n+R(r))\widetilde{a}_{ij}(r\theta)\theta_{i}\theta_{j}$

also satisfiy $|\psi_{i}(r\theta)|\leq c\,\omega(r)$ for $i=1,2$. Plugging
this into (51), we obtain an equation of the form $-\Delta v=F[v]$. We want to
apply $K$, convolution by the fundamental solution of the Laplacian, to solve
this, but there could be a problem: $F[v]$ may not be integrable at $x=0$.
However, from (51) we see that $\overline{F[v]}=0$ and multiplying $F[v](x)$
by $|x|$ makes it integrable, so we can use Proposition 1 in Appendix A.

Applying $K$ to both sides of (51), we obtain an equation for $v$ alone:

(53) $v+S_{1}v+S_{2}v=\phi(0)w,$

where

$w(x)=K_{y\to
x}\left(|y|^{-n}E_{+}(|y|)\,\left[\psi(y)-\overline{\psi}(|y|)\right]\right)$
$S_{1}v\\!=\\!-K_{y\to
x}\left(\\!|y|^{-n}E_{+}(|y|)\int_{0}^{|y|}\rho^{n-1}E_{-}(\rho)B[v](\rho)d\rho\left[\psi(y)-\overline{\psi}(|y|)\right]\\!\right)$
$S_{2}v=-K_{y\to
x}\left(B[v](|y|)\left[\widetilde{a}_{ij}\theta_{i}\theta_{j}-\overline{\widetilde{a}_{ij}\theta_{i}\theta_{j}}\right]+\widetilde{a}_{ij}\theta_{i}\theta_{j}-\overline{\widetilde{a}_{ij}\theta_{i}\theta_{j}}\right),$

and we have let $\psi:=\psi_{1}+\psi_{2}$. To find $v\in Y$ satisfying (53) we
need to show $w\in Y$ and $S_{i}:Y\to Y$ has small operator norm for $i=1,2$.

To show $w\in Y$, we must estimate $M_{2,p}(w,r)$ for $0<r<1$ and for $r>1$.
We apply Proposition 1 in Appendix A to
$f(x)=|x|^{-n}E_{+}(|x|)(\psi(x)-\overline{\psi}(|x|)$, which vanishes for
$|x|>1$, to conclude

(54)
$M_{2,p}(w,r)\\!\leq\\!c\left(\\!r^{1-n}\\!\int_{0}^{r}E_{+}(\rho)\omega(\rho)\,d\rho\\!+\\!r\\!\int_{r}^{1}E_{+}(\rho)\omega(\rho)\rho^{-n}\,d\rho\\!\right)\
\hbox{for}\ 0<r<1.$

We can use (43a) and the monotonicity of $\omega(r)$ to estimate

(55a) $\int_{0}^{r}E_{+}(\rho)\,\omega(\rho)\,d\rho\leq
E_{+}(r)\,\omega(r)\,r^{c\mu}\int_{0}^{r}\rho^{-c\mu}d\rho=c\,r\,E_{+}(r)\,\omega(r)$
and similarly (55b)
$\int_{r}^{1}E_{+}(\rho)\,\omega(\rho)\,\rho^{-n}\,d\rho\leq
c\,r^{1-n}\,E_{+}(r)\,\omega(r).$

Using these in (54), we obtain $M_{2,p}(w,r)\,r^{n-2}\leq
c\,\omega(r)\,E_{+}(r)\quad\hbox{for}\ 0<r<1,$ and we can invoke (45) to
replace $E_{+}(r)$ by $e^{I(r)}$ as required in the norm for $Y$. Meanwhile,
for $r>1$ we use (35) and $E_{+}(\rho)\leq\rho^{-c\,\mu}$ for $0<\rho<1$ to
conclude

$M_{2,p}(w,r)\leq c\,r^{1-n}\int_{0}^{1}E_{+}(\rho)\,\omega(\rho)\,d\rho\leq
c\,\mu\,r^{1-n}.$

Hence $M_{2,p}(w,r)\,r^{n-1}\leq c\,\mu\quad\hbox{for}\ r>1.$ These estimates
confirm that $w\in Y$.

Next let us show that $S_{1}:Y\to Y$ with small operator norm. We assume
$\|v\|_{Y}\leq 1$ and we want to estimate $M_{2,p}(S_{1}v,r)$ separately for
$0<r<1$ and $r>1$. Using (48) we see that the function

$f_{1}(y)=|y|^{-n}E_{+}(|y|)\int_{0}^{|y|}\rho^{n-1}E_{-}(\rho)\,B[v](\rho)\,d\rho\left(\psi(y)-\overline{\psi}(|y|)\right)$

satisfies $M_{p}(f_{1},r)\leq c\,\mu^{2}\,E_{+}(r)\,\omega(r)\,r^{-n}$ for
$0<r<1$ and $M_{p}(f_{1},r)=0$ for $r>1$. Since $S_{1}v=-Kf_{1}$, we can apply
Proposition 1 in Appendix 1 to obtain

$M_{2,p}(S_{1}v,r)\leq
c\,\mu^{2}\left(r^{1-n}\int_{0}^{r}E_{+}(\rho)\omega(\rho)\,d\rho+r\int_{r}^{1}E_{+}(\rho)\omega(\rho)\rho^{-n}\,d\rho\right).$

Using (55a) and (55b), we conclude $M_{2,p}(S_{1}v,r)\,r^{n-2}\leq
c\,\mu^{2}\,\omega(r)\,E_{+}(r)$ for $0<r<1$. On the other hand, for $r>1$,
Proposition 1 in Appendix 1 implies

$M_{2,p}(S_{1}v,r)\leq
c\,\mu^{2}\,r^{1-n}\int_{0}^{1}E_{+}(\rho)\,\omega(\rho)\,d\rho\leq
c\,\mu^{3}\,r^{1-n},$

so $M_{2,p}(S_{1}v,r)\,r^{n-1}\leq c\,\mu^{3}$ for $r>1$. Combining these
estimates, we see that $S_{1}:Y\to Y$ has small operator norm.

Finally, we show that $S_{2}:Y\to Y$ with small operator norm. Again we assume
$\|v\|_{Y}\leq 1$ and estimate $M_{2,p}(S_{2}v,r)$ separately for $0<r<1$ and
$r>1$. Notice that the function

$f_{2}=B[v]\,(\widetilde{a}_{ij}\theta_{i}\theta_{j}-\overline{\widetilde{a}_{ij}\theta_{i}\theta_{j}})$

satisfies

$M_{p}(f_{2},r)\leq\omega(r)M_{p}(B[v],r)\leq
c\,\omega^{3}(r)\,E_{+}(r)\,r^{-n}\quad\hbox{for}\ 0<r<1,$

where $c$ is independent of $v$, and $M_{p}(f_{2},r)=0$ for $r>1$. Similarly,
the function

$f_{3}=\widetilde{a}_{ij}\partial_{i}\partial_{j}v-\overline{\widetilde{a}_{ij}\partial_{i}\partial_{j}v}$

satisfies

$M_{p}(f_{3},r)\leq\omega(r)\,M_{p}(D^{2}v,r)\leq\omega^{2}(r)\,E_{+}(r)r^{-n}\quad\hbox{for}\
0<r<1,$

and $M_{p}(f_{3},r)=0$ for $r>1$. For $0<r<1$, we apply Proposition 1 in
Appendix 1 to $S_{2}v=-K(f_{2}+f_{3})$ to conclude

$M_{2,p}(S_{2}v,r)\leq
c\left(r^{1-n}\int_{0}^{r}\omega^{2}(\rho)\,E_{+}(\rho)\,d\rho+r\int_{r}^{1}\omega^{2}(\rho)\,E_{+}(\rho)\,\rho^{-n}\,d\rho\right).$

Using (35), (55a), and (55b), we conclude that $M_{2,p}(S_{2}v,r)\,r^{n-2}\leq
c\,\mu\,\omega(r)\,E_{+}(r)$ for $0<r<1$. Meanwhile, for $r>1$, we use (35)
and (43a) to estimate

$M_{2,p}(S_{2}v,r)\leq
c\,r^{1-n}\int_{0}^{1}\omega^{2}(\rho)E_{+}(\rho)\,d\rho\leq
c\,\mu^{2}\,r^{1-n}\int_{0}^{1}\rho^{-c\,\mu}\,d\rho\leq c\,\mu^{2}\,r^{1-n}.$

Consequently, $M_{2,p}(S_{2}v,r)\,r^{n-1}/\mu\leq c\,\mu$. These estimates
show that $S_{2}:Y\to Y$ has small operator norm.

Since $S_{1}+S_{2}$ has small operator norm on $Y$, we conclude that (53) has
a unique solution $v\in Y$, depending on the choice of the constant
$c_{*}=\phi(0)$. But once $c_{*}$ and $v$ are known, we find $g(r)$ from (49a)
and integrate to find $h(r)$:

(56)
$h(r)=\int_{r}^{\infty}s^{1-n}E_{+}(s)\left[c_{*}+\int_{0}^{s}\rho^{n-1}E_{-}(\rho)B[v](\rho)\,d\rho\right]ds+c_{2}$

where $c_{2}$ is an arbitrary constant. To show that the solution
$Z(x)=h(|x|)+v(x)$ is of the form (17), we choose $c_{*}$ to enable us to
replace $E_{+}(s)$ by $e^{I(s)}$; recalling (44a) we see that we should take
$c_{*}=A^{-1}$ and write $h(r)=h_{0}(r)+h_{1}(r)+c$ where

(57a) $h_{0}(r)=\int_{r}^{1}s^{1-n}\,e^{I(s)}\,ds$ and (recalling (44c)) (57b)
$h_{1}(r)=\int_{r}^{1}s^{1-n}e^{I(s)}\tau(s)\,ds+\int_{r}^{1}s^{1-n}E_{+}(s)\int_{0}^{s}\rho^{n-1}E_{-}(\rho)B[v](\rho)\,d\rho\,ds.$

Integrating by parts and using $|I^{\prime}(s)|\leq c\,s^{-1}\omega(s)$ we can
show

$\left|h_{0}(r)-\frac{r^{2-n}}{n-2}e^{I(r)}\right|\leq
c\,r^{2-n}\,e^{I(r)}\,\omega(r).$

Using $|\tau(s)|\leq c\,\sigma(s)$ and (50), we can also estimate

$|h_{1}(r)|\leq c\,r^{2-n}\,e^{I(r)}\max(\omega(r),\sigma(r))\quad\hbox{for}\
0<r<1.$

If we define

(58) $\zeta(r)=\frac{h_{1}(r)}{h_{0}(r)}\quad\hbox{for}\ 0<r<1,$

then we can easily estimate $|\zeta(r)|,|r\zeta^{\prime}(r)|\leq
c\,\max(\omega(r),\sigma(r))$. To estimate $\zeta^{\prime\prime}$, let us
write
$h_{0}\zeta^{\prime\prime}=h_{1}^{\prime\prime}-h_{0}^{\prime\prime}\zeta-2h_{0}^{\prime}\zeta^{\prime}$
where

(59a)
$h_{0}^{\prime\prime}(r)=(n-1)r^{-n}e^{I(r)}-r^{1-n}e^{I(r)}I^{\prime}(r)$ and
(59b) $\displaystyle h_{1}^{\prime\prime}(r)$
$\displaystyle=r^{-n}e^{I(r)}[(n-1)\tau(r)-rI^{\prime}(r)\tau(r)-r\tau^{\prime}(r)]$
$\displaystyle+r^{-n}E_{+}(r)(n-1+R(r))\int_{0}^{r}\rho^{n-1}\,R(\rho)\,B[v](\rho)\,d\rho-B[v](r).$

We can estimate $h_{0}^{\prime\prime}$ and $h_{0}^{\prime}$ pointwise, but
$h_{1}^{\prime\prime}$ needs to be estimated in $M_{p}$. However, using (41e)
and $v\in Y$, we may conclude $M_{p}(r^{2}\zeta^{\prime\prime},r)\leq
c\,\max(\omega(r),\sigma(r))$. Combining with the estimates of the lower-order
derivatives, we have shown (17b). Since (17c) follows from $v\in Y$, we have
proved the first part of Theorem 1.

Let us now turn to the second part of Theorem 1. Suppose $u\in
W^{2,p}_{{\ell\mathrm{o}c}}(\overline{B_{1}}\backslash\\{0\\})$ satisfies
${\mathcal{L}}u=0$ in $\overline{B_{1}}\backslash\\{0\\}$ and
$M_{2,p}(u,r)\leq c\,r^{1-n+\varepsilon_{0}}$ for some $\varepsilon_{0}>0$; we
want to show $u$ is of the form (19). We will use properties of the bounded
linear map

(60) ${\mathcal{L}}:W^{2,p}_{\delta_{0},\delta_{1}}({\mathbb{R}}^{n}_{o})\to
L^{p}_{\delta_{0}+2,\delta_{1}+2}({\mathbb{R}}^{n}_{o}),$

where ${\mathbb{R}}^{n}_{o}={\mathbb{R}}^{n}\backslash\\{0\\}$ and
$W^{2,p}_{\delta_{0},\delta_{1}}({\mathbb{R}}^{n}_{o}),L^{p}_{\delta_{0}+2,\delta_{1}+2}({\mathbb{R}}^{n}_{o})$
are the weighted Sobolev spaces that are defined in Appendix B. Since we are
interested in the behavior of functions at the origin, we fix
$\delta_{1}\in(-n/p,-2+n/p^{\prime})$ and allow $\delta_{0}$ to vary. Note
that (60) is a perturbation of (82), and the norm of the difference
${\mathcal{L}}-\Delta$ depends on the magnitude of

$\sup_{|x|<1}\left(|a_{ij}(x)-\delta_{ij}|+|x|\,|b_{k}(x)|\right).$

Thus, provided we take $\mu$ in (34) sufficiently small, we can arrange that
(60) and (82) are not only Fredholm for exactly the same values of
$\delta_{0}$ and $\delta_{1}$, but the nullity and deficiency of the two maps
agree. Hence, with $\delta_{1}$ satisfying (83), we find

  * (i)

(60) is an isomorphism for $\delta_{0}\in(-n/p,-2+n/p^{\prime})$;

  * (ii)

(60) is surjective with nullity $1$ for
$\delta_{0}\in(-2+n/p^{\prime},-1+n/p^{\prime})$;

  * (iii)

(60) is injective with deficiency $1$ for $\delta_{0}\in(-n/p-1,-n/p)$.

Introduce a cutoff function $\chi\in C_{0}^{\infty}(B_{1})$ with $\chi=1$ on
$B_{1/2}$, then $M_{2,p}(\chi u,r)\leq c\,r^{1-n+\varepsilon_{0}}$ implies
$\chi u\in W^{2,p}_{\delta_{0}}(B_{1})$ provided
$\delta_{0}>-1-\varepsilon_{0}+n/p^{\prime}$. Let us choose
$\delta^{+}_{0}\in(-1-\varepsilon_{0}+n/p^{\prime},-1+n/p^{\prime})$ and let
$f={\mathcal{L}}(\chi u)$. Since $f=0$ for $|x|<1/2$ and for $|x|>1$, we have
$f\in L^{p}_{\delta_{0}+2,\delta_{1}+2}({\mathbb{R}}^{n}_{o})$ for all
$\delta_{0}$ so let us choose $\delta_{0}\in(-n/p,-2+n/p^{\prime})$. By (i) we
can find $v={\mathcal{L}}^{-1}f\in
W_{\delta_{0},\delta_{1}}^{2,p}({\mathbb{R}}^{n}_{o})$. Then $\chi u-v\in
W^{2,p}_{\delta_{0}^{+},\delta_{1}}({\mathbb{R}}^{n}_{o})$ satisfies
${\mathcal{L}}(\chi u-v)=0$. By (ii),
${\mathcal{L}}:W^{2,p}_{\delta_{0}^{+},\delta_{1}}({\mathbb{R}}^{n}_{o})\to
L^{p}_{\delta_{0}^{+}+2,\delta_{1}+2}({\mathbb{R}}^{n}_{o})$ has nullity 1.
Since ${\mathcal{L}}Z=0$, the nullspace must be spanned by $Z$ and so $\chi
u-v=c_{0}\,Z$ for some constant $c_{0}$.

It only remains to show that $v=c_{1}+w$ where $M_{2,p}(w,r)\leq
c\,r^{1-\varepsilon_{1}}$ for any $\varepsilon_{1}>0$. Let us pick
$\delta_{0}^{-}\in(-1-n/p,-n/p)$ so that by (iii), the map

(61)
${\mathcal{L}}:W^{2,p}_{\delta^{-}_{0},\delta_{1}}({\mathbb{R}}^{n}_{o})\to
L^{p}_{\delta^{-}_{0}+2,\delta_{1}+2}({\mathbb{R}}^{n}_{o}),$

is injective with deficiency 1. Let $\zeta$ be a linear functional on
$L^{p}_{\delta^{-}_{0}+2,\delta_{1}+2}({\mathbb{R}}^{n}_{o})$ that vanishes on
the image of (61). Note that ${\mathcal{L}}\chi=0$ for $|x|<1/2$ and for
$|x|>1$, so ${\mathcal{L}}\chi\in
L^{p}_{\delta_{0}^{-}+2,\delta_{1}+2}({\mathbb{R}}^{n}_{o})$. But $\chi\not\in
W^{2,p}_{\delta_{0}^{-},\delta_{1}}({\mathbb{R}}^{n}_{o})$ since
$\delta_{0}^{-}<-n/p$, so ${\mathcal{L}}\chi$ is not in the image of (61), and
hence $\zeta[{\mathcal{L}}(\chi)]\not=0$. This enables us to find $c_{1}$ so
that

$\zeta[{\mathcal{L}}(\chi v)]=c_{1}\zeta[{\mathcal{L}}\chi].$

But this means $\zeta[{\mathcal{L}}(\chi v-c_{1}\chi)]=0$, i.e.
${\mathcal{L}}(\chi v-c_{1}\chi)$ is in the image of (61), so
${\mathcal{L}}w={\mathcal{L}}(\chi v-c_{1}\chi)$ for some $w\in
W^{2,p}_{\delta^{-}_{0},\delta_{1}}({\mathbb{R}}^{n}_{o})\subset
W^{2,p}_{\delta_{0},\delta_{1}}({\mathbb{R}}^{n}_{o})$. Since
$\delta_{0}\in(-n/p,-2+n/p^{\prime})$, by the isomorphism (i) we have $w=\chi
v-c_{1}\chi$. In other words, for $|x|<1/2$ we have $v=c_{1}+w$, where $w\in
W^{2,p}_{\delta^{-}_{0},\delta_{1}}({\mathbb{R}}^{n}_{o})$. For any
$\varepsilon_{1}\in(0,1)$ we can let $\delta_{0}^{-}=\varepsilon_{1}-1-n/p$
and conclude $M_{2,p}(w,r)\leq c\,r^{1-\varepsilon_{1}}$, as stated in Theorem
1. This completes the proof. $\Box$

Now let us combine Theorem 1 with a change of variables to treat a general
$y\in{\mathbb{R}}^{n}$ and we do not assume $a_{ij}(y)=\delta_{ij}$. We let
$B_{\varepsilon}(y)=\\{x:|x-y|<\varepsilon\\}$ and want to construct a
singular solution of

(62) ${\mathcal{L}}(x,\partial_{x})Z_{y}(x)=0\quad\hbox{for}\ x\in
B_{\varepsilon}(y)\backslash\\{y\\},$

provided $\varepsilon$ is sufficiently small. Since ${\bf A}_{y}$ is symmetric
and positive definite, we can define the symmetric matrix ${\bf
A}_{y}^{-1/2}$. This enables us to define the function $I_{y}(r)$ as in (30).

######  Corollary 1.

For $n\geq 3$, $p\in(1,\infty)$, and fixed $y\in U$, assume that ${\bf A}_{y}$
is positive definite and the coefficients $a_{ij}$ satisfy (29). If $y=0$, we
assume (8) but otherwise the $b_{k}$ are bounded on $U$. Then, for
$\varepsilon$ sufficiently small, there exists a solution $Z_{y}$ of (62) in
the form

(63a) $Z_{y}(x)=h_{y}(|{\bf A}_{y}^{-1/2}(x-y)|)+v(x),$ where (63b)
$h_{y}(r)=\int_{r}^{\varepsilon}s^{1-n}\,e^{I_{y}(s)}\,ds\,(1+\zeta_{y}(r))$
with $I_{y}(r)$ given by (30) and $\zeta_{y}$ satisfying (63c)
$M_{2,p}(\zeta_{y},r)\leq c\,\max(\omega(r),\sigma(r)),$ and $v$ satisfying
(63d) $M_{2,p}(v,r;y)\leq c\,r^{2-n}\,e^{I_{y}(r)}\,\omega(r).$

Moreover, for any $u\in
W^{2,p}_{\ell\mathrm{o}c}(\overline{B_{\varepsilon}(y)}\backslash\\{y\\})$
that is a strong solution of ${\mathcal{L}}(x,\partial_{x})u=0$ in
$\overline{B_{\varepsilon}(y)}\backslash\\{y\\}$ and subject to the growth
condition

$M_{2,p}(u,r;y)\leq c\,r^{1-n+\varepsilon_{0}}\quad\hbox{where}\
\varepsilon_{0}>0,$

there exist constants $c_{0},c_{1}$ such that

$u(x)=c_{0}\,Z_{y}(x)+c_{1}+w(x),$

where $M_{2,p}(w,r;y)\leq c\,r^{1-\varepsilon_{1}}$ for any
$\varepsilon_{1}>0$.

If we use integration by parts, we can write the solution of Corollary 1 as

(64) $Z_{y}(x)=\frac{\langle{\bf
A}_{y}^{-1}(x-y),(x-y)\rangle^{\frac{2-n}{2}}}{(n-2)}\,e^{I_{y}(\sqrt{\langle{\bf
A}_{y}^{-1}(x-y),(x-y)\rangle})}(1+\xi_{y}(x))$

where $M_{1,\infty}(\xi_{y},r;y)\leq c\,\max(\omega(r),\sigma(r))$ for
$0<r<\varepsilon$. This generalizes (20).

Proof of Corollary 1. First we consider $y=0$. Since ${\bf A}_{0}$ is positive
definite, we can define a symmetric matrix by ${\bf J}={\bf A}_{0}^{-1/2}$, so
${\bf J}\,{\bf A}_{0}\,{\bf J}={\bf I}$. If we introduce a change of variables
by $\tilde{x}={\bf J}\,x$ and new coefficients $\widetilde{a}_{ij}(\tilde{x})$
and $\widetilde{b}_{k}(\tilde{x})$ by

$\widetilde{\bf A}_{\tilde{x}}={\bf J}{\bf A}_{x}{\bf
J}\quad\hbox{and}\quad\widetilde{\bf B}_{\tilde{x}}={\bf B}_{x}\,{\bf J},$

then $\widetilde{\bf A}_{0}={\bf I}$ and

${\mathcal{L}}(x,\partial_{x})=a_{ij}\frac{\partial^{2}}{\partial
x_{i}\partial x_{j}}+b_{k}\frac{\partial}{\partial
x_{k}}=\widetilde{a}_{ij}\frac{\partial^{2}}{\partial\tilde{x}_{i}\partial\tilde{x}_{j}}+\widetilde{b}_{k}\frac{\partial}{\partial\tilde{x}_{k}}=\widetilde{\mathcal{L}}(\tilde{x},\partial_{\tilde{x}}).$

Hence we may apply Theorem 1 in the variables $\tilde{x}$ to conclude the
existence of a solution $\widetilde{Z}$ of
$\widetilde{\mathcal{L}}(\tilde{x},\partial_{\tilde{x}})\widetilde{Z}(\tilde{x})=0$
for $0<|\tilde{x}|<\varepsilon$ in the form
$\widetilde{Z}(\tilde{x})+\widetilde{h}(|\tilde{x}|)+\widetilde{v}(\tilde{x})$
where

$\widetilde{h}(r)=\int_{r}^{\varepsilon}s^{1-n}\,e^{\widetilde{I}(s)}\,ds\,(1+\zeta(r))$

with $M_{2,p}(\zeta,r)\leq c\,\max(\omega(r),\sigma(r))$,
$M_{2,p}(\widetilde{v},r)\leq c\,r^{2-n}\,e^{\widetilde{I}(r)}\,\omega(r)$,
and

$\widetilde{I}(r)=\frac{1}{|S^{n-1}|}\int_{r<|\tilde{x}|<\varepsilon}\left({\rm
tr}(\widetilde{\bf A}_{z})-n\,\frac{\langle\widetilde{\bf
A}_{z}z,z\rangle}{|z|^{2}}+\langle\widetilde{\bf
B}_{z},z\rangle\right)\,\frac{dz}{|z|^{n}}.$

In terms of the original variables, we see that $Z(x)=\widetilde{h}(|{\bf
J}x|)+\widetilde{v}({\bf J}x)$ satisfies ${\mathcal{L}}(x,\partial_{x})Z=0$
for $0<|{\bf J}x|<\varepsilon$, hence for $0<|x|<\varepsilon_{1}$ with
$\varepsilon_{1}$ sufficiently small.

Finally, if $y$ is a general point in $U$, then we use the change of variables
$\tilde{x}={\bf J}(x-y)$ with ${\bf J}={\bf A}_{y}^{-1/2}$ and let
$\widetilde{\bf A}_{\tilde{x}}=(\widetilde{a}_{ij}(\tilde{x}))={\bf J}{\bf
A}_{x}{\bf J}.$ Since $\tilde{x}=0$ corresponds to $x=y$, we have
$\widetilde{a}_{ij}(0)=\delta_{ij}$, so we can apply Theorem 1 to
$\widetilde{\mathcal{L}}(\tilde{x},\partial_{\tilde{x}})={\mathcal{L}}(x,\partial_{x})$
to obtain the solution $\widetilde{h}(|\tilde{x}|)+\widetilde{v}(\tilde{x})$.
We obtain the solution of (62) as

$Z_{y}(x)=\widetilde{h}(|{\bf J}(x-y)|)+\widetilde{v}({\bf J}(x-y)),$

where $\widetilde{h}(r)$ involves the above $\widetilde{I}(r)$. To transform
$\widetilde{I}(r)$ to the original variables, replace $\widetilde{\bf
A}_{\tilde{z}}$ by ${\bf A}_{z}$ and every occurence of $\tilde{z}$ by $x-y$;
we find $\widetilde{I}$ is of the desired form (30). Moreover, since
$\widetilde{v}$ satisfies $M_{2,p}(\widetilde{v},r)\leq
c\,r^{2-n}e^{I(r)}\omega(r)$, we find that $v(x)=\widetilde{v}({\bf J}(x-y))$
satisfies $M_{2,p}(v,r;y)\leq c\,r^{2-n}\,e^{I(r)}\omega(r)$, as desired.
$\Box$

##  2\. Finding the constant $C_{y}$ so that
$-{\mathcal{L}}Z(x)=C_{y}\delta(x-y)$

In this section we first prove Theorem 2, then state and prove Corollary 2,
which shows how the formulas change when we no longer assume $y=0$ and
$a_{ij}(y)=\delta_{ij}$. As in the proof of Theorem 1, we shall assume (34)
holds for $\mu$ sufficiently small and work in $B_{1}$ instead of
$B_{\varepsilon}$. So we want to determine when the singular solution $Z(x)$
of Theorem 1 satisfies

(65) $-{\mathcal{L}}(x,\partial_{x})\,Z(x)=C_{0}\,\delta(x)\quad\hbox{ in
$U=B_{1}$,}$

for some constant $C_{0}$. This was done for ${\mathcal{L}}_{0}$ in [11], and
many of the arguments here are the same, so we shall be brief; for more
details, consult [11].

From Theorem 1, we obtain estimates on $\partial Z$ and $\partial^{2}Z$ that
show ${\mathcal{L}}(x,\partial_{x})Z(x)$ can be regularized at $x=0$ as a
distribution ${\mathcal{F}}_{0}$ on $C^{\lambda}_{0}(U)$, the Hölder
continuous functions with compact support in $U$. So if we can define
${\mathcal{L}}(x,\partial_{x})Z(x)$ as a distribution ${\mathcal{F}}$, then it
must be supported at $x=0$ and only involve $\delta(x)$, not derivatives of
$\delta(x)$, i.e. satisfy (65) for some constant $C_{0}$.

The difficulty in defining ${\mathcal{L}}(x,\partial_{x})Z(x)$ as a
distribution comes from the lack of regularity of the coefficients, especially
$a_{ij}$. In particular, there is no difficulty in defining 2nd-order
distributional derivatives of $Z$ by

$\langle\partial_{i}\partial_{j}Z,\phi\rangle=-\int_{U}\partial_{j}Z(x)\,\partial_{i}\phi(x)\,dx\quad\hbox{for}\
\phi\in C^{1}_{0}(U),$

since the integral on the right converges. So let us try to define the
distribution ${\mathcal{L}}Z$ by

(66)
$\langle{\mathcal{L}}Z,\phi\rangle\\!=\\!\int_{U}\\!\left((a_{ij}-\delta_{ij})\partial_{i}\partial_{j}Z\phi-\partial_{i}Z\,\partial_{i}\phi+b_{k}\partial_{k}Z\,\phi\right)dx\quad\hbox{for}\
\phi\in C^{1}_{0}(U).$

This is an improper integral due to the singularities in
$\partial_{i}\partial_{j}Z$ and $b_{k}\partial_{k}Z$ at $x=0$, but if the
integral converges then we conclude (65) holds and we can compute $C_{0}$ from

(67) $-C_{0}=\lim_{\varepsilon\to
0}\int_{U}\left((a_{ij}-\delta_{ij})\,\partial_{i}\partial_{j}Z\,\phi_{\varepsilon}-\partial_{i}Z\,\partial_{i}\phi_{\varepsilon}+b_{k}\partial_{k}Z\,\phi_{\varepsilon}\right)dx,$

where $\phi_{\varepsilon}(|x|)=\chi(|x|/\varepsilon)$ with $\chi(r)$ being a
smooth cutoff function that is 1 for $0<r<1/4$ and vanishes for $r>1/2$. (We
may assume $\phi(x)=\phi(|x|)$ since we can write
$\phi(x)=\phi_{0}(|x|)+\phi_{1}(x)$ where
$|\phi_{1}(x)|+|x|\,|\nabla\phi_{1}(x)|\leq c\,|x|$ for $|x|<1$, which shows
that $\langle{\mathcal{L}}Z,\phi_{1}\rangle$ is well-defined as an integral
and contributes nothing to $C_{0}$.)

Proof of Theorem 2. Recall from the proof of Theorem 1 the decomposition
$Z(x)=h(|x|)+v(x)$ in (37a). We first show that $v$ makes no contribution to
determining the value of $C_{0}$. Since $v\in Y$, we have $M_{2,p}(v,r)\leq
c\,r^{2-n}\,\omega(r)\,e^{I(r)}$. But $I(r)$ is bounded above in both cases
(i) and (ii), so we have

(68) $M_{2,p}(v,r)\leq c\,r^{2-n}\,\omega(r)\quad\hbox{for}\ 0<r<1.$

Then, as $\varepsilon\to 0$, we have111In the following, integrals involving
$\int_{|x|<\varepsilon}$ should be interpreted as improper: $\lim_{\eta\to
0}\int_{\eta<|x|<\varepsilon}$.

$\left|\int_{|x|<\varepsilon}(a_{ij}-\delta_{ij})\,\partial_{i}\partial_{j}v\,\phi_{\varepsilon}\,dx\right|\leq
c\int_{0}^{\varepsilon}\omega^{2}(r)\frac{dr}{r}=c\,\sigma(\varepsilon)\to 0$
$\left|\int_{|x|<\varepsilon}\partial_{i}v\,\partial_{i}\phi_{\varepsilon}\,dx\right|\leq
c\,\varepsilon^{-1}\int_{0}^{\varepsilon}\omega(r)\,dr\leq
c\,\omega(\varepsilon)\to 0$

and

$\left|\int_{|x|<\varepsilon}b_{k}\,\partial_{k}v\,\phi_{\varepsilon}\,dx\right|\leq
c\int_{0}^{\varepsilon}\omega^{2}(r)\frac{dr}{r}=c\,\sigma(\varepsilon)\to 0.$

So $v$ makes no contribution to $C_{0}$.

Now we consider $h(r)$. In fact, from (56) we can write
$h(r)=h_{0}(r)+h_{1}(r)+c$, where $c$ is an arbitrary constant and

(69) $\displaystyle h_{0}(r)$
$\displaystyle=c_{*}\int_{r}^{1}s^{1-n}\,E_{+}(s)\,ds,$ $\displaystyle
h_{1}(r)$
$\displaystyle=\int_{r}^{1}s^{1-n}\,E_{+}(s)\int_{0}^{s}\rho^{n-1}\,E_{-}(\rho)\,B[v](\rho)\,d\rho\,ds,$

with $E_{+}$ defined in (42) and $c_{*}$ chosen as in the proof of Theorem 1
so that $c_{*}E_{+}(0)=e^{I(0)}$. (Note that this decomposition of $h(r)$ is
slightly different from (57).) Let us show that $h_{1}$ and $c$ do not
contribute to $C_{0}$. We compute

$\partial_{i}h_{1}=-x_{i}\,r^{-n}\,E_{+}(r)\int_{0}^{r}\rho^{n-1}\,E_{-}(\rho)\,B[v](\rho)\,d\rho,$
$\displaystyle\partial_{i}\partial_{j}h_{1}\\!=\\!-r^{-n}\,E_{+}(r)\\!\left(\delta_{ij}-(n+R(r))\frac{x_{i}x_{j}}{r^{2}}\right)\\!\int_{0}^{r}\rho^{n-1}E_{-}(\rho)B[v](\rho)\,d\rho\\!-\\!\frac{x_{i}x_{j}}{r^{2}}B[v](r),$

and hence

$\displaystyle\int_{|x|<\varepsilon}\left[(a_{ij}-\delta_{ij})\partial_{i}\partial_{j}h_{1}\phi_{\varepsilon}+b_{k}\partial_{k}h_{1}\phi_{\varepsilon}-\partial_{i}h_{1}\partial_{i}\phi_{\varepsilon}\right]dx$
$\displaystyle=\int_{0}^{\varepsilon}\left[-\frac{R(r)}{r}E_{+}(r)\chi\left(\frac{r}{\varepsilon}\right)+E_{+}\frac{d}{dr}\chi\left(\frac{r}{\varepsilon}\right)\right]$
$\displaystyle\int_{0}^{r}\rho^{n-1}E_{-}(\rho)B[v](\rho)\,d\rho\,dr$
$\displaystyle-\int_{0}^{\varepsilon}(\alpha(r)-1)r^{n-1}B[v](r)\chi\left(\frac{r}{\varepsilon}\right)\,dr$
$\displaystyle+\,\varepsilon^{-1}\int_{0}^{\varepsilon}E_{+}(r)\int_{0}^{r}\rho^{n-1}E_{-}(\rho)B[v](\rho)$
$\displaystyle\,d\rho\,\chi\left(\frac{r}{\varepsilon}\right)dr.$

Let us denote these three terms ${\mathcal{I}}_{1}(\varepsilon)$,
${\mathcal{I}}_{2}(\varepsilon)$, ${\mathcal{I}}_{3}(\varepsilon)$ and
estimate them separately. First,

$\displaystyle{\mathcal{I}}_{1}(\varepsilon)\\!=$
$\displaystyle\\!\int_{0}^{\varepsilon}\frac{d}{dr}\left[E_{+}(r)\chi\left(\frac{r}{\varepsilon}\right)\int_{0}^{r}\rho^{n-1}E_{-}(\rho)B[v](\rho)\,d\rho\right]dr\\!-\\!\int_{0}^{\varepsilon}r^{n-1}\chi\left(\frac{r}{\varepsilon}\right)B[v](r)\,dr$
$\displaystyle=E_{+}(\varepsilon)\int_{0}^{e}r^{n-1}E_{-}(r)B[v](r)\,dr-\int_{0}^{\varepsilon}r^{n-1}B[v](r)\chi\left(\frac{r}{\varepsilon}\right)\,dr.$

But recall from (50) that
$\int_{0}^{\varepsilon}r^{n-1}E_{-}(r)|B[v](r)|\,dr\leq
c\,\sigma(\varepsilon)$ and $I(r)$ is bounded implies $E_{\pm}(r)$ are bounded
above and below by positive constants, so
$|{\mathcal{I}}_{1}(\varepsilon)|\leq c\,\sigma(\varepsilon)\to 0$ as
$\varepsilon\to 0$. This also implies $|{\mathcal{I}}_{2}(\varepsilon)|\to 0$
as $\varepsilon\to 0$. As for ${\mathcal{I}}_{3}(\varepsilon)$, we also have

$|{\mathcal{I}}_{3}(\varepsilon)|\leq
c\,\varepsilon^{-1}\int_{0}^{\varepsilon}\sigma(r)\,dr\leq
c\,\sigma(\varepsilon)\to 0\quad\hbox{as}\ \varepsilon\to 0.$

We conclude that $h_{1}$ does not contribute to $C_{0}$.

Finally, we consider $h_{0}$. We can calculate
$\partial_{i}\,h_{0}=-c_{*}\,r^{-n}\,E_{+}(r)\,x_{i}$, and

$\partial_{i}\partial_{j}h_{0}=-c_{*}\,r^{-n}\,E_{+}(r)\left(\delta_{ij}-n\frac{x_{i}x_{j}}{r^{2}}-\frac{x_{i}x_{j}}{r^{2}}R(r)\right).$

It is easy to verify that

$(a_{ij}-\delta_{ij})\partial_{i}\partial_{j}h_{0}=-c_{*}\,r^{-n}\,E_{+}(r)\left(a_{ii}-n\frac{a_{ij}x_{i}x_{j}}{r^{2}}-\left(\frac{a_{ij}x_{i}x_{j}}{r^{2}}-1\right)R(r)\right).$

Notice that

$\int_{|x|<\varepsilon}r^{-n}\,E_{+}(r)\left(a_{ii}-n\frac{a_{ij}x_{i}x_{j}}{r^{2}}-\left(\frac{a_{ij}x_{i}x_{j}}{r^{2}}-1\right)R(r)\right)\phi_{\varepsilon}(|x|)\,dx$
$=|S^{n-1}|\int_{0}^{\varepsilon}E_{+}(r)\frac{\alpha_{n}(r)-n\alpha(r)-(\alpha(r)-1)R(r)}{r}\chi\left(\frac{r}{\varepsilon}\right)\,dr$
$=|S^{n-1}|\int_{0}^{\varepsilon}E_{+}(r)\frac{R(r)-\beta(r)}{r}\chi\left(\frac{r}{\varepsilon}\right)\,dr,$

since $\alpha(r)\,R(r)=\alpha_{n}(r)+\beta(r)-n\alpha(r)$. Similarly, we can
verify

$\int_{|x|<\varepsilon}b_{k}\partial_{k}h_{0}\,\phi_{\varepsilon}\,dx=-c_{*}\,|S^{n-1}|\int_{0}^{\varepsilon}E_{+}(r)\,\frac{\beta(r)}{r}\,\chi\left(\frac{r}{\varepsilon}\right)\,dr,$

so

$\langle{\mathcal{L}}h_{0},\phi_{\varepsilon}\rangle=c_{*}\,|S^{n-1}|\int_{0}^{\varepsilon}\left(-E_{+}(r)\frac{R(r)}{r}\chi\left(\frac{r}{\varepsilon}\right)+E_{+}(r)\frac{d}{dr}\left[\chi\left(\frac{r}{\varepsilon}\right)\right]\right)\,dr.$

Hence

(70)
$-\langle{\mathcal{L}}h_{0},\phi_{\varepsilon}\rangle=-c_{*}|S^{n-1}|\,\int_{0}^{\varepsilon}\frac{d}{dr}\left[E_{+}(r)\chi\left(\frac{r}{\varepsilon}\right)\right]\,dr=c_{*}\,|S^{n-1}|\,E_{+}(0).$

Letting $\varepsilon\to 0$, we conclude $Z$ satisfies (65) with
$C_{0}=|S^{n-1}|e^{I(0)}$. $\Box$

Now let us consider a general fixed $y\in U$ and try to solve

(71)
$-{\mathcal{L}}(x,\partial_{x})Z_{y}(x)=C_{y}\,\delta(x-y)\quad\hbox{for}\
x\in B_{1}(y)$

for some constant $C_{y}$. We replace (66) with

(72)
$\langle{\mathcal{L}}Z_{y},\phi\rangle\\!=\\!\int_{U}\\!\left((a_{ij}-a_{ij}(y))\partial_{i}\partial_{j}Z_{y}\phi-
a_{ij}(y)\,\partial_{i}Z_{y}\,\partial_{j}\phi+b_{k}\partial_{k}Z_{y}\,\phi\right)dx.$

######  Corollary 2.

Suppose the conditions of Corollary 1 hold and $Z_{y}$ is the singular
solution found there with $I_{y}(r)$ given by (30a).

  * (i)

If $I_{y}(0)=\lim_{r\to 0}I_{y}(r)$ exists and is finite, then we can solve
(71) in $B_{1}(y)$ with $C_{y}=|S^{n-1}|\,\sqrt{{\rm det}{\bf
A}_{y}}\,e^{I_{y}(0)}$.

  * (ii)

If $I_{y}(r)\to-\infty$ as $r\to 0$, then solving (10) in $B_{1}(y)$ yields
$C_{0}=0$, and so $Z$ solves $-{\mathcal{L}}(x,\partial_{x})Z_{y}(x)=0$ in
$B_{1}(y)$, despite its singularity at $x=y$.

If $I_{y}(0)$ exists and is finite, we see from (64) that the solution $F_{y}$
of $-{\mathcal{L}}(x,\partial_{x})F_{y}(x)=\delta(x-y)$ in
$B_{\varepsilon}(y)$ has the asymptotic behavior

(73) $F_{y}(x)=\frac{\langle{\bf
A}_{y}^{-1}(x-y),(x-y)\rangle^{\frac{2-n}{2}}}{(n-2)|S^{n-1}|\sqrt{{\rm
det}{\bf A}_{y}}}\,(1+o(1))\quad\hbox{as}\ x\to y.$

This generalizes (21) and establishes (31).

Proof of Corollary 2. We only need to show that

$\langle-{\mathcal{L}}Z_{y},\phi\rangle=|S^{n-1}|\,({\rm det}{\bf
A}_{y})^{-1/2}\,e^{I_{y}(0)}\phi(y)$

for some $\phi\in C^{\infty}_{0}(B_{\varepsilon_{y}}(y))$. Let us recall the
change of coordinates used in the proof of Theorem 1, namely $\tilde{x}={\bf
J}(x-y)$ where ${\bf J}={\bf A}_{y}^{-1/2}$, and let
$\widetilde{\phi}(\tilde{x})=\phi(x)$. Then

$-\int{\mathcal{L}}(x,\partial_{x})Z_{y}(x)\phi(x)\,dx=-({\det}{\bf
A}_{y})^{1/2}\int\widetilde{\mathcal{L}}(\tilde{x},\partial_{\tilde{x}})\widetilde{Z}_{0}(\tilde{x})\widetilde{\phi}(\tilde{x})\,d\tilde{x}.$

But Theorem 2 implies
$-\langle\widetilde{\mathcal{L}}\widetilde{Z}_{0},\widetilde{\phi}\rangle=|S^{n-1}|\,e^{I_{y}(0)}\widetilde{\phi}(0)$.
Since $\widetilde{\phi}(0)=\phi(y)$, we obtain the desired result. $\Box$

##  3\. Constructing the fundamental solution

Proof of Theorem 3. For each $y\in U$, denote the $\varepsilon$ in Corollary 1
by $\varepsilon_{y}$, and use Corollary 2 to calculate $C_{y}$, which is
positive since $I_{y}(0)$ is finite. We conclude that
$-{\mathcal{L}}(x,\partial_{x})Z_{y}(x)/C_{y}=\delta(x-y)$ for all $x,y\in U$
with $|x-y|<\varepsilon$. For fixed $y\in U$, let $\eta_{y}(r)$ be a smooth
cutoff function satisfying $\eta_{y}(r)=1$ for sufficiently small $r$ and
define

(74) $F(x,y)=\eta_{y}(|x-y|)Z_{y}(x)/C_{y}+v(x,y),$

where $v(x,y)$ is to be determined. But if we apply
$-{\mathcal{L}}(x,\partial_{x})$ we obtain

$-{\mathcal{L}}(x,\partial_{x})F(x,y)=\delta(x-y)+\psi(x,y)-{\mathcal{L}}(x,\partial_{x})v(x,y),$

where $\psi(\cdot,y)\in L^{p}(U)$ for $1<p<\infty$ and vanishes near $y$. We
may assume that $U$ has a smooth boundary (else we can embed it in such a
bounded domain and extend the coefficients, see [16]), so we may find $v(x,y)$
by solving the Dirichlet problem (for fixed $y\in U$):

(75) $\displaystyle-{\mathcal{L}}(x,\partial_{x})v(x,y)$
$\displaystyle=\psi(x,y)\quad\hbox{for}\ x\in U,$ $\displaystyle v(x,y)$
$\displaystyle=0\quad\hbox{for}\ x\in\partial U.$

When $b_{k}\in L^{\infty}(U)$, it is well-known (cf. Theorem 9.15 in [6]) that
(75) has a unique solution $v(\cdot,y)\in W^{2,p}(U)\cap W^{1,p}_{0}(U)$ for
$1<p<\infty$; associated with this unique solvability is the à priori
inequality $\|v\|_{W^{2,p}(U)}\leq c\,\|{\mathcal{L}}v\|_{L^{p}(U)}$ for $v\in
W^{2,p}(U)\cap W_{0}^{1,p}(U)$. But our assumption (8) enables us to write
${\mathcal{L}}={\mathcal{L}}_{1}+{\mathcal{L}}_{2}$, where ${\mathcal{L}}_{1}$
has bounded coefficients in $U$ and
${\mathcal{L}}_{2}=\tilde{b}_{k}\partial_{k}$ has coefficients supported in a
ball $B_{\sigma}$. For any $\varepsilon>0$, we can take $\sigma$ sufficiently
small that $\sup_{x\in
B_{\sigma}}|\tilde{b}_{k}(x)|\leq\varepsilon\,|x|^{-1}$. For $1<p<n$, by
Hardy’s inequality we have
$\|{\mathcal{L}}_{2}w\|_{L^{p}}\leq\varepsilon\,\|w\|_{W^{2,p}}$ for $w\in
W^{2,p}$ with support in $B_{\sigma}$. If we take $\varepsilon$ sufficiently
small, we can arrange that the à priori inequality $\|v\|_{W^{2,p}(U)}\leq
c\,\|{\mathcal{L}}v\|_{L^{p}(U)}$ holds for $v\in W^{2,p}(U)\cap
W_{0}^{1,p}(U)$. We conclude that (75) admits a unique solution $v(\cdot,y)\in
W^{2,p}(U)\cap W_{0}^{1,p}(U)$ for $1<p<n$. Taking $p\in(n/2,n)$, we see that
$v(x,y)$ is continuous in $x\in U$.

To confirm that $F(x,y)$ satisfies (31), we know from (74) and (64) that

$F(x,y)=\frac{\langle{\bf
A}_{y}^{-1}(x-y),(x-y)\rangle^{\frac{2-n}{2}}}{(n-2)|S^{n-1}|\sqrt{\det{\bf
A}_{y}}}\,e^{I_{y}(\sqrt{\langle{\bf
A}_{y}^{-1}(x-y),(x-y)\rangle})-I_{y}(0)}(1+\xi_{y}(x))$

where $M_{1,\infty}(\xi_{y},r;y)\leq c\,\max(\omega(r),\sigma(r))$ for
$0<r<\varepsilon$. Since the exponential term tends to 1 as $|x-y|\to 0$, we
obtain (31). $\Box$

##  Appendix A $L^{p}$-mean estimates for convolutions

In [11], we proved $L^{p}$-mean estimates for distribution solutions of

(76) $\Delta u=f\quad\hbox{in}\ {\mathbb{R}}^{n}\backslash\\{0\\},$

when $f$ has zero spherical mean, i.e. $\overline{f}(r)=0$. Let $K$ denote
convolution by the fundamental solution $\Gamma$ of the Laplacian. The
following appears as Corollary 1 in Section 1 of [11].

######  Proposition 1.

Suppose $n\geq 2$, $p\in(1,\infty)$, and $f\in
L^{p}_{{\ell\mathrm{o}c}}({\mathbb{R}}^{n}\backslash\\{0\\})$ satisfies

(77)
$\overline{f}(r)=0\quad\hbox{and}\quad\int_{|x|<1}|x|\,|f(x)|\,dx+\int_{|x|>1}|x|^{1-n}|f(x)|\,dx<\infty.$

Then $u=Kf=\Gamma\star f$ is a distribution solution of (76) that satisfies

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

##  Appendix B The Laplacian on weighted Sobolev spaces

Many authors have studied the mapping properties of elliptic operators such as
the Laplacian on weighted Sobolev spaces on ${\mathbb{R}}^{n}$ and other
noncompact manifolds with conical or cylindrical ends: cf. [10], [14], [15].
We will recall some of these results for punctured Euclidean space
${\mathbb{R}}^{n}_{o}:={\mathbb{R}}^{n}\backslash\\{0\\}$. Since we are mostly
concerned in this paper with singularities at $x=0$, let us first investigate
the weighted $L^{p}$-norm on the punctured unit ball
$B_{o}:=B_{1}\backslash\\{0\\}$. For $\delta\in{\mathbb{R}}$ and $1<p<\infty$,
define the Banach space $L^{p}_{\delta}(B_{o})$ by the norm

(79) $\|u\|^{p}_{L^{p}_{\delta}(B_{o})}:=\int_{|x|<1}|x|^{\delta
p}\,|u(x)|^{p}\,dx.$

For example, the constants are in $L^{p}_{\delta}(B_{o})$ if and only if
$\delta>-n/p$. To compare (79) with $L^{p}$-means, it is easy to see (cf.
[11]) that $M(u,r)\leq c\,r^{\varepsilon}$ for $0<r<1$ implies that $u\in
L_{\delta}^{p}(B_{o})$ provided $\varepsilon+\delta>-n/p$, and conversely
$u\in L_{\delta}^{p}(B_{o})$ implies $M_{p}(u,r)\leq
c_{\varepsilon}\,r^{\varepsilon}$ for $0<r<1$ if $\varepsilon=-\delta-n/p$.

Now we introduce the weighted $L^{p}$-norm for functions on
${\mathbb{R}}^{n}_{o}$ with separate weights at the origin and at infinity.
For $\delta_{0},\delta_{1}\in{\mathbb{R}}$, define

(80)
$\displaystyle\|u\|^{p}_{L^{p}_{\delta_{0},\delta_{1}}}:=\int_{|x|<1}|x|^{\delta_{0}p}|u(x)|^{p}\,dx+\int_{|x|>1}|x|^{\delta_{1}p}|u(x)|^{p}\,dx.$

We then define the weighted Sobolev space
$W^{2,p}_{\delta_{o},\delta_{1}}({\mathbb{R}}^{n}_{o})$ to be functions $u\in
W^{2,p}_{{\ell\mathrm{o}c}}({\mathbb{R}}^{n}_{o})$ for which

(81) $\|u\|_{W^{2,p}_{\delta_{o},\delta_{1}}}=\sum_{|\alpha|\leq
2}\||x|^{|\alpha|}\partial^{\alpha}u\|_{L^{p}_{\delta_{1},\delta_{2}}}$

is finite. It is clear that

(82) $\Delta:W^{2,p}_{\delta_{0},\delta_{1}}({\mathbb{R}}^{n}_{o})\to
L^{p}_{\delta_{0}+2,\delta_{1}+2}({\mathbb{R}}^{n}_{o})$

is a bounded linear operator, and (using the analysis of [10], [14], [15], for
example) it can be shown that (82) is an isomorphism for
$-n/p<\delta_{0},\delta_{1}<-2+n/p^{\prime}$, where $p^{\prime}=p/(p-1)$.
(Since $n\geq 3$, such $\delta_{0},\delta_{1}$ exist.) Moreover, provided
$\delta_{0},\delta_{1}$ do not take the values $-n/p-k$ or $-2+n/p^{\prime}+k$
where $k$ is a nonnegative integer, then (82) is a Fredholm operator whose
nullspace and cokernel are easily described in terms of harmonic polynomials.
Since we are principally interested in the behavior of functions at the
origin, we will fix

(83) $-n/p<\delta_{1}<-2+n/p^{\prime}$

and allow $\delta_{0}$ to vary. We only require a small range of values for
$\delta_{0}$.

######  Proposition 2.

Assume $1<p<\infty$ and (83). Then the map (82) is

  * (a)

an isomorphism for $-n/p<\delta_{0}<-2+n/p^{\prime}$;

  * (b)

surjective with nullspace spanned by $|x|^{2-n}$ if
$-2+n/p^{\prime}<\delta_{0}<-1+n/p^{\prime}$;

  * (c)

injective with cokernel spanned by $1$ if $-1-n/p<\delta_{0}<-n/p$.

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

## References

  * [1] R. Alvarado, D. Brigham, V. Maz’ya, M. Mitrea, E. Ziadé, Sharp geometric maximum principles for semi-elliptic operators with singular drift, Math. Res. Lett. 18 (2011), 10001-10007.
  * [2] M. Cranston, Z. Zhao, Conditional transformation of drift formula and potential theory for $\frac{1}{2}\Delta+b(\cdot)\cdot\nabla$, Comm. Math. Phys. 112 (1987), 613-625.
  * [3] H. Dong, S. Kim, On C1, C2, and weak type (1,1) estimates for linear elliptic operators, Comm. Partial Differential Equations, 42 (2017), 417-435.
  * [4] H. Dong, S. Kim, S. Lee, Note on Green’s functions of non-divergence elliptic operators with continuous coefficients, arXiv:2201.03764 .
  * [5] D. Gilbarg, J. Serrin, On isolated singularities of solutions of second-order elliptic equations, J. Analyse Math.  4 (1955/56), 309-340.
  * [6] D. Gilbarg, N. Trudinger, Elliptic Partial Differential Equations of Second Order, 2nd ed, Springer-Verlag, New York, 1983.
  * [7] S. Hwang, S. Kim, Green’s function for second order elliptic equations in non-divergence form, Potential Anal. 52 (2020), no. 1, 27–39.
  * [8] S. Kim, S. Lee, Estimates for Green’s functions of elliptic equations in non-divergence form with continuous coefficients,Ann. Appl. Math. 37 (2021), no. 2, 111–130.
  * [9] V.A. Kozlov, V.G. Maz’ya, Differential Equations with Operator Coefficients, Springer-Verlag, New York, 1999.
  * [10] R. Lockhart, R. McOwen, Elliptic differential operators on noncompact manifolds, Ann. Sc. Norm. Super. Pisa Cl. Sci. (4) 12 (1985), 409-447.
  * [11] V.G. Maz’ya, R. McOwen, On the fundamental solution of an elliptic equation in nondivergence form, AMS Translations: special volume dedicated to Nina Uraltseva 229 (2010), 145-172.
  * [12] V.G. Maz’ya, R. McOwen, Differentiablilty of solutions to second-order elliptic equations via dynamical systems, J. Differential Equations, 250 (2010), 1137-1168.
  * [13] V.G. Maz’ya, R. McOwen, Gilbarg-Serrin equation and Lipschitz regularity, J. Differential Equations 312 (2022), 45-64.
  * [14] V. Maz’ya, B. Plamenevski, Estimates in $L_{p}$ and Hölder classes and the Miranda-Agmon maximum principle solutions of elliptic boundary problems in domains with singular points on the boundary (Russian), Math. Nachr. 81 (1978), 25-82.
  * [15] R. McOwen, The behavior of the Laplacian on weighted Sobolev spaces, Comm. Pure Appl. Math. 32 (1979), 783-795.
  * [16] V. Mil’man, Extension of functions that preserve the modulus of continuity, Math. Notes 61 (1997), 193-200.
  * [17] C. Miranda, Partial Differential Equations of Elliptic Type, Springer-Verlag, Berlin, 1970.