111The research of J. Wei is partially supported by NSERC of Canada. Lei Zhang
is partially supported by a Simons Foundation Collaboration Grant

# Laplacian Vanishing Theorem for Quantized Singular Liouville Equation

Juncheng Wei  Department of Mathematics
University of British Columbia
Vancouver, BC V6T1Z2, Canada<EMAIL_ADDRESS>and  Lei Zhang  Department of
Mathematics
University of Florida
1400 Stadium Rd
Gainesville FL 32611<EMAIL_ADDRESS>
###### Abstract.

In this article we establish a vanishing theorem for singular Liouville
equation with quantized singular source. If a blowup sequence tends to
infinity near a quantized singular source and the blowup solutions violate the
spherical Harnack inequality around the singular source (non-simple blow-ups),
the Laplacian of a coefficient function must tend to zero. This seems to be
the first second order estimates for Liouville equation with quantized sources
and non-simple blow-ups. This result as well as the key ideas of the proof
would be extremely useful for various applications.

###### 1991 Mathematics Subject Classification:

35J75,35J61

##  1\. Introduction

This is the third article in our series to study blowup solutions of

(1.1) $\Delta u+|x|^{2N}\mathrm{H}(x)e^{u}=0,$

in a neighborhood of the origin in $\mathbb{R}^{2}$. Here $\mathrm{H}$ is a
positive smooth function and $N\in\mathbb{N}$ is a positive integer. Since the
analysis is local in nature we focus the discussion in a neighborhood of the
origin: Let $\mathfrak{u}_{k}$ be a sequence of solutions of

(1.2)
$\Delta\mathfrak{u}_{k}(x)+|x|^{2N}\mathrm{H}_{k}(x)e^{\mathfrak{u}_{k}}=0,\quad\mbox{in}\quad
B_{\tau}$

for some $\tau>0$ independent of $k$. $B_{\tau}$ is the ball centered at the
origin with radius $\tau$. In addition we postulate the usual assumptions on
$\mathfrak{u}_{k}$ and $\mathrm{H}_{k}$: For a positive constant $C$
independent of $k$, the following holds:

(1.3)
$\left\\{\begin{array}[]{ll}\|\mathrm{H}_{k}\|_{C^{3}(\bar{B}_{\tau})}\leq
C,\quad\frac{1}{C}\leq\mathrm{H}_{k}(x)\leq C,\quad x\in\bar{B}_{\tau},\\\ \\\
\int_{B_{\tau}}\mathrm{H}_{k}e^{\mathfrak{u}_{k}}\leq C,\\\ \\\
|\mathfrak{u}_{k}(x)-\mathfrak{u}_{k}(y)|\leq C,\quad\forall x,y\in\partial
B_{\tau},\end{array}\right.$

and since we study the asymptotic behavior of blowup solutions around the
singular source, we assume that there is no blowup point except at the origin:

(1.4) $\max_{K\subset\subset B_{\tau}\setminus\\{0\\}}\mathfrak{u}_{k}\leq
C(K).$

If a sequence of solutions $\\{u^{k}\\}_{k=1}^{\infty}$ of (1.1) satisfies

$\lim_{k\to\infty}u^{k}(x_{k})=\infty,\quad\mbox{ for some $\bar{x}\in
B_{\tau}$ and $x_{k}\to\bar{x}$,}$

we say $\\{u^{k}\\}$ is a sequence of bubbling solutions or blowup solutions,
$\bar{x}$ is called a blowup point. The question we consider in this work is
when $0$ is the only blowup point in a neighborhood of the origin, what
vanishing theorems will the coefficient functions $\mathrm{H}_{k}$ satisfy?

One indispensable assumption is that the blowup solutions violate the
spherical Harnack inequality around the origin:

(1.5) $\max_{x\in B_{\tau}}\mathfrak{u}_{k}(x)+2(1+N)\log|x|\to\infty,$

It is also mentioned in literature ( see [21, 26] ) that $0$ is called an non-
simple blowup point. The main result of this article is

######  Theorem 1.1.

Let $\\{\mathfrak{u}_{k}\\}$ be a sequence of solutions of (1.2) such that
(1.3),(1.4) hold and the spherical Harnack inequality is violated as in (1.5).
Then along a sub-sequence

$\lim_{k\to\infty}\Delta(\log\mathrm{H}_{k})(0)=0.$

Theorem 1.1 is a continuation of our previous result in [27]:

_Theorem A: Let $\\{\mathfrak{u}_{k}\\}$ be a sequence of solutions of (1.2)
such that (1.3),(1.4) and (1.5) hold. Then along a subsequence_

$\lim_{k\to\infty}\nabla(\log\mathrm{H}_{k}+\phi_{k})(0)=0$

_where $\phi_{k}$ is defined as_

(1.6) $\left\\{\begin{array}[]{ll}\Delta\phi_{k}(x)=0,\quad\mbox{in}\quad
B_{\tau},\\\ \\\
\phi_{k}(x)=\mathfrak{u}_{k}(x)-\frac{1}{2\pi\tau}\int_{\partial
B_{\tau}}\mathfrak{u}_{k}dS,\quad x\in\partial B_{\tau}.\end{array}\right.$

The equation (1.1) comes from its equivalent form

$\Delta v+\mathrm{H}e^{v}=4\pi N\delta_{0}$

by using a logarithmic function to eliminate the Dirac mass on the right hand
side. Since the strength of the Dirac mass is a multiple of $4\pi$, this type
of singularity is called “quantized”. An equation with a quantized singular
source is ubiquitous in literature. In particular the following mean field
equation defined on a Riemann surface $(M,g)$:

(1.7)
$\Delta_{g}u+\rho(\frac{h(x)e^{u(x)}}{\int_{M}he^{u}}-1)=4\pi\sum_{t=1}^{M}\alpha_{t}(\delta_{p_{t}}-1),$

represents a conformal metric with prescribed conic singularities (see [16,
24, 25]), where $h$ is a positive smooth function, $\rho>0$ is a constant and
the volume of $M$ is assumed to be $1$ for convenience, $\alpha_{j}>-1$ are
constants as well. If the singular source is quantized, the equation is
profoundly linked to Algebraic geometry, integrable system, number theory and
complex Monge-Ampere equations (see [13]). In Physics the main equation
reveals key features of mean field limits of point vortices in the Euler flow
[8, 9] and models in the Chern-Simons-Higgs theory [20] and in the electroweak
theory [2], etc.

So far the non-simple bubbling situation has been observed in Liouville
equation [21, 4], Liouville systems [18, 19, 28] and fourth order equations
[1]. The main theorem in this article would impact the study of these
equations as well as some well known open questions in Monge-Ampere equation
[26].

When compared with Theorem A, Theorem 1.1 is clearly more challenging in
analysis. In fact the proof of Theorem A is a special case of one step of the
proof of Theorem 1.1. However, their major difference is on applications.
Theorem 1.1 is significantly more influential for many reasons: First the main
motivation to study equation (1.1) is for equations or systems defined on
manifolds. Usually blowup analysis near a singular point needs to reflect the
curvature at the blowup point. In this respect Theorem 1.1 is directly related
to the Gauss curvature at the blowup point. Second, the harmonic function in
Theorem A causes inconvenience in application since it is generally hard to
identify what the harmonic function is. On the other hand Theorem 1.1 is only
involved with the Laplace of the coefficient function. This may lead to
substantial advances in applications: In many degree counting programs one
major difficulty is bubble-coalition, which means bubbling disks may collide
into one point. The formation of bubbling disks tending to one point is
accurately represented by equation (1.1). Theorem 1.1 and its proof could be
extremely useful to simplify blowup pictures. Third, the proof of Theorem 1.1
is also important for proving uniqueness of bubbling solutions, and the
results for Liouville equation with quantized singular sources is
inspirational to many equations and systems with similar singular poles.
Before our series of works most of the study of singular equations or systems
focuses on non-quantized singular situations. However it is the “quantized
situations” that manifest profound connections to different fields of
mathematics and Physics. Theorem 1.1 may be a starting point of multiple
directions of exciting adventures.

As a first application of Theorem 1.1 we present an advancement of the mean
field equation (1.7). Let $\Lambda$ be defined as

$\Lambda=\\{8\pi k+\sum_{j\in A}8\pi(1+\alpha_{j});\quad
k\in\mathbb{N}\cup\\{0\\},\quad A\subset\\{1,...,M\\}\quad\\}$

where $\mathbb{N}$ is the set of natural numbers. Through the works of
Bartoclucci-Tarantello [5, 6], Chen-Lin [13] etc, an a priori estimate holds
if $\rho\not\in\Lambda$. In other words, if $u^{k}$ is a sequence of blowup
solutions with parameters $\rho^{k}$, the limit of $\rho^{k}$ is in $\Lambda$.
Our second main theorem is

######  Theorem 1.2.

Let $u_{k}$ be a sequence of blowup solutions of (1.7) with parameters
$\rho^{k}\to\rho\in\Lambda$, $h$ be a positive smooth function,
$\alpha_{1},...,\alpha_{M}>-1$ are constants. If at each quantized blowup
point $p$ we have

(1.8) $\Delta\log h(p)-2K(p)-4\pi\sum_{t=1}^{M}\alpha_{t}+\rho\neq 0,$

where $K(p)$ is the Gauss curvature at $p$, all blowup points of $u_{k}$ are
simple blowup points.

The organization of the article is as follows. In section two we cite
preliminary results related to the proof of the main theorem. Then in section
three approximate the blowup solutions by a family of global solutions that
agree with the blowup solutions at one local maximum point. This is crucial
for our argument. Then we derive some intermediate estimates as preparation of
more precise analysis. In section four we prove the first order estimates that
cover the main result in [27]. This section proves stronger result than [27]
and provides more detail. Finally in section five we take advantage of the
result of the first order estimate and complete the proof of Theorem 1.1. The
proof of Theorem 1.2 is placed in section six. The final section is an
appendix that contains certain computations needed in the proof of Theorem
1.1.

Notation: We will use $B(x_{0},r)$ to denote a ball centered at $x_{0}$ with
radius $r$. If $x_{0}$ is the origin we use $B_{r}$. $C$ represents a positive
constant that may change from place to place.

##  2\. Preliminary discussions

In the first stage of the proof of Theorem 1.1 we set up some notations and
cite some preliminary results. For simple notation we set

(2.1) $u_{k}(x)=\mathfrak{u}_{k}(x)-\phi_{k}(x),\quad\mbox{and}$ (2.2)
$h_{k}(x)=\mathrm{H}_{k}(x)e^{\phi_{k}(x)}.$

to write the equation of $u_{k}$ as

(2.3) $\Delta u_{k}(x)+|x|^{2N}h_{k}(x)e^{u_{k}}=0,\quad\mbox{ in }\quad
B_{\tau}$

Without loss of generality we assume

(2.4) $\lim_{k\to\infty}h_{k}(0)=1.$

Obviously (1.5) is equivalent to

(2.5) $\max_{x\in B_{\tau}}u_{k}(x)+2(1+N)\log|x|\to\infty,$

It is well known [21, 4] that $u_{k}$ exhibits a non-simple blowup profile. It
is established in [21, 4] that there are $N+1$ local maximum points of
$u_{k}$: $p_{0}^{k}$,….,$p_{N}^{k}$ and they are evenly distributed on
$\mathbb{S}^{1}$ after scaling according to their magnitude: Suppose along a
subsequence

$\lim_{k\to\infty}p_{0}^{k}/|p_{0}^{k}|=e^{i\theta_{0}},$

then

$\lim_{k\to\infty}\frac{p_{l}^{k}}{|p_{0}^{k}|}=e^{i(\theta_{0}+\frac{2\pi
l}{N+1})},\quad l=1,...,N.$

For many reasons it is convenient to denote $|p_{0}^{k}|$ as $\delta_{k}$ and
define $\mu_{k}$ as follows:

(2.6) $\delta_{k}=|p_{0}^{k}|\quad\mbox{and
}\quad\mu_{k}=u_{k}(p_{0}^{k})+2(1+N)\log\delta_{k}.$

Also we use

$\epsilon_{k}=e^{-\frac{1}{2}\mu_{k}}$

to be the scaling factor most of the time. Since $p_{l}^{k}$’s are evenly
distributed around $\partial B_{\delta_{k}}$, standard results for Liouville
equations around a regular blowup point can be applied to have
$u_{k}(p_{l}^{k})=u_{k}(p_{0}^{k})+o(1)$. Also, (1.5) gives
$\mu_{k}\to\infty$. The interested readers may look into [21, 4] for more
detailed information.

Finally we shall use $E$ to denote a frequently appearing error term of the
size $O(\delta_{k}^{2})+O(\mu_{k}e^{-\mu_{k}})$.

##  3\. Approximating bubbling solutions by global solutions

We write $p_{0}^{k}$ as $p_{0}^{k}=\delta_{k}e^{i\theta_{k}}$ and define
$v_{k}$ as

(3.1)
$v_{k}(y)=u_{k}(\delta_{k}ye^{i\theta_{k}})+2(N+1)\log\delta_{k},\quad|y|<\tau\delta_{k}^{-1}.$

If we write out each component, (3.1) is

$v_{k}(y_{1},y_{2})=u_{k}(\delta_{k}(y_{1}\cos\theta_{k}-y_{2}\sin\theta_{k}),\delta_{k}(y_{1}\sin\theta_{k}+y_{2}\cos\theta_{k}))+2(1+N)\log\delta_{k}.$

Then it is standard to verify that $v_{k}$ solves

(3.2) $\Delta
v_{k}(y)+|y|^{2N}\mathfrak{h}_{k}(\delta_{k}y)e^{v_{k}(y)}=0,\quad|y|<\tau/\delta_{k},$

where

(3.3) $\mathfrak{h}_{k}(x)=h_{k}(xe^{i\theta_{k}}),\quad|x|<\tau.$

Thus the image of $p_{0}^{k}$ after scaling is $Q_{1}^{k}=e_{1}=(1,0)$. Let
$Q_{1}^{k}$, $Q_{2}^{k}$,…,$Q_{N}^{k}$ be the images of $p_{i}^{k}$
$(i=1,...,N)$ after the scaling:

$Q_{l}^{k}=\frac{p_{l}^{k}}{\delta_{k}}e^{-i\theta_{k}},\quad l=1,...,N.$

It is established by Kuo-Lin in [21] and independently by Bartolucci-
Tarantello in [4] that

(3.4)
$\lim_{k\to\infty}Q_{l}^{k}=\lim_{k\to\infty}p_{l}^{k}/\delta_{k}=e^{\frac{2l\pi
i}{N+1}},\quad l=0,....,N.$

Then it is proved in our previous work that ( see (3.13) in [26])

$Q_{l}^{k}-e^{\frac{2\pi
li}{N+1}}=O(\mu_{k}e^{-\mu_{k}})+O(|\nabla\log\mathfrak{h}_{k}(0)|\delta_{k}).$

Using the rate of $\nabla\mathfrak{h}_{k}(0)$ in [26] we have

(3.5) $Q_{l}^{k}-e^{\frac{2\pi
li}{N+1}}=O(\mu_{k}e^{-\mu_{k}})+O(\delta_{k}^{2}).$

Choosing $3\epsilon>0$ small and independent of $k$, we can make disks
centered at $Q_{l}^{k}$ with radius $3\epsilon$ (denoted as
$B(Q_{l}^{k},3\epsilon)$) mutually disjoint. Let

(3.6) $\mu_{k}=\max_{B(Q_{0}^{k},\epsilon)}v_{k}.$

Since $Q_{l}^{k}$ are evenly distributed around $\partial B_{1}$, it is easy
to use standard estimates for single Liouville equations ([30, 17, 12]) to
obtain

$\max_{B(Q_{l}^{k},\epsilon)}v_{k}=\mu_{k}+o(1),\quad l=1,...,N.$

Let

(3.7)
$V_{k}(x)=\log\frac{e^{\mu_{k}}}{(1+\frac{e^{\mu_{k}}\mathfrak{h}_{k}(\delta_{k}e_{1})}{8(1+N)^{2}}|y^{N+1}-e_{1}|^{2})^{2}}.$

Clearly $V_{k}$ is a solution of

(3.8) $\Delta
V_{k}+\mathfrak{h}_{k}(\delta_{k}e_{1})|y|^{2N}e^{V_{k}}=0,\quad\mbox{in}\quad\mathbb{R}^{2},\quad
V_{k}(e_{1})=\mu_{k}.$

This expression is based on the classification theorem of Prajapat-Tarantello
[23].

The estimate of $v_{k}(x)-V_{k}(x)$ is important for the main theorem of this
article. For convenience we use

$\beta_{l}=\frac{2\pi
l}{N+1},\quad\mbox{so}\,\,e_{1}=e^{i\beta_{0}}=Q_{0}^{k},\quad
e^{i\beta_{l}}=Q_{l}^{k}+E,\,\,\mbox{ for }\,\,l=1,...,N.$

##  4\. Vanishing of the first derivatives

Our first goal is to prove the following vanishing rate for
$\nabla\mathfrak{h}_{k}(0)$:

######  Theorem 4.1.

(4.1) $\nabla(\log\mathfrak{h}_{k})(0)=O(\delta_{k}\mu_{k})$

Proof of Theorem 4.1:

Note that we have proved in [26] that

$\nabla(\log\mathfrak{h}_{k})(0)=O(\delta_{k}^{-1}\mu_{k}e^{-\mu_{k}})+O(\delta_{k}).$

If $\delta_{k}\geq C\epsilon_{k}$, there is nothing to prove. So we assume
that

(4.2) $\delta_{k}=o(\epsilon_{k}).$

By way of contradiction we assume that

(4.3) $|\nabla\mathfrak{h}_{k}(0)|/(\delta_{k}\mu_{k})\to\infty.$

Another observation is that based on (3.5) we have

$\epsilon_{k}^{-1}|Q_{l}^{k}-e^{i\beta_{l}}|\leq
C\epsilon_{k}^{\epsilon},\quad l=0,...,N$

for some small $\epsilon>0$. Thus $\xi_{k}$ tends to $U$ after scaling. We
need this fact in our argument.

Under the assumption (4.2) we cite Proposition 3.1 of [27]:

_Proposition 3.1 of[27]: Let $l=0,...,N$ and $\delta$ be small so that
$B(e^{i\beta_{l}},\delta)\cap B(e^{i\beta_{s}},\delta)=\emptyset$ for $l\neq
s$. In each $B(e^{i\beta_{l}},\delta)$_

(4.4)
$|v_{k}(x)-V_{k}(x)|\leq\left\\{\begin{array}[]{ll}C\mu_{k}e^{-\mu_{k}/2},\quad|x-e^{i\beta_{l}}|\leq
Ce^{-\mu_{k}/2},\\\ \\\
C\frac{\mu_{k}e^{-\mu_{k}}}{|x-e^{i\beta_{l}}|}+O(\mu_{k}^{2}e^{-\mu_{k}}),\quad
Ce^{-\mu_{k}/2}\leq|x-e^{i\beta_{l}}|\leq\delta.\end{array}\right.$

######  Remark 4.1.

We only need a re-scaled version of the Proposition above:

(4.5)
$|v_{k}(e^{i\beta_{l}}+\epsilon_{k}y)-V_{k}(e^{i\beta_{l}}+\epsilon_{k}y)|\leq
C\epsilon_{k}^{\epsilon}(1+|y|)^{-1},\quad 0<|y|<\tau\epsilon_{k}^{-1}.$

for some small constants $\epsilon>0$ and $\tau>0$ both independent of $k$,

One major step in the proof of Theorem 4.1 is the following estimate:

######  Proposition 4.1.

Let $w_{k}=v_{k}-V_{k}$, then

$|w_{k}(y)|\leq C\tilde{\delta_{k}},\quad
y\in\Omega_{k}:=B(0,\tau\delta_{k}^{-1}),$

where
$\tilde{\delta_{k}}=|\nabla\mathfrak{h}_{k}(0)|\delta_{k}+\delta_{k}^{2}\mu_{k}$.

Proof of Proposition 4.1:

Obviously we can assume that
$|\nabla\mathfrak{h}_{k}(0)|\delta_{k}>2\delta_{k}^{2}\mu_{k}$ because
otherwise there is nothing to prove. Now we recall the equation for $v_{k}$ is
(3.2), $v_{k}$ is a constant on $\partial B(0,\tau\delta_{k}^{-1})$. Moreover
$v_{k}(e_{1})=\mu_{k}$. Recall that $V_{k}$ defined in (3.7) satisfies

$\Delta
V_{k}+\mathfrak{h}_{k}(\delta_{k}e_{1})|y|^{2N}e^{V_{k}}=0,\quad\mbox{in}\quad\mathbb{R}^{2},\quad\int_{\mathbb{R}^{2}}|y|^{2N}e^{V_{k}}<\infty,$

$V_{k}$ has its local maximums at $e^{i\beta_{l}}$ for $l=0,...,N$ and
$V_{k}(e_{1})=\mu_{k}$. For $|y|\sim\delta_{k}^{-1}$,

$V_{k}(y)=-\mu_{k}-4(N+1)\log\delta_{k}^{-1}+C+O(\delta_{k}^{N+1}).$

Let $\Omega_{k}=B(0,\tau\delta_{k}^{-1})$, we shall derive a precise, point-
wise estimate of $w_{k}$ in $B_{3}\setminus\cup_{l=1}^{N}B(Q_{l}^{k},\tau)$
where $\tau>0$ is a small number independent of $k$. Here we note that among
$N+1$ local maximum points, we already have $e_{1}$ as a common local maximum
point for both $v_{k}$ and $V_{k}$ and we shall prove that $w_{k}$ is very
small in $B_{3}$ if we exclude all bubbling disks except the one around
$e_{1}$. Before we carry out more specific computation we emphasize the
importance of

(4.6) $w_{k}(e_{1})=|\nabla w_{k}(e_{1})|=0.$

Now we write the equation of $w_{k}$ as

(4.7) $\Delta
w_{k}+\mathfrak{h}_{k}(\delta_{k}y)|y|^{2N}e^{\xi_{k}}w_{k}=(\mathfrak{h}_{k}(\delta_{k}e_{1})-\mathfrak{h}_{k}(\delta_{k}y))|y|^{2N}e^{V_{k}}$

in $\Omega_{k}$, where $\xi_{k}$ is obtained from the mean value theorem:

$e^{\xi_{k}(x)}=\left\\{\begin{array}[]{ll}\frac{e^{v_{k}(x)}-e^{V_{k}(x)}}{v_{k}(x)-V_{k}(x)},\quad\mbox{if}\quad
v_{k}(x)\neq V_{k}(x),\\\ \\\ e^{V_{k}(x)},\quad\mbox{if}\quad
v_{k}(x)=V_{k}(x).\end{array}\right.$

An equivalent form is

(4.8)
$e^{\xi_{k}(x)}=\int_{0}^{1}\frac{d}{dt}e^{tv_{k}(x)+(1-t)V_{k}(x)}dt=e^{V_{k}(x)}\big{(}1+\frac{1}{2}w_{k}(x)+O(w_{k}(x)^{2})\big{)}.$

For convenience we write the equation for $w_{k}$ as

(4.9) $\Delta
w_{k}+\mathfrak{h}_{k}(\delta_{k}y)|y|^{2N}e^{\xi_{k}}w_{k}=\delta_{k}\nabla\mathfrak{h}_{k}(\delta_{k}e_{1})\cdot(e_{1}-y)|y|^{2N}e^{V_{k}}+E_{1}$

where

$E_{1}=O(\delta_{k}^{2})|y-e_{1}|^{2}|y|^{2N}e^{V_{k}},\quad y\in\Omega_{k}.$

Note that the oscillation of $w_{k}$ on $\partial\Omega_{k}$ is
$O(\delta_{k}^{N+1})$, which all comes from the oscillation of $V_{k}$.

Let $M_{k}=\max_{x\in\bar{\Omega}_{k}}|w_{k}(x)|$. We shall get a
contradiction by assuming $M_{k}/\tilde{\delta_{k}}\to\infty$. This assumption
implies

(4.10) $M_{k}/(\delta_{k}^{2}\mu_{k})\to\infty.$

Set

$\tilde{w}_{k}(y)=w_{k}(y)/M_{k},\quad x\in\Omega_{k}.$

Clearly $\max_{x\in\Omega_{k}}|\tilde{w}_{k}(x)|=1$. The equation for
$\tilde{w}_{k}$ is

(4.11)
$\Delta\tilde{w}_{k}(y)+|y|^{2N}\mathfrak{h}_{k}(\delta_{k}e_{1})e^{\xi_{k}}\tilde{w}_{k}(y)=a_{k}\cdot(e_{1}-y)|y|^{2N}e^{V_{k}}+\tilde{E}_{1},$

in $\Omega_{k}$, where $a_{k}=\delta_{k}\nabla\mathfrak{h}_{k}(0)/M_{k}\to 0$,

(4.12) $\tilde{E}_{1}=o(1)|y-e_{1}|^{2}|y|^{2N}e^{V_{k}},\quad
y\in\Omega_{k}.$

Also on the boundary, since $M_{k}/\tilde{\delta}_{k}\to\infty$, we have

(4.13) $\tilde{w}_{k}=C+o(1/\mu_{k}),\quad\mbox{on}\quad\partial\Omega_{k}.$

By Proposition 3.1 of [27]

(4.14)
$\xi_{k}(e_{1}+\epsilon_{k}z)=V_{k}(e_{1}+\epsilon_{k}z)+O(\epsilon_{k}^{\epsilon})(1+|z|)^{-1}$

Since $V_{k}$ is not exactly symmetric around $e_{1}$, we shall replace the
re-scaled version of $V_{k}$ around $e_{1}$ by a radial function. Let $U_{k}$
be solutions of

(4.15) $\Delta
U_{k}+\mathfrak{h}_{k}(\delta_{k}e_{1})e^{U_{k}}=0,\quad\mbox{in}\quad\mathbb{R}^{2},\quad
U_{k}(0)=\max_{\mathbb{R}^{2}}U_{k}=0.$

By the classification theorem of Caffarelli-Gidas-Spruck [7] we have

$U_{k}(z)=\log\frac{1}{(1+\frac{\mathfrak{h}_{k}(\delta_{k}e_{1})}{8}|z|^{2})^{2}}$

and standard refined estimates yield (see [12, 30, 17])

(4.16)
$V_{k}(e_{1}+\epsilon_{k}z)+2\log\epsilon_{k}=U_{k}(z)+O(\epsilon_{k})|z|+O(\mu_{k}^{2}\epsilon_{k}^{2}).$

Also we observe that

(4.17) $\log|e_{1}+\epsilon_{k}z|=O(\epsilon_{k})|z|.$

Thus, the combination of (4.14), (4.16) and (4.17) gives

(4.18) $\displaystyle
2N\log|e_{1}+\epsilon_{k}z|+\xi_{k}(e_{1}+\epsilon_{k}z)+2\log\epsilon_{k}-U_{k}(z)$
$\displaystyle=$ $\displaystyle O(\epsilon_{k}^{\epsilon})(1+|z|)\quad
0\leq|z|<\delta_{0}\epsilon_{k}^{-1}.$

for a small $\epsilon>0$ independent of $k$. Since we shall use the re-scaled
version, based on (4.18) we have

(4.19)
$\epsilon_{k}^{2}|e_{1}+\epsilon_{k}z|^{2N}e^{\xi_{k}(e_{1}+\epsilon_{k}z)}=e^{U_{k}(z)}+O(\epsilon_{k}^{\epsilon})(1+|z|)^{-3}$

Here we note that the estimate in (4.18) is not optimal. In the following we
shall put the proof of Proposition 4.1 into a few estimates. In the first
estimate we prove

######  Lemma 4.1.

For $\delta>0$ small and independent of $k$,

(4.20) $\tilde{w}_{k}(y)=o(1),\quad\nabla\tilde{w}_{k}=o(1)\quad\mbox{in}\quad
B(e_{1},\delta)\setminus B(e_{1},\delta/8)$

where $B(e_{1},3\delta)$ does not include other blowup points.

Proof of Lemma 4.1:

If (4.20) is not true, we have, without loss of generality that
$\tilde{w}_{k}\to c>0$. This is based on the fact that $\tilde{w}_{k}$ tends
to a global harmonic function with removable singularity. So $\tilde{w}_{k}$
tends to constant. Here we assume $c>0$ but the argument for $c<0$ is the
same. Let

(4.21)
$W_{k}(z)=\tilde{w}_{k}(e_{1}+\epsilon_{k}z),\quad\epsilon_{k}=e^{-\frac{1}{2}\mu_{k}},$

then if we use $W$ to denote the limit of $W_{k}$, we have

$\Delta W+e^{U}W=0,\quad\mathbb{R}^{2},\quad|W|\leq 1,$

and $U$ is a solution of $\Delta U+e^{U}=0$ in $\mathbb{R}^{2}$ with
$\int_{\mathbb{R}^{2}}e^{U}<\infty$. Since $0$ is the local maximum of $U$,

$U(z)=\log\frac{1}{(1+\frac{1}{8}|z|^{2})^{2}}.$

Here we further claim that $W\equiv 0$ in $\mathbb{R}^{2}$ because
$W(0)=|\nabla W(0)|=0$, a fact well known based on the classification of the
kernel of the linearized operator. Going back to $W_{k}$, we have

$W_{k}(z)=o(1),\quad|z|\leq R_{k}\mbox{ for some }\quad R_{k}\to\infty.$

Based on the expression of $\tilde{w}_{k}$, (4.16) and (4.19) we write the
equation of $W_{k}$ as

(4.22) $\Delta
W_{k}(z)+\mathfrak{h}_{k}(\delta_{k}e_{1})e^{U_{k}(z)}W_{k}(z)=E_{2}^{k},$

for $|z|<\delta_{0}\epsilon_{k}^{-1}$ where a crude estimate of the error term
$E_{2}^{k}$ is

$E_{2}^{k}(z)=o(1)\epsilon_{k}^{\epsilon}(1+|z|)^{-3}.$

Let

(4.23) $g_{0}^{k}(r)=\frac{1}{2\pi}\int_{0}^{2\pi}W_{k}(r,\theta)d\theta.$

Then clearly $g_{0}^{k}(r)\to c>0$ for $r\sim\epsilon_{k}^{-1}$. The equation
for $g_{0}^{k}$ is

$\displaystyle\frac{d^{2}}{dr^{2}}g_{0}^{k}(r)+\frac{1}{r}\frac{d}{dr}g_{0}^{k}(r)+\mathfrak{h}_{k}(\delta_{k}e_{1})e^{U_{k}(r)}g_{0}^{k}(r)=\tilde{E}_{0}^{k}(r)$
$\displaystyle g_{0}^{k}(0)=\frac{d}{dr}g_{0}^{k}(0)=0.$

where $\tilde{E}_{0}^{k}(r)$ has the same upper bound as that of
$E_{2}^{k}(r)$:

$|\tilde{E}_{0}^{k}(r)|\leq o(1)\epsilon_{k}^{\epsilon}(1+r)^{-3}.$

For the homogeneous equation, the two fundamental solutions are known:
$g_{01}$, $g_{02}$, where

$g_{01}=\frac{1-c_{1}r^{2}}{1+c_{1}r^{2}},\quad
c_{1}=\frac{\mathfrak{h}_{k}(\delta_{k}e_{1})}{8}.$

By the standard reduction of order process, $g_{02}(r)=O(\log r)$ for $r>1$.
Then it is easy to obtain, assuming $|W_{k}(z)|\leq 1$, that

$\displaystyle|g_{0}(r)|\leq
C|g_{01}(r)|\int_{0}^{r}s|\tilde{E}_{0}^{k}(s)g_{02}(s)|ds+C|g_{02}(r)|\int_{0}^{r}s|g_{01}(s)\tilde{E}_{0}^{k}(s)|ds$
$\displaystyle\leq C\epsilon_{k}^{\epsilon}\log(2+r).\quad
0<r<\delta_{0}\epsilon_{k}^{-1}.$

Clearly this is a contradiction to (4.23). We have proved $c=0$, which means
$\tilde{w}_{k}=o(1)$ in $B(e_{1},\delta_{0})\setminus B(e_{1},\delta_{0}/8)$.
Then it is easy to use the equation for $\tilde{w}_{k}$ and standard Harnack
inequality to prove $\nabla\tilde{w}_{k}=o(1)$ in the same region. Lemma 4.1
is established. $\Box$

The second estimate is a more precise description of $\tilde{w}_{k}$ around
$e_{1}$:

######  Lemma 4.2.

For any given $\sigma\in(0,1)$ there exists $C>0$ such that

(4.24) $|\tilde{w}_{k}(e_{1}+\epsilon_{k}z)|\leq
C\epsilon_{k}^{\sigma}(1+|z|)^{\sigma},\quad 0<|z|<\tau\epsilon_{k}^{-1}.$

for some $\tau>0$.

Proof of Lemma 4.2: Let $W_{k}$ be defined as in (4.21). In order to obtain a
better estimate we need to write the equation of $W_{k}$ more precisely than
(4.22):

(4.25) $\Delta
W_{k}+\mathfrak{h}_{k}(\delta_{k}e_{1})e^{\Theta_{k}}W_{k}=E_{3}^{k}(z),\quad
z\in\Omega_{Wk}$

where $\Theta_{k}$ is defined by

$e^{\Theta_{k}(z)}=|e_{1}+\epsilon_{k}z|^{2N}e^{\xi_{k}(e_{1}+\epsilon_{k}z)+2\log\epsilon_{k}},$

$\Omega_{Wk}=B(0,\tau\epsilon_{k}^{-1})$ and $E_{3}^{k}(z)$ satisfies

$E_{3}^{k}(z)=O(\epsilon_{k})(1+|z|)^{-3},\quad z\in\Omega_{Wk}.$

Here we observe that by Lemma 4.1 $W_{k}=o(1)$ on $\partial\Omega_{Wk}$. Let

$\Lambda_{k}=\max_{z\in\Omega_{Wk}}\frac{|W_{k}(z)|}{\epsilon_{k}^{\sigma}(1+|z|)^{\sigma}}.$

If (4.24) does not hold, $\Lambda_{k}\to\infty$ and we use $z_{k}$ to denote
where $\Lambda_{k}$ is attained. Note that because of the smallness of $W_{k}$
on $\partial\Omega_{Wk}$, $z_{k}$ is an interior point. Let

$g_{k}(z)=\frac{W_{k}(z)}{\Lambda_{k}(1+|z_{k}|)^{\sigma}\epsilon_{k}^{\sigma}},\quad
z\in\Omega_{Wk},$

we see immediately that

(4.26)
$|g_{k}(z)|=\frac{|W_{k}(z)|}{\epsilon_{k}^{\sigma}\Lambda_{k}(1+|z|)^{\sigma}}\cdot\frac{(1+|z|)^{\sigma}}{(1+|z_{k}|)^{\sigma}}\leq\frac{(1+|z|)^{\sigma}}{(1+|z_{k}|)^{\sigma}}.$

Note that $\sigma$ can be as close to $1$ as needed. The equation of $g_{k}$
is

$\Delta
g_{k}(z)+\mathfrak{h}_{k}(\delta_{k}e_{1})e^{\Theta_{k}}g_{k}=o(\epsilon_{k}^{1-\sigma})\frac{(1+|z|)^{-3}}{(1+|z_{k}|)^{\sigma}},\quad\mbox{in}\quad\Omega_{Wk}.$

Then we can obtain a contradiction to $|g_{k}(z_{k})|=1$ as follows: If
$\lim_{k\to\infty}z_{k}=P\in\mathbb{R}^{2}$, this is not possible because that
fact that $g_{k}(0)=|\nabla g_{k}(0)|=0$ and the sub-linear growth of $g_{k}$
in (4.26) implies that $g_{k}\to 0$ over any compact subset of
$\mathbb{R}^{2}$ (see [12, 30]). So we have $|z_{k}|\to\infty$. But this would
lead to a contradiction again by using the Green’s representation of $g_{k}$:

(4.27) $\displaystyle\pm 1=g_{k}(z_{k})=g_{k}(z_{k})-g_{k}(0)$
$\displaystyle=\int_{\Omega_{k,1}}(G_{k}(z_{k},\eta)-G_{k}(0,\eta))(\mathfrak{h}_{k}(\delta_{k}e_{1})e^{\Theta_{k}}g_{k}(\eta)+o(\epsilon_{k}^{1-\sigma})\frac{(1+|\eta|)^{-3}}{(1+|z_{k}|)^{\sigma}})d\eta+o(1).$

where $G_{k}(y,\eta)$ is the Green’s function on $\Omega_{Wk}$ and $o(1)$ in
the equation above comes from the smallness of $W_{k}$ on
$\partial\Omega_{Wk}$. Let $L_{k}=\tau\epsilon_{k}^{-1}$, the expression of
$G_{k}$ is

$G_{k}(y,\eta)=-\frac{1}{2\pi}\log|y-\eta|+\frac{1}{2\pi}\log(\frac{|\eta|}{L_{k}}|\frac{L_{k}^{2}\eta}{|\eta|^{2}}-y|).$
$G_{k}(z_{k},\eta)-G_{k}(0,\eta)=-\frac{1}{2\pi}\log|z_{k}-\eta|+\frac{1}{2\pi}\log|\frac{z_{k}}{|z_{k}|}-\frac{\eta
z_{k}}{L_{k}^{2}}|+\frac{1}{2\pi}\log|\eta|.$

Using this expression in (4.27) we obtain from elementary computation that the
right hand side of (4.27) is $o(1)$, a contradiction to $|g_{k}(z_{k})|=1$.
Lemma 4.2 is established. $\Box$

The smallness of $\tilde{w}_{k}$ around $e_{1}$ can be used to obtain the
following third key estimate:

######  Lemma 4.3.

(4.28) $\tilde{w}_{k}=o(1)\quad\mbox{in}\quad B(e^{i\beta_{l}},\tau)\quad
l=1,..,N.$

Proof of Lemma 4.3: We abuse the notation $W_{k}$ by defining it as

$W_{k}(z)=\tilde{w}_{k}(e^{i\beta_{l}}+\epsilon_{k}z),\quad
z\in\Omega_{k,l}:=B(0,\tau\epsilon_{k}^{-1}).$

Here we point out that based on (3.5) and (4.2) we have
$\epsilon_{k}^{-1}|Q_{l}^{k}-e^{i\beta_{l}}|\to 0$. So the scaling around
$e^{i\beta_{l}}$ or $Q_{l}^{k}$ does not affect the limit function.

$\epsilon_{k}^{2}|e^{i\beta_{l}}+\epsilon_{k}z|^{2N}\mathfrak{h}_{k}(\delta_{k}e_{1})e^{\xi_{k}(e^{i\beta_{l}}+\epsilon_{k}z)}\to
e^{U(z)}$

where $U(z)$ is a solution of

$\Delta
U+e^{U}=0,\quad\mbox{in}\quad\mathbb{R}^{2},\quad\int_{\mathbb{R}^{2}}e^{U}<\infty.$

Here we recall that $\lim_{k\to\infty}\mathfrak{h}_{k}(\delta_{k}e_{1})=1$.
Since $W_{k}$ converges to a solution of the linearized equation:

$\Delta W+e^{U}W=0,\quad\mbox{in}\quad\mathbb{R}^{2}.$

$W$ can be written as a linear combination of three functions:

$W(x)=c_{0}\phi_{0}+c_{1}\phi_{1}+c_{2}\phi_{2},$

where

$\phi_{0}=\frac{1-\frac{1}{8}|x|^{2}}{1+\frac{1}{8}|x|^{2}}$
$\phi_{1}=\frac{x_{1}}{1+\frac{1}{8}|x|^{2}},\quad\phi_{2}=\frac{x_{2}}{1+\frac{1}{8}|x|^{2}}.$

The remaining part of the proof consisting of proving $c_{0}=0$ and
$c_{1}=c_{2}=0$. First we prove $c_{0}=0$.

Step one: $c_{0}=0$. First we write the equation for $W_{k}$ in a convenient
form. Since

$|e^{i\beta_{l}}+\epsilon_{k}z|^{2N}\mathfrak{h}_{k}(\delta_{k}e_{1})=\mathfrak{h}_{k}(\delta_{k}e_{1})+O(\epsilon_{k}z),$

and

$\epsilon_{k}^{2}e^{\xi_{k}(e^{i\beta_{l}}+\epsilon_{k}z)}=e^{U_{k}(z)}+O(\epsilon_{k}^{\epsilon})(1+|z|)^{-3}.$

Based on (4.11) we write the equation for $W_{k}$ as

(4.29) $\Delta
W_{k}(z)+\mathfrak{h}_{k}(\delta_{k}e_{1})e^{U_{k}}W_{k}=E_{l}^{k}(z)$

where

$E_{l}^{k}(z)=O(\epsilon_{k}^{\epsilon})(1+|z|)^{-3}\quad\mbox{in}\quad\Omega_{k,l}.$

In order to prove $c_{0}=0$, the key is to control the derivative of
$W_{0}^{k}(r)$ where

$W_{0}^{k}(r)=\frac{1}{2\pi r}\int_{\partial B_{r}}W_{k}(re^{i\theta})dS,\quad
0<r<\tau\epsilon_{k}^{-1}.$

To obtain a control of $\frac{d}{dr}W_{0}^{k}(r)$ we use $\phi_{0}^{k}(r)$ as
the radial solution of

$\Delta\phi_{0}^{k}+\mathfrak{h}_{k}(\delta_{k}e_{1})e^{U_{k}}\phi_{0}^{k}=0,\quad\mbox{in
}\quad\mathbb{R}^{2}.$

When $k\to\infty$, $\phi_{0}^{k}\to c_{0}\phi_{0}$. Thus using the equation
for $\phi_{0}^{k}$ and $W_{k}$, we have

(4.30) $\int_{\partial
B_{r}}(\partial_{\nu}W_{k}\phi_{0}^{k}-\partial_{\nu}\phi_{0}^{k}W_{k})=o(\epsilon_{k}^{\epsilon}).$

Thus from (4.30) we have

(4.31) $\frac{d}{dr}W_{0}^{k}(r)=\frac{1}{2\pi r}\int_{\partial
B_{r}}\partial_{\nu}W_{k}=o(\epsilon_{k}^{\epsilon})/r+O(1/r^{3}),\quad
1<r<\tau\epsilon_{k}^{-1}.$

Since we have known that

$W_{0}^{k}(\tau\epsilon_{k}^{-1})=o(1).$

By the fundamental theorem of calculus we have

$W_{0}^{k}(r)=W_{0}^{k}(\tau\epsilon_{k}^{-1})+\int_{\tau\epsilon_{k}^{-1}}^{r}(\frac{o(\epsilon_{k}^{\epsilon})}{s}+O(s^{-3}))ds=O(1/r^{2})+O(\epsilon_{k}^{\epsilon}\log\frac{1}{\epsilon_{k}})$

for $r\geq 1$. Thus $c_{0}=0$ because $W_{0}^{k}(r)\to c_{0}\phi_{0}$, which
means when $r$ is large, it is $-c_{0}+O(1/r^{2})$.

Step two $c_{1}=c_{2}=0$. We first observe that Lemma 4.3 follows from this.
Indeed, once we have proved $c_{1}=c_{2}=c_{0}=0$ around each
$e^{i\beta_{l}}$, it is easy to use maximum principle to prove
$\tilde{w}_{k}=o(1)$ in $B_{3}$ using $\tilde{w}_{k}=o(1)$ on $\partial B_{3}$
and the Green’s representation of $\tilde{w}_{k}$. The smallness of
$\tilde{w}_{k}$ immediately implies $\tilde{w}_{k}=o(1)$ in $B_{R}$ for any
fixed $R>>1$. Outside $B_{R}$, a crude estimate of $v_{k}$ is

$v_{k}(y)\leq-\mu_{k}-4(N+1)\log|y|+C,\quad 3<|y|<\tau\delta_{k}^{-1}.$

Using this and the Green’s representation of $w_{k}$ we can first observe that
the oscillation on each $\partial B_{r}$ is $o(1)$
($R<r<\tau\delta_{k}^{-1}/2$) and then by the Green’s representation of
$\tilde{w}_{k}$ and fast decay rate of $e^{V_{k}}$ we obtain
$\tilde{w}_{k}=o(1)$ in $\overline{B(0,\tau\delta_{k}^{-1})}$. A contradiction
to $\max|\tilde{w}_{k}|=1$.

There are $N+1$ local maximums with one of them being $e_{1}$. Correspondingly
there are $N+1$ global solutions $V_{l,k}$ that approximate $v_{k}$ accurately
near $Q_{l}^{k}$ for $l=0,...,N$. Note that $Q_{0}^{k}=e_{1}$. For $V_{l,k}$
the expression is

$V_{l,k}=\log\frac{e^{\mu_{l}^{k}}}{(1+\frac{e^{\mu_{l}^{k}}}{D_{l}^{k}}|y^{N+1}-(e_{1}+p_{l}^{k})|^{2})^{2}},\quad
l=0,...,N,$

where $p_{l}^{k}=E$ and

(4.32) $D_{l}^{k}=8(N+1)^{2}/\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k}).$

The equation that $V_{l,k}$ satisfies is

$\Delta
V_{l,k}+|y|^{2N}\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})e^{V_{l,k}}=0,\quad\mbox{in}\quad\mathbb{R}^{2}.$

Since $v_{k}$ and $V_{l,k}$ have the same common local maximum at $Q_{l}^{k}$,
it is easy to see that

(4.33)
$Q_{l}^{k}=e^{i\beta_{l}}+\frac{p_{l}^{k}e^{i\beta_{l}}}{N+1}+O(|p_{l}^{k}|^{2}),\quad\beta_{l}=\frac{2l\pi}{N+2}.$

Let $M_{l,k}$ be the maximum of $|v_{k}-V_{l,k}|$ and we claim that all these
$M_{l,k}$ are comparable:

(4.34) $M_{l,k}\sim M_{s,k},\quad\forall s\neq l.$

The proof of (4.34) is as follows: We use $L_{s,l}$ to denote the limit of
$(v_{k}-V_{l,k})/M_{l,k}$ around $Q_{s}^{k}$:

$\frac{(v_{k}-V_{l,k})(Q_{s}^{k}+\epsilon_{k}z)}{M_{l,k}}=L_{s,l}+o(1),\quad|z|\leq\tau\epsilon_{k}^{-1}$

where

$L_{s,l}=c_{1,s,l}\frac{z_{1}}{1+\frac{1}{8}|z|^{2}}+c_{2,s,l}\frac{z_{2}}{1+\frac{1}{8}|z|^{2}},\quad\mbox{and}\quad
L_{l,l}=0,\quad s=0,...,N.$

If all $c_{1,s,l}$ and $c_{2,s,l}$ are zero for a fixed $l$, we can obtain a
contradiction just like the beginning of step two. So at least one of them is
not zero. For each $s\neq l$, by Lemma 4.2 we have

(4.35)
$v_{k}(Q_{s}^{k}+\epsilon_{k}z)-V_{s,k}(Q_{s}^{k}+\epsilon_{k}z)=O(\epsilon_{k}^{\sigma})(1+|z|)^{\sigma}M_{s,k},\quad|z|<\tau\epsilon_{k}^{-1}.$

Let $M_{k}=\max_{i}M_{i,k}$ ($i=0,...,N$) and we suppose $M_{k}=M_{l,k}$. Then
to determine $L_{s,l}$ we see that

$\displaystyle\frac{v_{k}(Q_{s}^{k}+\epsilon_{k}z)-V_{l,k}(Q_{s}^{k}+\epsilon_{k}z)}{M_{k}}$
$\displaystyle=$ $\displaystyle
o(\epsilon_{k}^{\sigma})(1+|z|)^{\sigma}+\frac{V_{s,k}(Q_{s}^{k}+\epsilon_{k}z)-V_{l,k}(Q_{s}^{k}+\epsilon_{k}z)}{M_{k}}.$

This expression says that $L_{s,l}$ is mainly determined by the difference of
two global solutions $V_{s,k}$ and $V_{l,k}$. In order to obtain a
contradiction to our assumption we will put the difference in several terms.
The main idea in this part of the reasoning is that “first order terms” tell
us what the kernel functions should be, then the “second order terms” tell us
where the pathology is.

We write $V_{s,k}(y)-V_{l,k}(y)$ as

$V_{s,k}(y)-V_{l,k}(y)=\mu_{s}^{k}-\mu_{l}^{k}+2A-A^{2}+O(|A|^{3})$

where

$A(y)=\frac{\frac{e^{\mu_{l}^{k}}}{D_{l}^{k}}|y^{N+1}-e_{1}-p_{l}^{k}|^{2}-\frac{e^{\mu_{s}^{k}}}{D_{s}^{k}}|y^{N+1}-e_{1}-p_{s}^{k}|^{2}}{1+\frac{e^{\mu_{s}^{k}}}{D_{s}^{k}}|y^{N+1}-e_{1}-p_{s}^{k}|^{2}}.$

Here for convenience we abuse the notation $\epsilon_{k}$ by assuming
$\epsilon_{k}=e^{-\mu_{s}^{k}/2}$. Note that $\epsilon_{k}=e^{-\mu_{t}^{k}/2}$
for some $t$, but it does not matter which $t$ it is. From $A$ we claim that

(4.36) $\displaystyle
V_{s,k}(Q_{s}^{k}+\epsilon_{k}z)-V_{l,k}(Q_{s}^{k}+\epsilon_{k}z)$
$\displaystyle=$
$\displaystyle\phi_{1}+\phi_{2}+\phi_{3}+\phi_{4}+\mathfrak{R},$

where

$\displaystyle\phi_{1}=(\mu_{s}^{k}-\mu_{l}^{k})(1-\frac{(N+1)^{2}}{D_{s}^{k}}|z+O(\epsilon_{k})|z|^{2}|^{2})/B,$
$\displaystyle\phi_{2}=\frac{2(N+1)^{2}}{D_{s}^{k}}\delta_{k}\nabla\mathfrak{h}_{k}(\delta_{k}Q_{s}^{k})(Q_{l}^{k}-Q_{s}^{k})|z|^{2}/B$
$\displaystyle\phi_{3}=\frac{4(N+1)}{D_{s}^{k}B}Re((z+O(\epsilon_{k}|z|^{2}))(\frac{\bar{p}_{s}^{k}-\bar{p}_{l}^{k}}{\epsilon_{k}}e^{-i\beta_{s}}))$
$\displaystyle\phi_{4}=\frac{|p_{s}^{k}-p_{l}^{k}|^{2}}{\epsilon_{k}^{2}}\bigg{(}\frac{2}{D_{s}^{k}B}-\frac{2(N+1)^{2}|z|^{2}}{D_{s}^{2}B^{2}}-\frac{2(N+1)^{2}}{D_{s}^{2}B^{2}}|z|^{2}\cos(2\theta-2\theta_{st}-2\beta_{s})\bigg{)},$
$\displaystyle B=1+\frac{(N+1)^{2}}{D_{s}^{k}}|z+O(\epsilon_{k}|z|^{2})|^{2},$

and $\mathfrak{R}_{k}$ is the collections of other insignificant terms. Here
we briefly explain the roles of each term. $\phi_{1}$ corresponds to the
radial solution in the kernel of the linearized operator of the global
equation. In other words, $\phi_{1}^{k}/M_{k}$ should tend to zero because in
step one we have proved $c_{0}=0$. $\phi_{2}^{k}/M_{k}$ is the combination of
the two other functions in the kernel. $\phi_{4}$ is the second order term
which will play a leading role later. $\phi_{3}^{k}$ comes from the difference
of $\mathfrak{h}_{k}$ at $Q_{l}^{k}$ and $Q_{s}^{k}$. The derivation of (4.36)
is as follows: First by the expression of $Q_{s}^{k}$ in (4.33) we have

$y^{N+1}=1+p_{s}^{k}+(N+1)\epsilon_{k}ze^{-i\beta_{s}}+O(\epsilon_{k}^{2})|z|^{2},$

where $y=Q_{s}^{k}+\epsilon_{k}z$. Then

$|y^{N+1}-e_{1}-p_{s}^{k}|^{2}=(N+1)^{2}\epsilon_{k}^{2}|z+O(\epsilon_{k})z^{2}|^{2}+O(\epsilon_{k}^{3})|z|^{3}$
$|y^{N+1}-e_{1}-p_{l}^{k}|^{2}=(N+1)^{2}\epsilon_{k}^{2}|z+\frac{(p_{s}^{k}-p_{l}^{k})e^{i\beta_{s}}}{(N+1)\epsilon_{k}}+O(\epsilon_{k})|z|^{2}|^{2}+O(\epsilon_{k}^{3}|z|^{3}).$

Next by the definition of $D_{s}^{k}$ in (4.32)

$\frac{D_{s}^{k}-D_{l}^{k}}{D_{l}^{k}}=\delta_{k}\nabla(\log\mathfrak{h}_{k})(0)\cdot(Q_{l}^{k}-Q_{s}^{k})+O(\delta_{k}^{2}).$
$\displaystyle\frac{e^{\mu_{l}^{k}-\mu_{s}^{k}}}{D_{l}^{k}}$
$\displaystyle=\frac{1}{D_{s}^{k}}(1+\frac{D_{s}^{k}-D_{l}^{k}}{D_{l}^{k}}+\mu_{l}^{k}-\mu_{s}^{k}+O(\mu_{l}^{k}-\mu_{s}^{k})^{2}+O(\delta_{k}^{2})).$
(4.37)
$\displaystyle=\frac{1}{D_{s}^{k}}(1+\delta_{k}\nabla\log\mathfrak{h}_{k}(0)\cdot(Q_{l}^{k}-Q_{s}^{k})+\mu_{l}^{k}-\mu_{s}^{k}+O(\mu_{l}^{k}-\mu_{s}^{k})^{2}+O(\delta_{k}^{2})).$

Then the expression of $A$ is ( for simplicity we omit $k$ in some notations)

$\displaystyle
A=\bigg{(}\frac{e^{\mu_{l}-\mu_{s}}}{D_{l}}(N+1)^{2}\big{(}|z|^{2}+2Re(z\frac{\overline{p_{s}-p_{l}}}{\epsilon_{k}(N+1)}e^{-i\beta_{s}})+\frac{|p_{s}-p_{l}|^{2}}{(N+1)^{2}\epsilon_{k}^{2}}+O(\epsilon_{k}|z|^{3})\big{)}$
$\displaystyle-\frac{(N+1)^{2}}{D_{s}}(|z|^{2}+O(\epsilon_{k})|z|^{3})\bigg{)}/B$

After using (4) we have

(4.38) $\displaystyle
A=\bigg{(}\frac{1}{D_{s}}(\delta_{k}\nabla(\log\mathfrak{h}_{k}(0)(Q_{l}-Q_{s})+\mu_{l}-\mu_{s}+O(\mu_{l}-\mu_{s})^{2})(N+1)^{2}|z|^{2}$
$\displaystyle+2Re(z\frac{\bar{p}_{s}-\bar{p}_{l}}{\epsilon_{k}}e^{-i\beta_{s}})(N+1)\frac{1}{D_{s}}+\frac{|p_{s}-p_{l}|^{2}}{\epsilon_{k}^{2}D_{s}}+O(\epsilon_{k})|z|^{3}+O(\delta_{k}^{2})|z|^{2}\bigg{)}/B.$

(4.39)
$A^{2}=(\frac{(N+1)^{2}}{D_{s}^{2}}4(Re(z\frac{\bar{p}_{s}-\bar{p}_{l}}{\epsilon_{k}}e^{-i\beta_{s}})^{2})/B^{2}+\mbox{other
terms}.$

The numerator of $A^{2}$ has the following leading term:

$\frac{(N+1)^{2}}{D_{s}^{2}}\bigg{(}2|z|^{2}(\frac{|p_{s}-p_{l}|}{\epsilon_{k}})^{2}\big{(}1+2\cos(2\theta-2\theta_{sl})\big{)}\bigg{)}$

where $z=|z|e^{i\theta}$, $p_{s}-p_{l}=|p_{s}-p_{l}|e^{i\theta_{sl}}$. Using
these expressions we can obtain (4.36) by direct computation. Here $\phi_{1}$,
$\phi_{3}$ correspond to solutions to the linearized operator. Here we note
that if we set $\epsilon_{l,k}=e^{-\mu_{l}^{k}/2}$, there is no essential
difference between $\epsilon_{l,k}$ and
$\epsilon_{k}=e^{-\frac{1}{2}\mu_{1,k}}$ because
$\epsilon_{l,k}=\epsilon_{k}+O(\epsilon_{k}E)$. If
$|\mu_{s,k}-\mu_{l,k}|/M_{k}\geq C$ there is no way to obtain a limit in the
form of $L_{s,l}$ mentioned before. Thus we must have
$|\mu_{s,k}-\mu_{l,k}|/M_{k}\to 0$. After simplification (see $\phi_{3}$ of
(4.36)) we have

(4.40) $\displaystyle
c_{1,s,l}=\lim_{k\to\infty}\frac{|p_{s}^{k}-p_{l}^{k}|}{2(N+1)M_{k}\epsilon_{k}}\cos(\beta_{s}+\theta_{sl}),$
$\displaystyle
c_{2,s,l}=\lim_{k\to\infty}\frac{|p_{s}^{k}-p_{l}^{k}|}{2(N+1)\epsilon_{k}M_{k}}sin(\beta_{s}+\theta_{sl})$

We omit $k$ for convenience. It is also important to observe that even if
$M_{k}=o(\epsilon_{k})$ we still have
$M_{k}\sim\max_{s}|p_{s}^{k}-p_{l}^{k}|/\epsilon_{k}$. Since each
$|p_{l}^{k}|=E$, an upper bound for $M_{k}$ is

(4.41) $M_{k}\leq C\mu_{k}\epsilon_{k}+C\delta_{k}^{2}\epsilon_{k}^{-1}.$

Equation (4.40) gives us a key observation:
$|c_{1,s,l}|+|c_{2,s,l}|\sim|p_{s}^{k}-p_{l}^{k}|/(\epsilon_{k}M_{k})$. So
whenever $|c_{1,s,l}|+|c_{2,s,l}|\neq 0$ we have
$\frac{|p_{s}^{k}-p_{l}^{k}|}{\epsilon_{k}}\sim M_{k}$. In other words for
each $l$, $M_{l,k}\sim\max_{t\neq
l}\frac{|p_{t}^{k}-p_{l}^{k}|}{\epsilon_{k}}$. Hence for any $t$, if
$\frac{|p_{t}^{k}-p_{l}^{k}|}{\epsilon_{k}}\sim M_{k}$, let $M_{t,k}$ be the
maximum of $|v_{k}-V_{t,k}|$, we have $M_{t,k}\sim M_{k}$. If all
$\frac{|p_{t}^{k}-p_{l}^{k}|}{\epsilon_{k}}\sim M_{k}$ (4.34) is proved. So we
prove that even if some $p_{t}^{k}$ is very close to $p_{l}^{k}$, $M_{t}^{k}$
is still comparable to $M_{k}$. The reason is there exists $q$ such that
$\frac{|p_{l}^{k}-p_{q}^{k}|}{\epsilon_{k}}\sim M_{k}$, if
$\frac{|p_{t}^{k}-p_{l}^{k}|}{\epsilon_{k}}=o(1)M_{k}$,

$|p_{t}^{k}-p_{q}^{k}|\geq|p_{l}^{k}-p_{q}^{k}|-|p_{t}^{k}-p_{l}^{k}|\geq\frac{1}{2}|p_{l}^{k}-p_{q}^{k}|.$

Thus $\frac{|p_{t}^{k}-p_{q}^{k}|}{\epsilon_{k}}\sim M_{k}$ and $M_{t}^{k}\sim
M_{k}$. (4.34) is established. From now on for convenience we shall just use
$M_{k}$. Since $M_{k}\sim\max_{s,t}|p_{s}^{k}-p_{t}^{k}|/\epsilon_{k}$, an
upper bound of $M_{k}$ is

(4.42) $M_{k}\leq C\mu_{k}\epsilon_{k}.$

Set $w_{l,k}=(v_{k}-V_{l,k})$, then we have $w_{l,k}(Q_{l}^{k})=|\nabla
w_{l,k}(Q_{l}^{k})|=0$. Correspondingly we set

$\tilde{w}_{l,k}=w_{l,k}/M_{k}.$

The equation of $w_{l,k}$ can be written as

(4.43) $\displaystyle\Delta
w_{l,k}+|y|^{2N}\mathfrak{h}_{k}(\delta_{k}Q_{l})e^{\xi_{l}}w_{l,k}$
$\displaystyle=$
$\displaystyle-\delta_{k}\nabla\mathfrak{h}_{k}(\delta_{k}Q_{l})(y-Q_{l})|y|^{2N}e^{V_{l,k}}-\delta_{k}^{2}\sum_{|\alpha|=2}\frac{\partial^{\alpha}\mathfrak{h}_{k}(\delta_{k}Q_{l})}{\alpha!}(y-Q_{l})^{\alpha}|y|^{2N}e^{V_{k}}$
$\displaystyle+O(\delta_{k}^{3})|y-Q_{l}|^{3}|y|^{2N}e^{V_{k}}$

where we omitted $k$ in $Q_{l}$ and $\xi_{l}$. $\xi_{l}$ comes from the Mean
Value Theorem and satisfies

(4.44) $e^{\xi_{l}}=e^{V_{l,k}}(1+\frac{1}{2}w_{l,k}+O(w_{l,k}^{2})).$

The function $\tilde{w}_{l,k}$ satisfies

(4.45)
$\lim_{k\to\infty}\tilde{w}_{l,k}(Q_{s}^{k}+\epsilon_{k}z)=\frac{c_{1,s,l}z_{1}+c_{2,s,l}z_{2}}{1+\frac{1}{8}|z|^{2}}$

and around each $Q_{s}^{k}$ (4.35) holds with $M_{s,k}$ replaced by $M_{k}$.

Now for $|y|\sim 1$, we use $w_{l,k}(Q_{l}^{k})=0$ to write $w_{l,k}(y)$ as

$\displaystyle w_{l,k}(y)$
$\displaystyle=\int_{\Omega_{k}}(G_{k}(y,\eta)-G_{k}(Q_{l},\eta))\bigg{(}\mathfrak{h}_{k}(\delta_{k}Q_{l})|\eta|^{2N}e^{\xi_{l}}w_{l,k}(\eta)$
$\displaystyle+\delta_{k}\nabla\mathfrak{h}_{k}(\delta_{k}Q_{l})(\eta-
Q_{l})|\eta|^{2N}e^{V_{l,k}}$
$\displaystyle+\delta_{k}^{2}\sum_{|\alpha|=2}\frac{\partial^{\alpha}\mathfrak{h}_{k}(\delta_{k}Q_{l})}{\alpha!}(\eta-
Q_{l})^{\alpha}|y|^{2N}e^{V_{k}}\bigg{)}+O(\delta_{k}^{N+2}).$

Note that the last term is $O(\delta_{k}^{N+2})$ because it comes from the
oscillation of $w_{l,k}$ on $\partial\Omega_{k}$. The harmonic function
defined by the boundary value of $w_{l,k}$ has an oscillation of
$O(\delta_{k}^{N+1})$ on $\partial\Omega_{k}$. The oscillation of this
harmonic function in $B_{R}$ (for any fixed $R>1$) is $O(\delta_{k}^{N+2})$.
The regular part of the Green’s function brings little error in the
computation, so we have

$\displaystyle\tilde{w}_{l,k}(y)$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi}\int_{\Omega_{k}}\log\frac{|y-\eta|}{|Q_{l}^{k}-\eta|}\bigg{(}\tilde{w}_{l,k}(\eta)\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})|\eta|^{2N}e^{\xi_{l}}$
$\displaystyle+\sigma_{k}\nabla\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})(\eta-
Q_{l}^{k})|\eta|^{2N}e^{V_{l,k}}+\frac{\delta_{k}^{2}}{M_{k}}\sum_{|\alpha|=2}\frac{\partial^{\alpha}\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})}{\alpha!}(\eta-
Q_{l}^{k})^{\alpha}|\eta|^{2N}e^{V_{l,k}}\bigg{)}d\eta$
$\displaystyle+o(\epsilon_{k})$

for $|y|\sim 1$.

Around each $Q_{s}^{k}$ the $e^{\xi_{l}}$ can be replaced by $e^{V_{s,k}}$
with controllable error (based on Lemma 4.2 and (4.42)). In order to evaluate
the expression of $\tilde{w}_{l,k}$ we need the following identity based on
(4.36)

(4.46)
$\int_{B(Q_{s}^{k},\tau)}(\tilde{w}_{l,k}(\eta)\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})|\eta|^{2N}e^{V_{s,k}}+\sigma_{k}\nabla\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})(\eta-
Q_{l})|\eta|^{2N}e^{V_{l,k}})d\eta=O(\epsilon_{k}^{\sigma})$

Note that $e^{V_{l,k}}$ in the first term was replaced by $e^{V_{s,k}}$ but in
the second term above this replacement is not necessary. (4.36) is mainly used
in the evaluation of the first term. The proof of (4.46) can be found in the
appendix.

Equation (4.46) also leads to a more accurate estimate of $\tilde{w}_{l,k}$ in
regions between bubbling disks. By the Green’s representation formula of
$\tilde{w}_{l,k}$ it is easy to have, for $|y|\sim 1$,

Writing the logarithmic term above as

$\log\frac{|y-\eta|}{|Q_{l}^{k}-\eta|}=\log\frac{|y-Q_{s}^{k}|}{|Q_{l}^{k}-Q_{s}^{k}|}+(\log\frac{|y-\eta|}{|Q_{l}^{k}-\eta|}-\log\frac{|y-Q_{s}^{k}|}{|Q_{l}^{k}-Q_{s}^{k}|}),$

then we see that the integration related to the second term is
$O(\epsilon_{k})$. The integration involving the first term is
$O(\epsilon_{k}^{\sigma})+o(\frac{\delta_{k}^{2}}{M_{k}})$ by (4.46) and the
definition of $M_{k}$. Therefore

$|\tilde{w}_{l,k}(y)|=o(1/\mu_{k}),\quad y\in
B_{3}\setminus\cup_{s=0}^{N}B(Q_{s}^{k},\tau)$

for $\sigma\in(0,1)$. Thus this extra control of $\tilde{w}_{l,k}$ away from
bubbling disks gives a better estimate than (4.35) around $Q_{l}^{k}$: Using
the same argument for Lemma 4.2 we have

(4.47) $|\tilde{w}_{l,k}(Q_{l}^{k}+\epsilon_{k}z)|\leq
o(\epsilon_{k})\frac{(1+|z|)}{\log(2+|z|)},\quad|z|<\tau\epsilon_{k}^{-1}.$

From the decomposition in (4.36) we can now estimate the integral of
$\tilde{w}_{l,k}$ more precisely.

(4.48)
$\displaystyle\int_{B(Q_{s}^{k},\tau)}\bigg{(}\tilde{w}_{l,k}(\eta)\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})|\eta|^{2N}e^{\xi_{l}}+\sigma_{k}\nabla\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})(\eta-
Q_{l}^{k})|\eta|^{2N}e^{V_{l,k}}$
$\displaystyle+\frac{\delta_{k}^{2}}{M_{k}}\sum_{|\alpha|=2}\frac{\partial^{\alpha}\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})}{\alpha!}(\eta-
Q_{l}^{k})^{\alpha}|\eta|^{2N})d\eta$
$\displaystyle=\int_{B(Q_{s}^{k},\tau)}\tilde{w}_{l,k}(\eta)\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})|\eta|^{2N}e^{V_{s,k}}d\eta+8\pi\sigma_{k}\nabla\log\mathfrak{h}_{k}(\delta_{k}Q_{s}^{k})(Q_{s}^{k}-Q_{l}^{k})$
$\displaystyle\quad+2\pi\frac{\delta_{k}^{2}}{M_{k}}\frac{\Delta\mathfrak{h}_{k}(\delta_{k}Q_{s}^{k})}{\mathfrak{h}_{k}(\delta_{k}Q_{s}^{k})}|Q_{s}^{k}-Q_{l}^{k}|^{2}+o(\epsilon_{k})$
$\displaystyle=D_{s,l}^{k}+o(\epsilon_{k}).$

where

$D_{s,l}^{k}=\frac{\pi}{(N+1)^{2}}\frac{|p_{s}^{k}-p_{l}^{k}|^{2}}{\epsilon_{k}^{2}M_{k}^{2}}M_{k}+2\pi\frac{\delta_{k}^{2}}{M_{k}}\frac{\Delta\mathfrak{h}_{k}(\delta_{k}Q_{s}^{k})}{\mathfrak{h}_{k}(\delta_{k}Q_{s}^{k})}|Q_{s}^{k}-Q_{l}^{k}|^{2}$

Let

$H_{y,l}(\eta)=\frac{1}{2\pi}\log\frac{|y-\eta|}{|Q_{l}^{k}-\eta|}.$

Then

$\displaystyle\tilde{w}_{l,k}(y)$ $\displaystyle=-\sum_{s\neq
l}H_{y,l}(Q_{s})D_{s,l}^{k}$ $\displaystyle-\sum_{s\neq
l}\int_{B(Q_{s},\tau)}\bigg{(}\partial_{1}H_{y,l}(Q_{s})\eta_{1}+\partial_{2}H_{y,l}(Q_{s})\eta_{2})(\mathfrak{h}_{k}(\delta_{k}Q_{l})|\eta|^{2N}e^{\xi_{l}}\tilde{w}_{l,k}(\eta)$
$\displaystyle+$
$\displaystyle\sigma_{k}\nabla\mathfrak{h}_{k}(\delta_{k}Q_{l})(\eta-
Q_{l})|\eta|^{2N}e^{V_{l,k}}+O(\delta_{k}^{2}M_{k}^{-1})|\eta-
Q_{l}|^{2}|\eta|^{2N}e^{V_{l,k}}\bigg{)}d\eta+o(\epsilon_{k}).$

After evaluation we have

$\displaystyle\tilde{w}_{l,k}(y)$ $\displaystyle=-\frac{1}{2\pi}\sum_{s\neq
l}\log\frac{|y-Q_{s}^{k}|}{|Q_{l}^{k}-Q_{s}^{k}|}D_{s,l}^{k}$
$\displaystyle-\sum_{s\neq
l}8(\frac{y_{1}-Q_{s}^{1}}{|y-Q_{s}|^{2}}-\frac{Q_{l}^{1}-Q_{s}^{1}}{|Q_{l}-Q_{s}|^{2}})c_{1,s,l}\epsilon_{k}$
$\displaystyle-8(\frac{y_{2}-Q_{s}^{2}}{|y-Q_{s}|^{2}}-\frac{Q_{l}^{2}-Q_{s}^{2}}{|Q_{s}-Q_{l}|^{2}})c_{2,s,l}\epsilon_{k}+o(\epsilon_{k}).$

where we used

$\int_{\mathbb{R}^{2}}\frac{z_{1}^{2}}{(1+\frac{1}{8}|z|^{2})^{3}}dz=\int_{\mathbb{R}^{2}}\frac{z_{2}^{2}}{(1+\frac{1}{8}|z|^{2})^{3}}dz=16\pi.$

Recall that

$c_{1,s,l}=\frac{|p_{s}-p_{l}|}{2(N+1)M_{k}\epsilon_{k}}\cos(\beta_{s}+\theta_{sl})$
$c_{2,s,l}=\frac{|p_{s}-p_{l}|}{2(N+1)M_{k}\epsilon_{k}}\sin(\beta_{s}+\theta_{sl})$

For $|y|\sim 1$ but away from the $N+1$ bubbling disks, we have, for $l\neq
s$,

$v_{k}(y)=V_{l,k}(y)+M_{k}\tilde{w}_{l,k}(y)$

and

$v_{k}(y)=V_{s,k}(y)+M_{k}\tilde{w}_{s,k}(y).$

Thus for $s\neq l$ we have

(4.49)
$\frac{V_{s,k}(y)-V_{l,k}(y)}{M_{k}}=\tilde{w}_{l,k}(y)-\tilde{w}_{s,k}(y).$

For $|y|\sim 1$ away from bubbling disks, we have

$V_{s,k}(y)-V_{l,k}(y)=4\log|y^{N+1}-e_{1}-p_{l}|-4\log|y^{N+1}-e_{1}-p_{s}|+o(\epsilon_{k}M_{k}).$
$\displaystyle|y^{N+1}-e_{1}-p_{l}|^{2}$ $\displaystyle=$
$\displaystyle|y^{N+1}-e_{1}-p_{s}|^{2}+2Re((y^{N+1}-e_{1}-p_{s})(\bar{p}_{s}-\bar{p}_{l}))+|p_{l}-p_{s}|^{2}.$

Thus

$\displaystyle\frac{V_{s,k}(y)-V_{l,k}(y)}{M_{k}}$ $\displaystyle=$
$\displaystyle
4Re(\frac{y^{N+1}-1}{|y^{N+1}-1|^{2}}\frac{\bar{p}_{l}-\bar{p}_{s}}{M_{k}\epsilon_{k}})\epsilon_{k}+o(\epsilon_{k}).$

We are going to derive a contradiction based on (4.49). For this purpose we
choose $s=0$ in (4.49), which means $\tilde{w}_{s,k}=o(\epsilon_{k})$, then we
have

(4.50) $\displaystyle
4Re(\frac{y^{N+1}-1}{|y^{N+1}-1|^{2}}\frac{\bar{p}_{l}}{\epsilon_{k}M_{k}})\epsilon_{k}+o(\epsilon_{k})$
$\displaystyle=$ $\displaystyle-\frac{1}{2\pi}\sum_{s\neq
l}\log\frac{|y-Q_{s}^{k}|}{|Q_{l}^{k}-Q_{s}^{k}|}D_{s,l}^{k}$
$\displaystyle+\sum_{s\neq
l}8\bigg{(}(\frac{y_{1}-Q_{l}^{1}}{|y-Q_{l}|^{2}}-\frac{Q_{l}^{1}-Q_{s}^{1}}{|Q_{l}-Q_{s}|^{2}})(\frac{|p_{s}-p_{l}|}{2(N+1)M_{k}\epsilon_{k}})\cos(\beta_{s}+\theta_{sl})$
$\displaystyle+(\frac{y_{2}-Q_{l}^{2}}{|y-Q_{l}|^{2}}-\frac{Q_{l}^{2}-Q_{s}^{2}}{|Q_{l}-Q_{s}|^{2}})\frac{|p_{s}-p_{l}|}{2(N+1)\epsilon_{k}M_{k}}\sin(\beta_{s}+\theta_{sl})\bigg{)}\epsilon_{k},\quad\forall
l.$

If $M_{k}\geq C\epsilon_{k}$, the first term on the right hand side will
dominate all other terms when $|y|>>1$, this violates the equality in (4.50).
But when we consider $M_{k}=o(\epsilon_{k})$, the equality still cannot hold
if choose $|y|>>1$ becase the last two terms of (4.50) will majorize the left-
hand-side. Lemma 4.3 is established. $\Box$

Proposition 4.1 is an immediate consequence of Lemma 4.3. $\Box$.

Now we finish the proof of Theorem 4.1.

Let $\hat{w}_{k}=w_{k}/\tilde{\delta}_{k}$. (Recall that
$\tilde{\delta}_{k}=\delta_{k}|\nabla\mathfrak{h}_{k}(0)|+\delta_{k}^{2}\mu_{k}$).
If $|\nabla\mathfrak{h}_{k}(0)|/(\delta_{k}\mu_{k})\to\infty$, we see that in
this case
$\tilde{\delta}_{k}\sim\delta_{k}\mu_{k}|\nabla\mathfrak{h}_{k}(0)|$. The
equation of $\hat{w}_{k}$ is

(4.51)
$\Delta\hat{w}_{k}+|y|^{2N}e^{\xi_{k}}\hat{w}_{k}=a_{k}\cdot(e_{1}-y)|y|^{2N}e^{V_{k}}+b_{k}e^{V_{k}}|y-e_{1}|^{2}|y|^{2N},$

in $\Omega_{k}$, where
$a_{k}=\delta_{k}\nabla\mathfrak{h}_{k}(0)/\tilde{\delta}_{k}$,
$b_{k}=o(1/\mu_{k})$. By Proposition 4.1, $|\hat{w}_{k}(y)|\leq C$. Before we
carry out the remaining part of the proof we observe that $\hat{w}_{k}$
converges to a harmonic function in $\mathbb{R}^{2}$ minus finite singular
points. Since $\hat{w}_{k}$ is bounded, all these singularities are removable.
Thus $\hat{w}_{k}$ converges to a constant. Based on the information around
$e_{1}$, we shall prove that this constant is $0$. However, looking at the
right hand side the equation,

$(e_{1}-y)|y|^{2N}e^{V_{k}}\rightharpoonup\sum_{l=1}^{N}8\pi(e_{1}-e^{i\beta_{l}})\delta_{e^{i\beta_{l}}}.$

we will get a contradiction by comparing the Pohozaev identities of $v_{k}$
and $V_{k}$, respectively.

Now we use the notation $W_{k}$ again and use Proposition 4.1 to rewrite the
equation for $W_{k}$. Let

$W_{k}(z)=\hat{w}_{k}(e_{1}+\epsilon_{k}z),\quad|z|<\delta_{0}\epsilon_{k}^{-1}$

for $\delta_{0}>0$ small. Then from Proposition 4.1 we have

(4.52)
$\mathfrak{h}_{k}(\delta_{k}y)=\mathfrak{h}_{k}(\delta_{k}e_{1})+\delta_{k}\nabla\mathfrak{h}_{k}(\delta_{k}e_{1})(y-e_{1})+O(\delta_{k}^{2})|y-e_{1}|^{2},$
(4.53) $|y|^{2N}=|e_{1}+\epsilon_{k}z|^{2N}=1+O(\epsilon_{k})|z|,$ (4.54)
$V_{k}(e_{1}+\epsilon_{k}z)+2\log\epsilon_{k}=U_{k}(z)+O(\epsilon_{k})|z|+O(\epsilon_{k}^{2})(\log(1+|z|))^{2}$

and

(4.55)
$\xi_{k}(e_{1}+\epsilon_{k}z)+2\log\epsilon_{k}=U_{k}(z)+O(\epsilon_{k})(1+|z|).$

Using (4.52),(4.53),(4.54) and (4.55) in (4.51) we write the equation of
$W_{k}$ as

(4.56) $\Delta
W_{k}+\mathfrak{h}_{k}(\delta_{k}e_{1})e^{U_{k}(z)}W_{k}=-\epsilon_{k}a_{k}\cdot
ze^{U_{k}(z)}+E_{w},\quad 0<|z|<\delta_{0}\epsilon_{k}^{-1}$

where

(4.57)
$E_{w}(z)=O(\epsilon_{k})(1+|z|)^{-3},\quad|z|<\delta_{0}\epsilon_{k}^{-1}.$

Since $\hat{w}_{k}$ obviously converges to a global harmonic function with
removable singularity, we have $\hat{w}_{k}\to\bar{c}$ for some
$\bar{c}\in\mathbb{R}$. Then we claim that

######  Lemma 4.4.

$\bar{c}=0$.

Proof of Lemma 4.4:

If $\bar{c}\neq 0$, we use $W_{k}(z)=\bar{c}+o(1)$ on
$B(0,\delta_{0}\epsilon_{k}^{-1})\setminus
B(0,\frac{1}{2}\delta_{0}\epsilon_{k}^{-1})$ and consider the projection of
$W_{k}$ on $1$:

$g_{0}(r)=\frac{1}{2\pi}\int_{0}^{2\pi}W_{k}(re^{i\theta})d\theta.$

If we use $F_{0}$ to denote the projection to $1$ of the right hand side we
have, using the rough estimate of $E_{w}$ in (4.57)

$g_{0}^{\prime\prime}(r)+\frac{1}{r}g_{0}^{\prime}(r)+\mathfrak{h}_{k}(\delta_{k}e_{1})e^{U_{k}(r)}g_{0}(r)=F_{0},\quad
0<r<\delta_{0}\epsilon_{k}^{-1}$

where

$F_{0}(r)=O(\epsilon_{k})(1+|z|)^{-3}.$

In addition we also have

$\lim_{k\to\infty}g_{0}(\delta_{0}\epsilon_{k}^{-1})=\bar{c}+o(1).$

For simplicity we omit $k$ in some notations. By the same argument as in Lemma
4.1, we have

$g_{0}(r)=O(\epsilon_{k})\log(2+r),\quad 0<r<\delta_{0}\epsilon_{k}^{-1}.$

Thus $\bar{c}=0$. Lemma 4.4 is established. $\Box$

Based on Lemma 4.4 and standard Harnack inequality for elliptic equations we
have

(4.58) $\tilde{w}_{k}(x)=o(1),\,\,\nabla\tilde{w}_{k}(x)=o(1),\,\,x\in
B_{3}\setminus(\cup_{l=1}^{N}(B(e^{i\beta_{l}},\delta_{0})\setminus
B(e^{i\beta_{l}},\delta_{0}/8))).$

Equation (4.58) is equivalent to $w_{k}=o(\tilde{\delta}_{k})$ and $\nabla
w_{k}=o(\tilde{\delta}_{k})$ in the same region.

In the next step we consider the difference between two Pohozaev identities.
For $s=1,...,N$ we consider the Pohozaev identity around $Q_{s}^{k}$. Let
$\Omega_{s,k}=B(Q_{s}^{k},r)$ for small $r>0$. For $v_{k}$ we have

(4.59)
$\displaystyle\int_{\Omega_{s,k}}\partial_{\xi}(|y|^{2N}\mathfrak{h}_{k}(\delta_{k}y))e^{v_{k}}-\int_{\partial\Omega_{s,k}}e^{v_{k}}|y|^{2N}\mathfrak{h}_{k}(\delta_{k}y)(\xi\cdot\nu)$
$\displaystyle=\int_{\partial\Omega_{s,k}}(\partial_{\nu}v_{k}\partial_{\xi}v_{k}-\frac{1}{2}|\nabla
v_{k}|^{2}(\xi\cdot\nu))dS.$

where $\xi$ is an arbitrary unit vector. Correspondingly the Pohozaev identity
for $V_{k}$ is

(4.60)
$\displaystyle\int_{\Omega_{s,k}}\partial_{\xi}(|y|^{2N}\mathfrak{h}_{k}(\delta_{k}e_{1}))e^{V_{k}}-\int_{\partial\Omega_{s,k}}e^{V_{k}}|y|^{2N}\mathfrak{h}_{k}(\delta_{k}e_{1})(\xi\cdot\nu)$
$\displaystyle=\int_{\partial\Omega_{s,k}}(\partial_{\nu}V_{k}\partial_{\xi}V_{k}-\frac{1}{2}|\nabla
V_{k}|^{2}(\xi\cdot\nu))dS.$

Using $w_{k}=v_{k}-V_{k}$ and $|w_{k}(y)|\leq C\tilde{\delta}_{k}$ we have

$\displaystyle\int_{\partial\Omega_{s,k}}(\partial_{\nu}v_{k}\partial_{\xi}v_{k}-\frac{1}{2}|\nabla
v_{k}|^{2}(\xi\cdot\nu))dS$ $\displaystyle=$
$\displaystyle\int_{\partial\Omega_{s,k}}(\partial_{\nu}V_{k}\partial_{\xi}V_{k}-\frac{1}{2}|\nabla
V_{k}|^{2}(\xi\cdot\nu))dS$
$\displaystyle+\int_{\partial\Omega_{s,k}}(\partial_{\nu}V_{k}\partial_{\xi}w_{k}+\partial_{\nu}w_{k}\partial_{\xi}V_{k}-(\nabla
V_{k}\cdot\nabla w_{k})(\xi\cdot\nu))dS+o(\tilde{\delta}_{k}).$

If we just use crude estimate: $\nabla w_{k}=o(\tilde{\delta}_{k})$, we have

$\displaystyle\int_{\partial\Omega_{s,k}}(\partial_{\nu}v_{k}\partial_{\xi}v_{k}-\frac{1}{2}|\nabla
v_{k}|^{2}(\xi\cdot\nu))dS$ $\displaystyle-$
$\displaystyle\int_{\partial\Omega_{s,k}}(\partial_{\nu}V_{k}\partial_{\xi}V_{k}-\frac{1}{2}|\nabla
V_{k}|^{2}(\xi\cdot\nu))dS=o(\tilde{\delta}_{k}).$

The difference on the second terms is minor: If we use the expansion of
$v_{k}=V_{k}+w_{k}$ and that of $\mathfrak{h}_{k}(\delta_{k}y)$ around
$e_{1}$, it is easy to obtain

$\int_{\partial\Omega_{s,k}}e^{v_{k}}|y|^{2N}\mathfrak{h}_{k}(\delta_{k}y)(\xi\cdot\nu)-\int_{\partial\Omega_{s,k}}e^{V_{k}}|y|^{2N}\mathfrak{h}_{k}(\delta_{k}e_{1})(\xi\cdot\nu)=o(\tilde{\delta}_{k}).$

To evaluate the first term, we use

(4.61)
$\displaystyle\partial_{\xi}(|y|^{2N}\mathfrak{h}_{k}(\delta_{k}y))e^{v_{k}}$
$\displaystyle=$
$\displaystyle\partial_{\xi}(|y|^{2N}\mathfrak{h}_{k}(\delta_{k}e_{1})+|y|^{2N}\delta_{k}\nabla\mathfrak{h}_{k}(\delta_{k}e_{1})(y-e_{1})+O(\delta_{k}^{2}))e^{V_{k}}(1+w_{k}+O(\delta_{k}^{2}\mu_{k}))$
$\displaystyle=$
$\displaystyle\partial_{\xi}(|y|^{2N})\mathfrak{h}_{k}(\delta_{k}e_{1})e^{V_{k}}+\delta_{k}\partial_{\xi}(|y|^{2N}\nabla\mathfrak{h}_{k}(\delta_{k}e_{1})(y-e_{1}))e^{V_{k}}$
$\displaystyle+\partial_{\xi}(|y|^{2N}\mathfrak{h}_{k}(\delta_{k}e_{1}))e^{V_{k}}w_{k}+O(\delta_{k}^{2}\mu_{k})e^{V_{k}}.$

For the third term on the right hand side of (4.61) we use the equation for
$w_{k}$:

$\Delta
w_{k}+\mathfrak{h}_{k}(\delta_{k}e_{1})e^{V_{k}}|y|^{2N}w_{k}=-\delta_{k}\nabla\mathfrak{h}_{k}(\delta_{k}e_{1})\cdot(y-e_{1})|y|^{2N}e^{V_{k}}+O(\delta_{k}^{2})e^{V_{k}}|y|^{2N}.$

From integration by parts we have

$\displaystyle\int_{\Omega_{s,k}}\partial_{\xi}(|y|^{2N})\mathfrak{h}_{k}(\delta_{k}e_{1})e^{V_{k}}w_{k}$
$\displaystyle=$ $\displaystyle
2N\int_{\Omega_{s,k}}|y|^{2N-2}y_{\xi}\mathfrak{h}_{k}(\delta_{k}e_{1})e^{V_{k}}w_{k}$
$\displaystyle=$ $\displaystyle
2N\int_{\Omega_{s,k}}\frac{y_{\xi}}{|y|^{2}}(-\Delta
w_{k}-\delta_{k}\nabla\mathfrak{h}_{k}(\delta_{k}e_{1})(y-e_{1})|y|^{2N}e^{V_{k}}+O(\delta_{k}^{2})e^{V_{k}}|y|^{2N})$
$\displaystyle=$
$\displaystyle-2N\delta_{k}\int_{\Omega_{s,k}}\frac{y_{\xi}}{|y|^{2}}\nabla\mathfrak{h}_{k}(\delta_{k}e_{1})(y-e_{1})|y|^{2N}e^{V_{k}}$
$\displaystyle+2N\int_{\partial\Omega_{s,k}}(\partial_{\nu}(\frac{y_{\xi}}{|y|^{2}})w_{k}-\partial_{\nu}w_{k}\frac{y_{\xi}}{|y|^{2}})+o(\tilde{\delta}_{k})$
(4.62) $\displaystyle=$
$\displaystyle\nabla\mathfrak{h}_{k}(\delta_{k}e_{1})\bigg{(}-16N\delta_{k}\pi(e^{i\beta_{s}}\cdot\xi)(e^{i\beta_{s}}-e_{1})+O(\mu_{k}\epsilon_{k}^{2})\bigg{)}+o(\tilde{\delta}_{k}),$

where we have used $\nabla w_{k},w_{k}=o(\tilde{\delta}_{k})$ on
$\partial\Omega_{s,k}$. For the second term on the right hand side of (4.61),
we have

(4.63)
$\displaystyle\int_{\Omega_{s,k}}\delta_{k}\partial_{\xi}(|y|^{2N}\nabla\mathfrak{h}_{k}(\delta_{k}e_{1})(y-e_{1}))e^{V_{k}}$
$\displaystyle=$ $\displaystyle
2N\delta_{k}\int_{\Omega_{s,k}}y_{\xi}|y|^{2N-2}\nabla\mathfrak{h}_{k}(\delta_{k}e_{1})(y-e_{1})e^{V_{k}}+\delta_{k}\int_{\Omega_{s,k}}|y|^{2N}\partial_{\xi}\mathfrak{h}_{k}(\delta_{k}e_{1})e^{V_{k}}$
$\displaystyle=$
$\displaystyle\nabla\mathfrak{h}_{k}(\delta_{k}e_{1})\big{(}16N\pi\delta_{k}(e^{i\beta_{s}}\cdot\xi)(e^{i\beta_{s}}-e_{1})+O(\mu_{k}\epsilon_{k}^{2})\big{)}$
$\displaystyle+\delta_{k}\partial_{\xi}\mathfrak{h}_{k}(\delta_{k}e_{1})(8\pi+O(\mu_{k}\epsilon_{k}^{2}))+o(\tilde{\delta}_{k}).$

Using (4) and (4.63) in the difference between (4.59) and (4.60), we have

$\delta_{k}\partial_{\xi}\mathfrak{h}_{k}(\delta_{k}e_{1})(1+O(\mu_{k}\epsilon_{k}^{2}))=o(\tilde{\delta}_{k}).$

Thus $\nabla\mathfrak{h}_{k}(\delta_{k}e_{1})=O(\delta_{k}\mu_{k})$. Theorem
4.1 is established. $\Box$

##  5\. Proof of Theorem 1.1

First we prove the case $N\geq 2$. In [26] we have already proved that

$\Delta(\log\mathfrak{h}_{k})(0)=O(\delta_{k}^{-2}\mu_{k}e^{-\mu_{k}})+O(\delta_{k}).$

So if $\delta_{k}/(\mu_{k}^{\frac{1}{2}}\epsilon_{k})\to\infty$ there is
nothing to prove. So we only consider the case that $\delta_{k}\leq
C\mu_{k}^{\frac{1}{2}}\epsilon_{k}$. By in this case
$\epsilon_{k}^{-1}\delta_{k}^{2}\leq C\epsilon_{k}^{\epsilon}$ for some
$\epsilon\in(0,1)$. The whole argument of Proposition 4.1 can be employed to
prove

(5.1) $|w_{k}(y)|\leq C\delta_{k}^{2}\mu_{k}^{\frac{7}{4}}$

In order to employ the same strategy of proof, one needs to have three things:
first $\epsilon_{k}^{-1}\delta_{k}^{2}=O(\epsilon_{k}^{\epsilon})$. This is
clear from the definition of $\delta_{k}$. Second, in the proof of Lemma 4.3
we need

$O(\delta_{k}^{N+2}/M_{k})=o(\epsilon_{k}),$

where $M_{k}>\delta_{k}^{2}\mu_{k}^{7/4}$. Since $\delta_{k}\leq
C\mu_{k}^{\frac{1}{2}}\epsilon_{k}$ and $N\geq 1$, the required inequality
holds. Thirdly, we need to have

$\frac{\delta_{k}^{3}}{M_{k}}=o(\epsilon_{k}).$

This is used in (7.1). From the requirement on $\delta_{k}$ and the definition
of $M_{k}$ This clearly also holds. The proof of Proposition 4.1 follows. Thus
for $N\geq 2$ we also have (5.1).

The precise upper bound of $w_{k}$ in (5.1) leads to the vanishing rate of the
Laplacian estimate for $N\geq 2$ and some cases of $N=1$: If we use

$W_{k}(z)=w_{k}(e_{l}+\epsilon_{k}z)/(\delta_{k}^{2}\mu_{k}^{\frac{7}{4}}),\quad|z|<\tau\epsilon_{k}^{-1}$

where $e_{l}\neq e_{1}$. We shall show that the projection of $W_{k}$ over $1$
is not bounded when $|z|\sim\epsilon_{k}^{-1}$, which gives the desired
contradiction.

We write the equation of $w_{k}$ as

$\Delta
w_{k}+|y|^{2N}e^{\xi_{k}}w_{k}=(\mathfrak{h}_{k}(\delta_{k}e_{1})-\mathfrak{h}_{k}(\delta_{k}y))|y|^{2N}e^{V_{k}}.$

Then for $l\neq 1$,

$\displaystyle\Delta W_{k}(z)+e^{U_{k}}W_{k}(z)$ $\displaystyle=$
$\displaystyle
a_{0}e^{U_{k}}+a_{1}ze^{U_{k}}+\frac{1}{2\mu_{k}^{7/4}}\Delta\mathfrak{h}_{k}(0)|z|^{2}e^{U_{k}}+\frac{1}{\mu_{k}^{7/4}}R_{2}(\theta)|z|^{2}e^{U_{k}}+O(\epsilon_{k}^{\epsilon}(1+|z|)^{-3}).$

where

$a_{0}=(\mathfrak{h}_{k}(\delta_{k}e_{1})-\mathfrak{h}_{k}(\delta_{k}e_{l}))/(\delta_{k}^{2}\mu_{k}^{7/4}),$
$a_{1}=-\nabla\mathfrak{h}_{k}(\delta_{k}e_{l})/(\delta_{k}\mu_{k}^{7/4}),$

$R_{2}$ is the collection of spherical harmonic functions of degree $2$. Note
that the assumption $l\neq 1$ means there is no appearance of $\epsilon_{k}$
or $\epsilon_{k}^{2}$ in the equation for $W_{k}$.

Let $g_{k}(r)$ be the projection of $W_{k}$ on $1$, by the same ODE analysis
as before, we see that $g_{k}$ satisfies

$g_{k}^{\prime\prime}+\frac{1}{r}g_{k}^{\prime}(r)+e^{U_{k}}g_{k}=E_{k}$

where

$E_{k}(r)=O(\epsilon_{k}^{\epsilon})(1+r)^{-3}+\frac{1}{2\mu_{k}^{7/4}}\Delta(\log\mathfrak{h}_{k})(0)r^{2}e^{U_{k}}.$

Using the same argument as in Lemma 4.1, we have

$g_{k}(r)\sim\Delta(\log\mathfrak{h}_{k})(0)(\log r)^{2}\mu_{k}^{-7/4},\quad
r>10.$

Clearly if $\Delta(\log\mathfrak{h}_{k}(0)\neq 0$ we obtain a violation of the
bound of $w_{k}$ for $r\sim\epsilon_{k}^{-1}$. Theorem 1.1 for $N\geq 2$ is
proved under the assumption

(5.2)
$\epsilon_{k}^{-1}|Q_{s}^{k}-e^{i\beta_{s}}|\leq\epsilon_{k}^{\epsilon},\quad
s=1,...,N.$

We need this assumption because the $\xi_{k}$ function that comes from the
equation of $w_{k}$ needs to tend to $U$ after scaling. From (3.13) in [26],
$|Q_{s}^{k}-e^{i\beta_{s}}|=O(\delta_{k}^{2})+O(\mu_{k}e^{-\mu_{k}})$. If
$\delta_{k}^{2}\epsilon_{k}^{-1}\geq C$, the argument in Theorem 4.1 cannot be
used because either $\xi_{k}$ does not tend $U$ or $c_{0}=0$ cannot be proved.
For $N\geq 2$, this is not a problem because we only consider $\delta_{k}\leq
C\mu_{k}^{\frac{1}{2}}\epsilon_{k}$.

Next we prove Theorem 1.1 for $N=1$ and $\delta_{k}\leq\mu_{k}\epsilon_{k}$.
The reader could see immediately that the same proof for the case $N\geq 2$
still works. So the only remaining case is

Proof of Theorem 1.1 for $N=1$ and $\delta_{k}\geq\mu_{k}\epsilon_{k}$.

In this case we write the equation of $w_{k}$ as

$\Delta
w_{k}+|y|^{2}\mathfrak{h}_{k}(\delta_{k}y)e^{v_{k}}-|y|^{2}\mathfrak{h}_{k}(\delta_{k}e_{1})e^{V_{k}}=0.$

From $0=\nabla w_{k}(e_{1})$ we have

(5.3)
$0=\int_{\Omega_{k}}\nabla_{1}G_{k}(e_{1},\eta)|\eta|^{2}(\mathfrak{h}_{k}(\delta_{k}\eta)e^{v_{k}}-\mathfrak{h}_{k}(\delta_{k}e_{1})e^{V_{k}})d\eta+O(\delta_{k}^{3})$

Note that $v_{k}$ is close to another global solution $\bar{V}_{k}$ which
matches with a local maximum of $v_{k}$ at $Q_{2}^{k}$. Evaluating the right
hand side of (5.3) we have

$\nabla_{1}G_{k}(e_{1},Q_{2}^{k})-\nabla_{1}G_{k}(e_{1},e^{i\pi})=O(\epsilon_{k}^{2}\mu_{k})+O(\delta_{k}^{3}).$

This expression gives

$Q_{2}^{k}-e^{i\pi}=O(\delta_{k}^{3})+O(\mu_{k}\epsilon_{k}^{2}).$

This estimate will lead to a better estimate of $w_{k}$ outside the two
bubbling disks. From the Green’s representation for $w_{k}$ we now obtain

$w_{k}(y)=\int_{\Omega_{k}}(G_{k}(y,\eta)-G_{k}(e_{1},\eta))|\eta|^{2}(\mathfrak{h}_{k}(\delta_{k}\eta)e^{v_{k}(y)}-\mathfrak{h}_{k}(\delta_{k}e_{1})e^{V_{k}})d\eta+O(\delta_{k}^{2})$

where the last term $O(\delta_{k}^{2})$ comes from the oscillation of $w_{k}$
on $\partial\Omega_{k}$. Then we have

$\displaystyle w_{k}(y)$
$\displaystyle=-\frac{1}{2\pi}\int_{\Omega_{k}}\log\frac{|y-\eta|}{|e_{1}-\eta|}|\eta|^{2}(\mathfrak{h}_{k}(\delta_{k}\eta)-\mathfrak{h}_{k}(\delta_{k}e_{1})e^{V_{k}})+O(\delta_{k}^{2})$
$\displaystyle=-4\log\frac{|y-Q_{2}^{k}|}{|e_{1}-Q_{2}^{k}|}+4\log\frac{|y-e^{i\pi}|}{2}+O(\delta_{k}^{2}\mu_{k}).$

By $|Q_{2}^{k}-e^{i\pi}|=O(\delta_{k}^{2})$ we see that
$w_{k}(y)=O(\delta_{k}^{2})$ on $|y-e^{i\pi}|=\tau$.

The standard point-wise estimate for singular equation ( see [30, 17] ) gives

$\displaystyle v_{k}(Q_{2}^{k}+\epsilon_{k}z)+2\log\epsilon_{k}$
$\displaystyle=$
$\displaystyle\log\frac{e^{\mu_{k}}}{(1+\frac{e^{\mu_{k}}}{8\mathfrak{h}_{k}(\delta_{k}Q_{k})}|z|^{2})^{2}}+\phi_{1}^{k}+C\delta_{k}^{2}\Delta(\log\mathfrak{h}_{k})(0)(\log(1+|z|))^{2},\quad|z|\sim\epsilon_{k}^{-1}.$

$\displaystyle V_{k}(e^{i\pi}+\epsilon_{k}z)+2\log\epsilon_{k}$
$\displaystyle=$
$\displaystyle\log\frac{e^{\mu_{k}}}{(1+\frac{e^{\mu_{k}}}{8\mathfrak{h}_{k}(\delta_{k}e_{1})}|z|^{2})^{2}}+\phi_{2}^{k}+O(\epsilon_{k}^{2}(\log\epsilon_{k})^{2}),\quad|z|\sim\epsilon_{k}^{-1}.$

Thus

$w_{k}(Q_{2}^{k}+\epsilon_{k}z)=O(\epsilon_{k}^{2}(\log\epsilon_{k})^{2})+\phi_{1}^{k}-\phi_{2}^{k}+C\Delta(\log\mathfrak{h}_{k})(0)\delta_{k}^{2}(\log(1+|z|)^{2}),$

for $|z|\sim\epsilon_{k}^{-1}.$ Taking the average around the origin, the
spherical averages of the two harmonic functions are zero and
$O(\delta_{k}^{2})$ respectively, since they take zero at the origin and a
point at most $O(\delta_{k}^{2})$ from the origin. So the spherical average of
$w_{k}$ is comparable to

$\Delta(\log\mathfrak{h}_{k})(0)\delta_{k}^{2}(\log\epsilon_{k})^{2}$

for $|z|\sim\epsilon_{k}^{-1}$. Thus we know
$\Delta(\log\mathfrak{h}_{k})(0)=o(1)$ because
$w_{k}=O(\delta_{k}^{2}\mu_{k})$ in this region, Theorem 1.1 is established
for all the cases. $\Box$

##  6\. Singular Mean field equation

In this section we prove Theorem 1.2. First it is well known that if $p$ is a
blowup point that has a non-quantized singular source ($\alpha_{p}=0$ or
$\alpha_{p}\not\in\mathbb{N}$), the profile of the bubbling solutions around
$p$ is a simple blwoup (see [30, 31]). So we only need to focus on the case
that $\alpha_{p}\in\mathbb{N}$. Let $G(\cdot,\cdot)$ be the Green’s function
corresponding to $-\Delta_{g}$:

$-\Delta_{y}G(p,y)=\delta_{p}-1,\quad\int_{M}G(p,y)dV_{g}=0.$

By setting

$G_{1}(y)=4\pi\sum_{t=1}^{M}\alpha_{t}G(p_{t},y),$

we have

$-\Delta G_{1}=4\pi\sum_{t}\alpha_{t}(\delta_{p_{t}}-1).$

Then the function $v_{k}=u_{k}+G_{1}$ satisfies

$\Delta_{g}v_{k}+\rho^{k}(\frac{he^{v_{k}}e^{-G_{1}}}{\int_{M}he^{u_{k}}}-1)=0.$

Let $p_{1}$ be a quantized singular source, which means $\alpha_{p_{1}}\in
4\pi\mathbb{N}$, in the neighborhood of $p_{1}$ we have

$\Delta_{g}v_{k}+|y-p_{1}|^{2\alpha_{p_{1}}}H_{k}e^{v_{k}}=\rho_{k}$

where

$H_{k}=\frac{-\rho_{k}he^{4\pi\alpha_{1}\gamma(p_{1},y)-4\pi\sum_{t\neq
1}\alpha_{t}G(p_{t},y)}}{\int_{M}he^{u_{k}}},$

where $\gamma$ is the regular part of the Green’s function. In local
coordinates around $p_{1}$, the equation can be written as

$\Delta v_{k}+|x|^{2\alpha_{p_{1}}}He^{\phi}e^{v_{k}}=\rho e^{\phi}.$

where $\phi(0)=|\nabla\phi(0)|=0$ and $\Delta\phi(0)=-2K(p_{1})$. Finally we
use $f$ to remove the right hand side:

$\Delta f=\rho^{k}e^{\phi},\quad f(0)=0,\quad f=\mbox{constant on }\partial
B_{\tau}$

for $\tau>0$ small. When we consider $v^{k}-f$ as the blowup solutions, we
have

$\Delta(v_{k}-f)+|y|^{2\alpha_{p_{1}}}H_{k}e^{f}e^{v_{k}-f}=0.$

It is a standard result that $H_{k}$ is uniformly bounded above and below.
From the defintion of $H_{k}$ we have

$\Delta H_{k}(0)=\Delta h(p_{1})-\sum_{t=1}^{M}4\pi\alpha_{t}.$

Using Theorem 1.1 we would have

$\Delta\log H_{k}(0)+\Delta\phi+\Delta f=o(1).$

if non-simple blowup happens at $p_{1}$, which is

(6.1) $\Delta\log
h(p_{1})-2K(p_{1})-4\pi\sum_{t=1}^{M}\alpha_{t}+\rho^{k}=o(1).$

Since $\rho^{k}\to\rho\in\Lambda$, we see from (1.8) that (6.1) cannot hold.
Theorem 1.2 is established. $\Box$

##  7\. Appendix

In this section we prove (4.46). The terms of $\phi_{1}$ and $\phi_{3}$ lead
to $o(\epsilon_{k})$, the integration involving $\phi_{2}$ cancels with the
second term of (4.46). The computation of $\phi_{2}$ is based on this
equation:

$\int_{\mathbb{R}^{2}}\frac{\frac{\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})}{4}\sigma_{k}\nabla\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})(Q_{l}^{k}-Q_{s}^{k})|z|^{2}}{(1+\frac{\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})}{8}|z|^{2})^{3}}dz=8\pi\sigma_{k}\nabla(\log\mathfrak{h}_{k})(\delta_{k}Q_{l}^{k})(Q_{l}^{k}-Q_{s}^{k}),$

and by (4.2)

(7.1)
$\nabla\log\mathfrak{h}_{k}(\delta_{k}Q_{l}^{k})-\nabla\log\mathfrak{h}_{k}(\delta_{k}Q_{s}^{k})=O(\delta_{k})=o(\epsilon_{k}).$

The integration involving $\phi_{4}$ provides the leading term. More detailed
information is the following: First for a global solution

$V_{\mu,p}=\log\frac{e^{\mu}}{(1+\frac{e^{\mu}}{\lambda}|z^{N+1}-p|^{2})^{2}}$

of

$\Delta
V_{\mu,p}+\frac{8(N+1)^{2}}{\lambda}|z|^{2N}e^{V_{\mu,p}}=0,\quad\mbox{in
}\quad\mathbb{R}^{2},$

by differentiation with respect to $\mu$ we have

$\Delta(\partial_{\mu}V_{\mu,p})+\frac{8(N+1)^{2}}{\lambda}|z|^{2N}e^{V_{\mu,p}}\partial_{\mu}V_{\mu,p}=0,\quad\mbox{in}\quad\mathbb{R}^{2}.$

By the expression of $V_{\mu,p}$ we see that

$\partial_{r}\bigg{(}\partial_{\mu}V_{\mu,p}\bigg{)}(x)=O(|x|^{-2N-3}).$

Thus we have

(7.2)
$\int_{\mathbb{R}^{2}}\partial_{\mu}V_{\mu,p}|z|^{2N}e^{V_{\mu,p}}=\int_{\mathbb{R}^{2}}\frac{(1-\frac{e^{\mu}}{\lambda}|z^{N+1}-P|^{2})|z|^{2N}}{(1+\frac{e^{\mu}}{\lambda}|z^{N+1}-P|^{2})^{3}}dz=0.$

From $V_{\mu,p}$ we also have

$\int_{\mathbb{R}^{2}}\partial_{P}V_{\mu,p}|y|^{2N}e^{V_{\mu,p}}=\int_{\mathbb{R}^{2}}\partial_{\bar{P}}V_{\mu,p}|y|^{2N}e^{V_{\mu,p}}=0,$

which gives

(7.3)
$\int_{\mathbb{R}^{2}}\frac{\frac{e^{\mu}}{\lambda}(\bar{z}^{N+1}-\bar{P})|z|^{2N}}{(1+\frac{e^{\mu}}{\lambda}|z^{N+1}-P|^{2})^{3}}=\int_{\mathbb{R}^{2}}\frac{\frac{e^{\mu}}{\lambda}(z^{N+1}-P)|z|^{2N}}{(1+\frac{e^{\mu}}{\lambda}|z^{N+1}-P|^{2})^{3}}=0.$

Now we need more precise expressions of $\phi_{1}$, $\phi_{3}$ and $B$:

$\displaystyle\phi_{1}=(\mu_{s}^{k}-\mu_{l}^{k})(1-\frac{(N+1)^{2}}{D_{s}^{k}}|z+\frac{N}{2}\epsilon_{k}z^{2}e^{-i\beta_{s}}|^{2}/B,$
$\displaystyle\phi_{3}=\frac{4(N+1)}{D_{s}^{k}B}Re((z+\frac{N}{2}\epsilon_{k}e^{-i\beta_{s}}z^{2}))(\frac{\bar{p}_{s}^{k}-\bar{p}_{l}^{k}}{\epsilon_{k}}e^{-i\beta_{s}}))$
$\displaystyle
B=1+\frac{(N+1)^{2}}{D_{s}^{k}}|z+\frac{N}{2}z^{2}e^{-i\beta_{s}}\epsilon_{k}|^{2},$

From here we use scaling and cancellation to have

$\int_{B(0,\tau\epsilon_{k}^{-1})}\frac{\phi_{1}}{M_{k}}B^{-2}=o(\epsilon_{k}),\quad\int_{B(0,\tau\epsilon_{k}^{-1})}\frac{\phi_{3}}{M_{k}}B^{-2}=o(\epsilon_{k}).$

Thus (4.46) holds.

## References

  * [1] M. Ahmedou, L. Wu, Lei Zhang, Classification and a priori estimate for singular prescribing Q-curvature equations on 4-manifolds, preprint, 2022, https://arxiv.org/abs/2012.11785
  * [2] J. Ambjorn and P. Olesen, A condensate solution of the electroweak theory which interpolates between the broken and the symmetry phase, Nucl. Phys. B 330 (1990), 193-204.
  * [3] Bartolucci, Daniele; Jevnikar, Aleks; Lee, Youngae; Yang, Wen Uniqueness of bubbling solutions of mean field equations. J. Math. Pures Appl.  (9) 123 (2019), 78–126.
  * [4] D. Bartolucci, G. Tarantello, Asymptotic Blow-up Analysis for Singular Liouville type Equations with Applications, J. Diff. Eq.  262, n. 7 (2017), 3887-3931
  * [5] Bartolucci, ; Tarantello, G. The Liouville equation with singular data: a concentration-compactness principle via a local representation formula. J. Differential Equations 185 (2002), no. 1, 161–180.
  * [6] D. Bartolucci, G. Tarantello, Liouville type equations with singular data and their applications to periodic multivortices for the Electroweak theory. Comm. Math. Physics  229 (2002), 3-47,
  * [7] Caffarelli, Luis A.; Gidas, Basilis; Spruck, Joel Asymptotic symmetry and local behavior of semilinear elliptic equations with critical Sobolev growth. Comm. Pure Appl. Math. 42 (1989), no. 3, 271–297.
  * [8] E. Caglioti, P.L. Lions, C. Marchioro , M. Pulvirenti, A special class of stationary flows for two-dimensional Euler equations: A statistical mechanics description, Comm. Math. Phys. 143 (1992), 501-525.
  * [9] E. Caglioti, P.L. Lions, C. Marchioro , M. Pulvirenti, A special class of stationary flows for two-dimensional Euler equations: A statistical mechanics description, part II, Comm. Math. Phys. 174 (1995), 229-260.
  * [10] Chanillo, Sagun; Kiessling, Michael, Rotational symmetry of solutions of some nonlinear problems in statistical mechanics and in geometry. Comm. Math. Phys.  160 (1994), no. 2, 217-238.
  * [11] Li, Yan Yan Harnack type inequality: the method of moving planes. Comm. Math. Phys. 200 (1999), no. 2, 421–444.
  * [12] C. C. Chen, C. S. Lin, Sharp estimates for solutions of multi-bubbles in compact Riemann surfaces. Comm. Pure Appl. Math.  55 (2002), no. 6, 728-771.
  * [13] Chen, Chiun-Chuan; Lin, Chang-Shou Mean field equation of Liouville type with singular data: topological degree. Comm. Pure Appl. Math. 68 (2015), no. 6, 887-947.
  * [14] M. Del Pino, P. Esposito , M. Musso, Nondegeneracy of entire solutions of a singular Liouville equation, Proc. of AMS  140 (2012), no. 2 581-588.
  * [15] M. del Pino, P. Esposito, M. Musso, Two-dimensional Euler flows with concentrated vorticities, Trans. Amer. Math. Soc. 362 (2010), no. 12, 6381-6395.
  * [16] A. Eremenko, A. Gabrielov, V. Tarasov, Metrics with conic singularities and spherical polygons, Illinois J. Math, vol 58, no3, 2014, 739-755.
  * [17] M. Gluck, Asymptotic behavior of blow up solutions to a class of prescribing Gauss curvature equations. Nonlinear Anal. 75 (2012), no. 15, 5787-5796.
  * [18] Y. Gu, L. Zhang, Degree counting theorems for singular Liouville systems, to appear on Annali della Scuola Normale Superiore di Pisa, Classe di Scienze (5) Vol. XXI (2020), 1103-1135.
  * [19] Y. Gu, L. Zhang, Structure of bubbling solutions of Liouville systems with negative singular sources, preprint 2022, https://arxiv.org/abs/2112.10031
  * [20] R. Jackiw , E.J. Weinberg, Self-Dual Chern Simons vortices, Phys. Rev. Lett. 64 (1990), 2234-2237.
  * [21] Kuo, Ting-Jung; Lin, Chang-Shou Estimates of the mean field equations with integer singular sources: non-simple blowup. J. Differential Geom. 103 (2016), no. 3, 377-424.
  * [22] C. S. Lin, S. Yan, On the mean field type bubbling solutions for Chern–Simons–Higgs equation, Adv. Math. 338 (2018) 1141–1188.
  * [23] Prajapat, J.; Tarantello, G. On a class of elliptic problems in R2: symmetry and uniqueness results. Proc. Roy. Soc. Edinburgh Sect. A 131 (2001), no. 4, 967-985.
  * [24] M. Troyanov, Prescribing curvature on compact surfaces with conical singularities, Trans. AMS, 324 (1991) 793-821.
  * [25] M. Troyanov, Metrics of constant curvature on a sphere with two conical singularities, Lect. Notes Math., 1410, Springer, NY, 1989, 296-308.
  * [26] Juncheng Wei, Lei Zhang, Estimates for Liouville equation with quantized singularities. Adv. Math.  380 (2021), 107606.
  * [27] Juncheng Wei, Lei Zhang, Vanishing estimates for Liouville equation with quantized singularities, Proceedings of the London Mathematical Society. (3) 2022;1-26
  * [28] Lina Wu, Lei Zhang, Uniqueness of bubbling solutions of mean field equations with non-quantized singularities, Communications in Contemporary Mathematics,23 (2021), no. 4, 2050038, 39 pp.
  * [29] Y. Yang, Solitons in Field Theory and Nonlinear Analysis, Springer-Verlag, 2001.
  * [30] L. Zhang, Blowup solutions of some nonlinear elliptic equations involving exponential nonlinearities. Comm. Math. Phys. 268 (2006), no. 1, 105–133.
  * [31] L. Zhang, Asymptotic behavior of blowup solutions for elliptic equations with exponential nonlinearity and singular data. Commun. Contemp. Math.  11 (2009), no. 3, 395–411.
  * [32] L. Zhang, Vanishing estimates for fully bubbling solutions of $SU(n+1)$ Toda Systems at a singular source, Int. Math. Res. Not.  2020, no. 18, 5774-–5795.