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

# On Non-simple blowup solutions of Liouville equation

Teresa D’Aprile  Dipartimento di Matematica, Università di Roma “Tor Vergata”,
via della Ricerca Scientifica 1, 00133 Roma . Italy<EMAIL_ADDRESS>,  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.

For Liouville equation with quantized singular sources, the non-simple blowup
phenomenon has been a major difficulty for years. It was conjectured by the
first two authors that the non-simple blowup phenomenon does not occur if the
equation is defined on the unit ball with Dirichlet boundary condition. In
this article we not only completely settle this conjecture in its entirety,
but also extend our result to cover any bounded domain. Since the main theorem
in this article rules out the non-simple phenomenon in commonly observed
applications, it may pave the way for advances in degree counting programs,
uniqueness of blowup solutions and construction of solutions, etc.

###### Key words and phrases:

Non-simple blowup, quantized singularity, Liouville equation, blowup
solutions, spherical Harnack inequality, Phozaev identity, Fourier Analysis

###### 1991 Mathematics Subject Classification:

35J75,35J61

##  1\. Introduction

In this article we study the following Liouville equation with finite singular
sources: let $\Omega$ be an open, bounded and connected subset of
$\mathbb{R}^{2}$ with smooth boundary $\partial\Omega$ and let $u$ be defined
as

(1.1) $\Delta
u+e^{u}=\sum_{t=1}^{M}4\pi\gamma_{t}\delta_{p_{t}}\quad\mbox{in}\quad\Omega$

where $p_{1},...,p_{M}$ are points in $\Omega$, $\gamma_{1},...,\gamma_{M}$
are constants greater than $-1$. If a $\gamma_{t}$ is a positive integer
($\gamma_{t}\in\mathbb{N}$, $\mathbb{N}$ is the set of natural numbers), we
say $p_{t}$ is a _quantized singular source_. Studying the asymptotic behavior
of blowup solutions is particularly important for a number of applications.
Let $\mathfrak{u}_{k}$ be a sequence of solutions of (1.1); we say
$\\{\mathfrak{u}_{k}\\}$ is a sequence of blowup solutions if, for a point
$p\in\Omega$, there exist $x_{k}\to p$ such that
$\mathfrak{u}_{k}(x_{k})-2\gamma_{p}\log|x_{k}-p|\to\infty$. Here we set
$\gamma_{p}=0$ if $p$ is not a singular source. For a sequence of blowup
solutions $\mathfrak{u}_{k}$, it is standard to assume a uniform bound on the
total integration and boundary oscillation: there exists $C>0$ independent of
$k$ such that

(1.2) $\int_{\Omega}e^{\mathfrak{u}_{k}}\leq C$

and

(1.3) $|\mathfrak{u}_{k}(x)-\mathfrak{u}_{k}(y)|\leq C,\quad\forall
x,y\in\partial\Omega.$

If a blowup point $p$ is either a regular point or a “non-quantized” singular
source, the asymptotic behavior of $\mathfrak{u}_{k}$ around $p$ is well
understood (see [3, 7, 13, 14, 21, 26, 27, 32, 45, 44]). As a matter of fact,
$\mathfrak{u}_{k}$ satisfies the spherical Harnack inequality around $p$,
which implies that, after scaling, the sequence $\mathfrak{u}_{k}$ behaves as
a single bubble around the maximum point. However, if $p$ happens to be a
quantized singular source, the so-called “non-simple” blowup phenomenon does
happen (see [25, 39, 40, 41]), which is equivalent to stating that
$\mathfrak{u}_{k}$ violates the spherical Harnack inequality around $p$. The
study of non-simple blowup solutions has been a major challenge for Liouville
equations and its research laid dormant for years. Recently significant
progress has been made by Kuo-Lin, Bartolucci-Tarantello and other authors [5,
6, 17, 25, 39, 40, 41]. In particular it is established in [6] and [25] that
there are $\gamma_{p}+1$ local maximum points and they are evenly distributed
on $\mathbb{S}^{1}$ after scaling according to their magnitude.

The case $\gamma_{p}\in\mathbb{N}$ is more difficult to treat, and at the same
time the most relevant to physical applications. Indeed, in vortex theory the
number $\gamma_{p}$ represents vortex multiplicity, so that in that context
the most interesting case is precisely when it is a positive integer. The
difference between the case $\gamma_{p}\in\mathbb{N}$ and
$\gamma_{p}\not\in\mathbb{N}$ is also analytically essential. Indeed, as usual
in problems involving concentration phenomena like (1.1), after suitable
rescaling of the blowing-up around a concentration point one sees a limiting
equation which, in this case, takes the form of the planar singular Liouville
equation:

$\Delta U+e^{U}=4\pi\gamma\delta_{p}\;\;\hbox{ in
}\mathbb{R}^{2},\quad\int_{\mathbb{R}^{2}}e^{U}dx<\infty;$

only if $\gamma_{p}\in\mathbb{N}$ the above limiting equation admits non-
radial solutions around $p$ since the family of all solutions extends to one
carrying an extra parameter (see [35]). This suggests that if
$\gamma_{p}\in\mathbb{N}$ and the blow-up point happens to be the singular
source, then solutions of (1.1) may exhibit non-simple blow-up phenomenon.

So, from analytical viewpoints the study of non-simple blowup solutions is far
more challenging than simple blowup solutions, but the impact of this study
may be even more significant because non-simple blowup solutions represent
certain situations in the blowup analysis of systems. If local maximums of
blowup solutions in a system tend to one point, the profile of solutions can
be described by a Liouville equation with quantized singular source. It is
desirable to know exactly when non-simple blowup phenomenon happens. In [17]
the first two authors studied the following equation:

(1.4) $\displaystyle\Delta u+\lambda e^{u}$
$\displaystyle=\sum_{t=1}^{M}4\pi\gamma_{t}\delta_{p_{t}}\quad\mbox{in}\quad\Omega\subset\mathbb{R}^{2},$
$\displaystyle u$ $\displaystyle=0\quad\mbox{on}\quad\partial\Omega,$

where $\Omega$ is an open and bounded subset of $\mathbb{R}^{2}$,
$p_{1},...,p_{M}\in\Omega$, $\partial\Omega$ is smooth, $\lambda>0$ and
$\gamma_{t}\in\mathbb{N}$. We note equation (1.4) transforms into a Liouville
equation with no boundary condition of the type (1.1) by using the presence of
the free parameter $\lambda>0$ under the transformation
$\mathfrak{u}=u(x)+\log\lambda$; then, we say that $u_{k}$ is a sequence of
blow-up solutions of (1.4) with parameter $\lambda_{k}\to 0^{+}$ if, for some
point $p\in\Omega$, there exists $x_{k}\to p$ such that
$u_{k}(x)+\log\lambda_{k}-2\gamma_{p}\log|x_{k}-p|\to+\infty$ and there is a
uniform bound on its mass:

(1.5) $\lambda_{k}\int_{\Omega}e^{u_{k}}<C.$

In [17] it was conjectured that if $\Omega=B_{1}$, $M=1$ and $p_{1}=0$, then
there is no non-simple blowup sequence in $B_{1}$. This conjecture of D’Aprile
and Wei was considered audacious since there was no restriction of number of
bubbles inside $\Omega$ and only the Dirichlet boundary condition is placed.
The purpose of this conjecture is to claim that the boundary data seem to have
a great influence on the profile of blowup solutions inside. In comparison,
the series of works of Wei-Zhang [39, 40, 41] focuses on the vanishing rate of
coefficient functions. In this respect, the conjecture of D’Aprile and Wei
seems more natural and useful for application. In the first main result we
completely settle this conjecture in a far more general form:

######  Theorem 1.1.

Let $u_{k}$ be a sequence of blowup solutions of (1.4) with parameter
$\lambda_{k}$ that satisfies (1.5). Then $u_{k}$ is simple around any blowup
point in $\Omega$.

Theorem 1.1 does not impose any symmetry condition on $\Omega$. As long as
$\Omega$ is a bounded open set with smooth boundary, the conclusion of Theorem
1.1 holds under a natural uniform bound on its mass
$\lambda\int_{\Omega}e^{u}$. This comes as a complete surprise. (Note that by
a nice result of Battaglia [8], the mass $\lambda\int_{\Omega}e^{u}$ is
uniformly bounded as long as $\Omega$ is simply connected and $M=1$. In
particular when $\Omega=B_{1},M=1$, the mass is uniformly bounded. As a result
we have completely solved the original conjecture of the first two authors.)

If the boundary condition is the usual oscillation finiteness assumption
(1.3), we can also rule out non-simple blowup in a surprising way. The second
main result is:

######  Theorem 1.2.

Let $\mathfrak{u}_{k}$ be a sequence of blowup solutions of (1.1) such that
(1.2) and (1.3) hold. If there are at least two blowup points in $\Omega$,
each blowup point is a simple blowup point.

Theorem 1.2 is also unexpected. First there is no requirement on the data of
$\mathfrak{u}_{k}$ on $\partial\Omega$ except that the oscillation of
$\mathfrak{u}_{k}$ is finite. There is also no specific requirement of what
$\Omega$ has to be. As long as there are at least two blowup points in
$\Omega$, each one of them has to be simple. Thus non-simple blow-up, if
happens, can only be single. This greatly simplifies the blowup analysis in
many applications.

The conclusions of Theorem 1.2 and Theorem 1.1 rule out non-simple blowup
phenomenon in several applicable situations. They seem to suggest that the
only case that non-simple blowup solutions occur is when the profile of blowup
solutions is very close to global solutions in the classification theorem of
Prajapat-Tarantello [35]. The proofs of the main results should lead to
advances in multiple related problems. Even though we study only one equation
in this article, it represents certain situations in systems. For example, in
the blowup analysis of Toda systems, which has ties with conformal geometry,
algebraic geometry, integrable system and complex analysis [28, 36, 37, 42,
46], one always needs to compare blowup speeds of different components. If
different components all tend to infinity in a neighborhood of one blowup
point, the behavior of the “fast” component is similar to a _quantized_
singular source to “slow” components. Therefore the asymptotic behavior of
blowup solutions of Liouville equation with quantized singular source provides
crucial information for systems.

It is also important to point out that the study of singular equation with
“quantized” singular source is ubiquitous in mathematical literature; the
nonsimple blow-up phenomenon also appears in the research of Liouville system
[23], prescribing Q curvature equation [1], Monge-Ampere equation [14] and the
vortices in a planar model of Euler flows [19], etc.

The proof of Theorem 1.2 comes naturally from Theorem 2.1, which contains most
of the key ideas. The proof of Theorem 2.1 is by way of contradiction. If non-
simple blowup happens, the blowup solutions (denoted $v_{k}$) would have $N+1$
local maximums evenly distributed around the unit disk. We shall use $N+1$
global solutions, each is very close to $v_{k}$ near a local maximum. Even
though these $N+1$ global solutions are close to one another, their mutual
difference can be captured as kernels of linearized Liouville equation. The
mutual locations of the local maximum points plays a key role in our argument.
One way to understand this argument is that $v_{k}$ cannot be very close to
different global solutions at the same time. In order to obtain a
contradiction, the point-wise estimate has to be very precise. In this article
we use Fourier analysis to obtain precise pointwise estimates around one “easy
point” first. Then this accuracy can be passed to other regions by Harnack
inequality. In the final step a contradiction can be obtained by comparing
Pohozaev identities of $v_{k}$ and other global solutions. Another key
ingredient in our proof is that the behavior of $v_{k}$ on the boundary of its
domain is significantly different from all the approximating global solutions.
This difference can be turned into a contrast on coefficient functions, which
leads to a contradiction from Pohozaev identities. This set of ideas has
turned out to be not only successful for single equations [39, 40, 41], but
also in Toda systems [38]. Many related projects such as construction of
blowup solutions, the uniqueness of blowup sequence, etc will also be
influenced by the proof in this article. The research in these directions will
be carried out soon.

The organization of the article is as follows. From Chapter 2 to Chapter 4 we
present the key proposition needed for the proof of the main theorem. The
proof of Theorem 2.1 in these chapters contains the key ingredients, which are
inspired by a series of works of the second and third authors [39, 40, 41].
Finally we place some computations in the appendix, which consists of the last
three chapters.

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\. Local Approximation

In the section we study the following local equation defined in the unit disk
$B_{1}$ in $\mathbb{R}^{2}$:

$\Delta\mathfrak{u}+e^{\mathfrak{u}}=4\pi N\delta_{0},\quad\mbox{in}\quad
B_{1}$

where $N$ is a positive integer. Since
$\Delta(\frac{1}{2\pi}\log|x|)=\delta_{0}$ we can use this function to write
the equation above as

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

if we set $u(x)=\mathfrak{u}(x)-2N\log|x|$. The purpose of this section is to
study blowup solutions of (2.1).

Let $u_{k}$ be a sequence of solutions of (2.1):

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

We say $u_{k}$ is a sequence of blowup solutions with blowup point at the
origin, if there exists $x_{k}\to 0$ such that $u_{k}(x_{k})\to\infty$ as
$k\to\infty$. Suppose the oscillation of $u_{k}$ on the boundary of $B_{1}$ is
finite:

(2.3) $|u_{k}(x)-u_{k}(y)|\leq C,\quad\forall x,y\in\partial B_{1}$

for some $C>0$ independent of $k$, and there is a uniform bound on the
integration of $|x|^{2N}e^{u_{k}}$:

(2.4) $\int_{B_{1}}|x|^{2N}e^{u_{k}}<C.$

Our goal is to study the asymptotic behavior of $u_{k}$ near the origin and
its relation with the oscillation of $u_{k}$ on $\partial B_{1}$ (the boundary
of $B_{1}$). For this purpose we set

$\Phi_{k}(x)=u_{k}(x)-\frac{1}{2\pi}\int_{\partial B_{1}}u_{k},\quad x\in
B_{1}.$

Since $u_{k}$ has bounded oscillation on $\partial B_{1}$, $\Phi_{k}(0)=0$ and
all the derivatives of $\Phi_{k}$ are uniformly bounded in $B_{1/2}$. Let
$\Phi$ be the limit of $\Phi_{k}$ over any fixed compact subset of $B_{1}$.

Then our assumption of $\Phi_{k}$ is

(2.5) $\mbox{Either }\Phi\neq 0\quad\mbox{or}\quad\Phi_{k}\equiv 0.$

######  Theorem 2.1.

Let $u_{k}$ be a sequence of blowup solutions of (2.2) that takes $0$ as its
only blowup point in $B_{1}$. Suppose (2.5) and (2.4) hold. Then $u_{k}$ is a
simple blowup sequence:

$u_{k}(x)+2(1+N)\log|x|\leq C$

for some $C>0$.

##  3\. Proof of Theorem 2.1

Suppose non-simple blowup does happen. It is well known [25, 5] that there are
exactly $N+1$ local maximums forming a circle around the origin. We use
$p_{0}^{k}$,…,$p_{N}^{k}$ 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 $\bar{\mu}_{k}$ as follows:

(3.1) $\delta_{k}=|p_{0}^{k}|\quad\mbox{and
}\quad\bar{\mu}_{k}=u_{k}(p_{0}^{k})+2(1+N)\log\delta_{k}.$

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,
$\bar{\mu}_{k}\to\infty$. The interested readers may look into [25, 5] for
more detailed information.

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

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

If we write out each component, (3.2) 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.3) $\Delta v_{k}(y)+|y|^{2N}e^{v_{k}(y)}=0,\quad|y|<\delta_{k}^{-1},$

Thus the image of $p_{0}^{k}$ after scaling is $Q_{0}^{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=0,...,N.$

It is established by Kuo-Lin in [25] and independently by Bartolucci-
Tarantello in [5] 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 [39])

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

Choosing $3\tau>0$ small and independent of $k$, we can make disks centered at
$Q_{l}^{k}$ with radius $3\tau$ (denoted as $B(Q_{l}^{k},3\tau)$) mutually
disjoint. The $\bar{\mu}_{k}$ in (3.1) is

$\bar{\mu}_{k}=\max_{B(Q_{0}^{k},\tau)}v_{k}.$

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

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

Let

(3.6)
$V_{k}(x)=\log\frac{e^{\bar{\mu}_{k}}}{(1+\frac{e^{\bar{\mu}_{k}}}{8(1+N)^{2}}|y^{N+1}-e_{1}|^{2})^{2}}.$

Clearly $V_{k}$ is a solution of

(3.7) $\Delta
V_{k}+|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
[35]. Now we use the following expansion of $V_{k}$ for $|y|=L_{k}$
($L_{k}=\delta_{k}^{-1}$)

(3.8) $\displaystyle V_{k}(y)$
$\displaystyle=-\bar{\mu}_{k}+2\log(8(N+1)^{2})-4(N+1)\log
L_{k}+\frac{2}{L_{k}^{2N+2}}$
$\displaystyle+\frac{4\cos((N+1)\theta)}{L_{k}^{N+1}}+\frac{4}{L_{k}^{2N+2}}\cos((2N+2)\theta)$
$\displaystyle+O(L_{k}^{-3N-3})+O(e^{-\bar{\mu}_{k}}L_{k}^{-2N-2}).$

To eliminate the main oscillation of $v_{k}-V_{k}$ on $\partial\Omega_{k}$, we
set $\phi_{v,k}(\delta_{k}\cdot)$ be the harmonic function on
$\partial\Omega_{k}$ such that $V_{k}(y)-\phi_{v,k}(\delta_{k}y)$ has no
oscillation. The reason we use $\phi_{v,k}(\delta_{k}\cdot)$ is because
$\phi_{v,k}$ is the harmonic function that eliminates the oscillation of
$\tilde{V}_{k}$ on $\partial B_{1}$. $\tilde{V}_{k}$ is a re-scaled $V_{k}$:

$\tilde{V}_{k}(\delta_{k}y)+2(1+N)\log\delta_{k}=V_{k}(y).$

on $\partial B_{1}$. From the expression of $V_{k}$ we see that $\phi_{v,k}\to
0$ uniformly on $B_{1}$. It is easy to see that the leading terms of
$\phi_{v,k}(\delta_{k}y)$ are

$\phi_{v,k}(\delta_{k}y)=\frac{4}{L_{k}^{2N+2}}r^{N+1}\cos((N+1)\theta)+\frac{4}{L_{k}^{4N+4}}r^{2N+2}\cos((2N+2)\theta)+....$

Naturally we define $\phi_{0,k}(y)$, $h_{0,k}$ and $v_{0,k}$ as

(3.9)
$\displaystyle\phi_{0,k}(y)=\Phi_{k}(\delta_{k}y)-\phi_{v,k}(\delta_{k}y)$
$\displaystyle=\Phi_{k}(\delta_{k}y)-4\delta_{k}^{2N+2}r^{N+1}\cos((N+1)\theta)-\delta_{k}^{4N+4}r^{2N+2}\cos((2N+2)\theta)+...$
(3.10) $h_{0,k}=e^{\phi_{0,k}}.$ $v_{0,k}=v_{k}-\phi_{0,k},$

and we write the equation of $v_{0,k}$ as

(3.11) $\Delta
v_{0,k}+h_{0,k}|y|^{2N}e^{v_{0,k}}=0,\quad\mbox{in}\quad\Omega_{k}$

Based on the definition of $h_{0,k}$ in (3.10) we prove

######  Lemma 3.1.

There exist and integer $L>0$ ,$\delta_{k}^{*}\in(\delta_{k}^{L},\delta_{k})$
an integer $0\leq s\leq N$ such that

(3.12) $\nabla h_{0,k}(e^{\frac{2\pi is}{N+1}})/\delta_{k}^{*}\neq 0.$

Proof of Lemma 3.1: The proof of (3.12) is based on (2.5). We first write the
expansion of $\Phi_{k}(x)$ as

$\Phi_{k}(x)=\sum_{n=1}^{\infty}r^{n}(a_{n}^{k}\cos(n\theta)+b_{n}^{k}\sin(n\theta))$

where $x=re^{i\theta}$ and the two sequences $|a_{n}^{k}|$ and $|b_{n}^{k}|$
are uniformly bounded from above. They could tend to zero for some $n$. The
expansion of $\phi_{v,k}(x)$ can be written as

$\phi_{v,k}(x)=\sum_{n=1}^{\infty}r^{n}(a_{n,v}^{k}\cos(n\theta)+b_{n,v}^{k}\sin(n\theta))$

with $a_{n,v}^{k}$ and $b_{n,v}^{k}$ both tending to zero as $k\to\infty$ for
each fixed $n$. Then clearly

$\displaystyle\phi_{0,k}(y)$
$\displaystyle=\Phi_{k}(\delta_{k}y)-\phi_{v,k}(\delta_{k}y)$
$\displaystyle=\sum_{n=1}^{\infty}\delta_{k}^{n}r^{n}\bigg{(}(a_{n}^{k}-a_{n,v}^{k})\cos
n\theta+(b_{n}^{k}-b_{n,v}^{k})\sin n\theta)\bigg{)}.$

Obviously we use $y=re^{i\theta}$ in this case. If $\Phi\neq 0$ (see 2.5)
there exists an integer $L$ such that $|a_{L}^{k}|+|b_{L}^{k}|\geq C>0$ for
some $C>0$ independent of $k$. Now we set $\delta_{k}^{*}$ as

$\delta_{k}^{*}:=\sum_{n=1}^{L}\delta_{k}^{n}(|a_{n}^{k}-a_{n,v}^{k}|+|b_{n}^{k}-b_{n,v}^{k}|).$

We make three trivial remarks about $\delta_{k}^{*}$. First

$\delta_{k}^{*}\geq C\delta_{k}^{L}$

because $|a_{L}^{n}|+|b_{L}^{n}|\geq C$ for some $C$ independent of $k$.
Second, since it is only a finite sum in the definition of $\delta_{k}^{*}$,
$\delta_{k}^{*}$ is comparable to

$\max_{1\leq n\leq
L}(|a_{n}^{k}-a_{n,v}^{k}|+|b_{n}^{k}-b_{n,v}^{k}|)\delta_{k}^{n}.$

The third observation is

$\sum_{n=L+1}^{\infty}(|a_{n}^{k}-a_{n,v}^{k}|+|b_{n}^{k}-b_{n,v}^{k}|)\delta_{k}^{n}=O(\delta_{k}^{L+1}).$

In other words, the quantity is far less than $\delta_{k}^{*}$. Our goal is to
compute $\nabla\phi_{0,k}$ at points close to $\partial B_{1}$. We use the
following formula:

$|\nabla\phi_{0,k}(x)|^{2}=|\partial_{r}\phi_{0,k}|^{2}+\frac{1}{r^{2}}|\partial_{\theta}\phi_{0,k}|^{2},\quad|x|\sim
1.$

For convenience we use

$(a_{n}^{k}-a_{n,v}^{k})\delta_{k}^{n}=\delta_{k}^{*}|c_{n}^{k}|\cos(\beta_{n}^{k}),\quad(b_{n}-b_{n,v})\delta_{k}^{n}=\delta_{k}^{*}|c_{n}^{k}|\sin(\beta_{n}^{k}).$

Note that at least one $|c_{n}^{k}|\sim 1$. In fact we consider those
comparable to $1$ in the summation from $1$ to $L$. Then

$\partial_{r}\phi_{0,k}=\delta_{k}^{*}\sum_{n=1}^{L}nr^{n-1}|c_{n}^{k}|\cos(\beta_{n}^{k}-n\theta)+O(\delta_{k}^{L+1})$
$\displaystyle\partial_{\theta}\phi_{0,k}$ $\displaystyle=$
$\displaystyle\delta_{k}^{*}\sum_{n=1}^{L}r^{n}(-n|c_{n}^{k}|\cos(\beta_{n}^{k})\sin(n\theta)+n|c_{n}^{k}|\sin(\beta_{n}^{k})\cos(n\theta))+O(\delta_{k}^{L+1})$
$\displaystyle=$
$\displaystyle\delta_{k}^{*}\sum_{n=1}^{N}nr^{n}|c_{n}^{k}|\sin(\beta_{n}^{k}-n\theta)+O(\delta_{k}^{L+1}).$
(3.13)
$\displaystyle|\nabla\phi_{0,k}|^{2}=|\partial_{r}\phi_{0,k}|^{2}+\frac{1}{r^{2}}|\partial_{\theta}\phi_{0,k}|^{2}$
$\displaystyle=$
$\displaystyle(\delta_{k}^{*})^{2}(\sum_{n=1}^{L}n^{2}r^{2n-2}|c_{n}^{k}|^{2}+2\sum_{s<t}str^{s+t-2}|c_{s}^{k}||c_{t}^{k}|\cos(\beta_{s}^{k}-\beta_{t}^{k}-(s-t)\theta)$
$\displaystyle+O(\delta_{k}^{L+1}).$

Now we have the freedom to choose $e^{\frac{2\pi li}{N+1}}$. Since all these
points are on the unit circle, $r=1$ in the evaluation of (3.13). By Cauchy’s
inequality we see that if $|\nabla\phi_{0,k}|=O(\delta_{k}^{L+1})$ at one of
$e^{\frac{2\pi il}{N+1}}$ it is obviously not the case on another such point,
which means at that point, say $e^{\frac{2\pi il}{N+1}}$,
$|\nabla\phi_{0,k}(e^{\frac{2\pi il}{N+1}})|$ is comparable to
$\delta_{k}^{*}$. Lemma 3.1 is established if $\Phi\neq 0$.

Finally if $\Phi_{k}\equiv 0$, based on the expansion of $V_{k}$ we take
$\delta_{k}^{*}=\delta_{k}^{2N+2}$. Lemma 3.1 is established. $\Box$

If we let $Q_{0}^{k}$ denote the local maximum point of $v_{0,k}$, the
difference between $Q_{0}^{k}$ and $e_{1}$ is
$O(\delta_{k}^{*}e^{-\bar{\mu}_{k}})$ by the non-degeneracy of $v_{k}$ around
$e_{1}$ (Indeed, if we use $\bar{\epsilon}_{k}=e^{-\bar{\mu}_{k}/2}$ to be the
scaling factor, after scaling, it is easy to see that the location of the
critical point is $O(\bar{\epsilon}_{k}\delta_{k}^{*})$ for the scaled
function. In the original setting before the scaling, the location of the
local maximum is $O(\delta_{k}^{*}e^{-\bar{\mu}_{k}})$ away from $e_{1}$. The
maximum of $v_{0,k}$ is $\bar{\mu}_{k}+O(\delta_{k}^{*})$.

Let $V_{0,k}$ be the global solution of that agrees with $v_{0,k}$ at
$Q_{0}^{k}$:

(3.14) $v_{0,k}(Q_{0}^{k})=V_{0,k}(Q_{0}^{k}),\quad\nabla
v_{0,k}(Q_{0}^{k})=\nabla V_{0,k}(Q_{0}^{k})=0$

and the equation for $V_{0,k}$ is

(3.15) $\Delta
V_{0,k}+|y|^{2N}h_{0,k}(Q_{0}^{k})e^{V_{0,k}}=0,\quad\mbox{in}\quad\mathbb{R}^{2},\quad\int_{\mathbb{R}^{2}}|y|^{2N}e^{V_{0,k}}<\infty.$

Of course $V_{0,k}$ is a small perturbation of the previous $V_{k}$, and we
see that the oscillation of $v_{0,k}-V_{0,k}$ is
$O(\delta_{k}^{*}\delta_{k}^{N+1}e^{-\bar{\mu}_{k}})$ on $\partial\Omega_{k}$.
At this moment we set the scaling factor to be

$\epsilon_{k}=e^{-\frac{1}{2}\mu_{k}}.\quad\mu_{k}=v_{0,k}(Q_{0}^{k})$

Indeed, the expression of $V_{0,k}$ is

$V_{0,k}(y)=\log\frac{e^{\mu_{k}}}{(1+\frac{h_{0,k}(Q_{0}^{k})e^{\mu_{k}}}{D}|y^{N+1}-1-p_{k}|^{2})^{2}}$

where $p_{k}=O(e^{-\mu_{k}}\delta_{k}^{*})$. From here we see that the
perturbation of $v_{0,k}-V_{0,k}$ on $\partial\Omega_{k}$ is
$O(\delta_{k}^{*}\delta_{k}^{N+1}e^{-\mu_{k}})$:

(3.16) $|(v_{0,k}-V_{0,k})(x)-(v_{0,k}-V_{0,k})(y)|\leq
C\delta_{k}^{*}\delta_{k}^{N+1}e^{-\mu_{k}},\quad\forall
x,y\in\partial\Omega_{k}.$

###  3.1. Point-wise estimate for $v_{0,k}-V_{0,k}$

In the appendix we establish the closeness of local maximum points of $v_{k}$
with $e^{\frac{2i\pi l}{N+1}}$ for $l=0,...,N$. As mentioned before the local
maximums of $v_{0,k}$ is only $O(\delta_{k}^{*}\epsilon_{k}^{2})$ perturbation
of the corresponding local maximum points of $v_{k}$.

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.

Now we cite Proposition 3.1 of [40]:

_Proposition 3.1 of[40]: 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)$_

(3.17)
$|v_{0,k}(x)-V_{0,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 3.1.

We only need a re-scaled version of Proposition 3.1 of [40]:

(3.18)
$|v_{0,k}(e^{i\beta_{l}}+\epsilon_{k}y)-V_{0,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$.

######  Remark 3.2.

The main idea of the proof of Theorem 2.1 can be observed from the equations
of $v_{0,k}$ (3.11) and $V_{0,k}$ (3.15). While (3.15) has a constant
coefficient $h_{0,k}(Q_{0}^{k})$, (3.11) has a function $h_{0,k}$ whose
derivative is not zero at some $Q_{s}^{k}$. We shall obtain a precise
pointwise estimate of $v_{0,k}-V_{0,k}$ and the difference on these
coefficient functions will lead to a contradiction.

######  Proposition 3.1.

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

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

Proof of Proposition 3.1: Let $M_{k}$ be the maximum of $|w_{0,k}|$ on
$\Omega_{k}$. By way of contradiction we assume that

(3.19) $M_{k}/(\delta_{k}^{*})\to\infty.$

Now we recall the equation for $v_{0,k}$ is (3.3), $v_{0,k}-V_{0,k}$ is very
close to a constant on $\partial B(0,\delta_{k}^{-1})$. Moreover

(3.20) $w_{0,k}(Q_{0}^{k})=|\nabla w_{0,k}(Q_{0}^{k})|=0.$

Recall that $V_{0,k}$ defined in (3.6) satisfies (3.15).

We shall derive a precise, point-wise estimate of $w_{0,k}$ in
$B_{3}\setminus\cup_{l=1}^{N}B(Q_{l}^{k},\tau)$ where $\tau>0$ is a small
number independent of $k$. We shall prove that $w_{0,k}$ is very small in
$B_{3}$ if we exclude all bubbling disks except the one around $e_{1}$.

Now based on (3.11) and (3.15) we write the equation of $w_{0,k}$ as

(3.21) $\Delta
w_{0,k}+|y|^{2N}h_{0,k}(y)e^{\xi_{k}}w_{0,k}=|y|^{2N}(h_{0,k}(Q_{0}^{k})-h_{0,k}(y))e^{V_{0,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_{0,k}(x)}-e^{V_{0,k}(x)}}{v_{0,k}(x)-V_{0,k}(x)},\quad\mbox{if}\quad
v_{0,k}(x)\neq V_{0,k}(x),\\\ \\\ e^{V_{0,k}(x)},\quad\mbox{if}\quad
v_{0,k}(x)=V_{0,k}(x).\end{array}\right.$

An equivalent form is

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

Note that the oscillation of $w_{0,k}$ on $\partial\Omega_{k}$ is
$O(\delta_{k}^{*}\delta_{k}^{N+1}e^{-\mu_{k}})$. A trivial observation is that
the right hand side of (3.21) is zero if $y=Q_{0}^{k}$. This simple fact
determines a pathway for our argument. We will prove smallness of $w_{k}$
around $Q_{0}^{k}$ first and then pass it to other places.

By normalizing $w_{0,k}$ we shall study the following function:

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

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

(3.23)
$\Delta\tilde{w}_{k}(y)+|y|^{2N}h_{0,k}(y)e^{\xi_{k}}\tilde{w}_{k}(y)=O(\sigma_{k})|y|^{2N}(y-Q_{0}^{k})e^{V_{0,k}}+O(\sigma_{k})|y-Q_{0}^{k}|^{2}e^{V_{0,k}}$

in $\Omega_{k}$. Here $\delta_{k}^{*}/M_{k}=\sigma_{k}\to 0$. Also on the
boundary, the oscillation of $\tilde{w}_{k}$ is
$o(\delta_{k}^{N+1}e^{-\mu_{k}})$. By Proposition 3.1 of [40]

(3.24)
$\xi_{k}(Q_{0}^{k}+\epsilon_{k}z)=V_{0,k}(Q_{0}^{k}+\epsilon_{k}z)+O(\epsilon_{k}^{\epsilon})(1+|z|)^{-1}$

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

(3.25) $\Delta
U_{k}+h_{0,k}(Q_{0}^{k})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 [9] we have

$U_{k}(z)=\log\frac{1}{(1+\frac{h_{0,k}(Q_{0}^{k})}{8}|z|^{2})^{2}}$

and standard refined estimates yield (see [15, 45, 21])

(3.26)
$V_{0,k}(Q_{0}^{k}+\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

(3.27) $\log|Q_{0}^{k}+\epsilon_{k}z|=O(\epsilon_{k})(1+|z|).$

Thus, the combination of (3.24), (3.26) and (3.27) gives

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

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

(3.29)
$\epsilon_{k}^{2}|Q_{0}^{k}+\epsilon_{k}z|^{2N}h_{0,k}(Q_{0}^{k}+\epsilon_{k}z)e^{\xi_{k}(Q_{0}^{k}+\epsilon_{k}z)}=h_{0,k}(Q_{0}^{k})e^{U_{k}(z)}+O(\epsilon_{k}^{\epsilon})(1+|z|)^{-3}$

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

######  Lemma 3.2.

(3.30) $|\tilde{w}_{k}(Q_{0}^{k}+\epsilon_{k}z)|\leq
o(\epsilon_{k})(1+|z|),\quad 0<|z|<\tau\epsilon_{k}^{-1}.$

for some $\tau>0$.

Proof of Lemma 3.2: Step one: In this step we prove the following statement:
For $\delta>0$ small and independent of $k$,

(3.31) $\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.

If (3.31) is not true, we have, without loss of generality that
$\tilde{w}_{k}\to c\neq 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. Let

(3.32)
$W_{k}(z)=\tilde{w}_{k}(Q_{0}^{k}+\epsilon_{k}z),\quad\epsilon_{k}=e^{-\frac{1}{2}v_{0,k}(Q_{0}^{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}$, (3.26) and (3.29) we write the
equation of $W_{k}$ as

(3.33) $\Delta W_{k}(z)+h_{0,k}(Q_{0}^{k})e^{U_{k}(z)}W_{k}(z)=E_{k},$

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

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

Let

(3.34) $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)+h_{0,k}(Q_{0}^{k})e^{U_{k}(r)}g_{0}^{k}(r)=E_{0}^{k}(r)$
$\displaystyle g_{0}^{k}(0)=\frac{d}{dr}g_{0}^{k}(0)=0.$

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

$|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}^{k}$, $g_{02}^{k}$, where

$g_{01}^{k}=\frac{1-c_{k}r^{2}}{1+c_{k}r^{2}},\quad
c_{k}=\frac{h_{0,k}(Q_{0}^{k})}{8}.$

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

$\displaystyle|g_{0}^{k}(r)|\leq
C|g_{01}^{k}(r)|\int_{0}^{r}s|E_{0}^{k}(s)g_{02}^{k}(s)|ds+C|g_{02}^{k}(r)|\int_{0}^{r}s|g_{01}^{k}(s)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 (3.34). 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. (3.31) is
established.

Step two: Now we extend the estimate to the whole neighborhood of $e_{1}$. The
estimate will be obtained in a progressive way. Let $W_{k}$ be defined as in
(3.32). In order to obtain a better estimate we need to write the equation of
$W_{k}$ more precisely than (3.33):

(3.35) $\Delta
W_{k}+h_{0,k}(Q_{0}^{k})e^{\Theta_{k}}W_{k}=o(\sigma_{k})\epsilon_{k}y(1+|y|)^{-4},\quad
z\in B(0,\tau\epsilon_{k}^{-1})$

where $\Theta_{k}$ is defined by

$e^{\Theta_{k}(z)}=|Q_{0}^{k}+\epsilon_{k}z|^{2N}e^{\xi_{k}(Q_{0}^{k}+\epsilon_{k}z)+2\log\epsilon_{k}}.$

Here we observe that by step one $W_{k}=o(1)$ on $\partial
B(0,\tau\epsilon_{k}^{-1})$. In the computation of (3.35) we also used

$h_{0,k}(Q_{0}^{k}+\epsilon_{k}y)=h_{0,k}(Q_{0}^{k})+O(\epsilon_{k}y).$

To replace $\Theta_{k}$ by $U_{k}$ we have an extra error that depends on the
bound of $W_{k}$:

(3.36) $\Delta W_{k}+h_{0,k}(Q_{0}^{k})e^{U_{k}}W_{k}=E_{k},\quad z\in
B(0,\tau\epsilon_{k}^{-1})$

where $E_{k}=E_{1}+E_{2}$ and $E_{1}$ is the right hand side of (3.35). Thus

$|E_{1}|\leq C\sigma_{k}\epsilon_{k}y(1+|y|)^{-4}\quad|E_{2}|\leq
C\sigma_{k}\epsilon_{k}^{\epsilon}(1+|z|)^{-3}.$

Note that in the crude bound of $E_{2}$, we used $|W_{k}|\leq C$. If the bound
of $|W_{k}|$ is better, the estimate of $E_{2}$ will improve as well.

The proof is by considering the projection of $W_{k}$ in different nodes. Let
$g_{0}^{k}$ be the projection on $1$, then we have

$\displaystyle\frac{d^{2}}{dr^{2}}g_{0}^{k}(r)+\frac{1}{r}\frac{d}{dr}g_{0}^{k}(r)+h_{0,k}(Q_{0}^{k})e^{U_{k}}g_{0}^{k}=E_{k,0}$
$\displaystyle g_{0}^{k}(0)=\frac{d}{dr}g_{0}^{k}(0)=0$

where $E_{k,0}$ is the projection of $E_{k}$ to $1$. Just like in the proof of
step one, the estimate of $g_{0}^{k}$ is

$\displaystyle|g_{0}^{k}(r)|$ $\displaystyle\leq
C|g_{01}(r)|\int_{0}^{r}s|E_{k,0}(s)||g_{02}^{k}(s)|ds$
$\displaystyle+C|g_{02}^{k}(r)|\int_{0}^{r}s|E_{k,0}(s)||g_{01}^{k}(s)|ds$

where

$g_{01}^{k}(r)=\frac{1-c_{k}r^{2}}{1+c_{k}r^{2}},\quad
c_{k}=\frac{h_{0k}(Q_{0}^{k})}{8}.$

and $g_{02}^{k}=O(\log r)$ for $r>1$ and for $r$ close to $0$. If we just use

$|E_{k,0}(s)|\leq o(1)\epsilon_{k}^{\epsilon}(1+s)^{-3},$

we have

$|g_{0}^{k}(r)|=o(1)\epsilon_{k}^{\epsilon}\log(2+r).$

Let $g_{1}^{k}$ be the radial part of the projection of $W_{k}$ on
$e^{i\theta}$. Then $g_{1}^{k}$ satisfies

$\frac{d^{2}}{dr^{2}}g_{1}^{k}(r)+\frac{1}{r}\frac{d}{dr}g_{1}^{k}(r)+(h_{0,k}(Q_{0}^{k})e^{U_{k}}-\frac{1}{r^{2}})g_{1}^{k}(r)=E_{k,1}$

and we use the following crude upper bound for $E_{k,1}$

$E_{k,1}(r)=o(1)\epsilon_{k}^{\epsilon}(1+r)^{-3}$

and $g_{1}^{k}(0)=\frac{d}{dr}g_{1}^{k}(0)=0$. Then the two fundamental
solutions are

$g_{11}^{k}(r)=\frac{r}{1+c_{k}r^{2}},$

$g_{12}^{k}(r)$ behaves like $O(r^{-1})$ near $0$ and behaves like $O(r)$ near
infinity.

$\displaystyle|g_{1}^{k}(r)|\leq$ $\displaystyle
o(1)|g_{11}^{k}(r)|\int_{0}^{r}s\epsilon_{k}^{\epsilon}(1+s)^{-3}|g_{12}^{k}(s)|ds$
$\displaystyle+o(1)|g_{12}^{k}(r)|\int_{0}^{r}s\epsilon_{k}^{\epsilon}(1+s)^{-3}|g_{11}^{k}(s)|ds$
$\displaystyle=o(1)\epsilon_{k}^{\epsilon}(1+r)$

For $l\geq 2$, let $g_{l}^{k}$ be the radial part of the projection on
$e^{il\theta}$, then the equation is

$\frac{d}{dr^{2}}g_{l}^{k}(r)+\frac{1}{r}\frac{d}{dr}g_{l}^{k}(r)+(h_{0,k}(Q_{0}^{k})e^{U_{k}}-\frac{l^{2}}{r^{2}})g_{l}^{k}(r)=o(1)\epsilon_{k}^{\epsilon}(1+r)^{-3}$

with $g_{l}^{k}(0)=0$ and $g_{l}^{k}(L_{k})=s_{l,k}$ for some $s_{l,k}=o(1)$
as a result of step one. For the corresponding homogeneous equation, the two
fundamental solutions $g_{l1}^{k}$, $g_{l2}^{k}$ can be chosen to behave like

$g_{l1}^{k}\sim r^{l}\mbox{near $0$ and $\infty$},$ $g_{l2}^{k}\sim
r^{-l}\mbox{near $0$ and $\infty$}$

with bounds independent of $k$. Then

$\displaystyle|g_{l}^{k}(r)|\leq|c_{1,l}^{k}g_{l1}^{k}(r)|+o(\epsilon_{k}^{\epsilon})|g_{l1}^{k}(r)|\int_{r}^{\infty}\frac{s}{l}(1+s)^{-3}|g_{l2}^{k}(s)|ds$
$\displaystyle+o(\epsilon_{k}^{\epsilon})|g_{l2}^{k}(r)|\int_{0}^{r}\frac{s}{l}|g_{l1}^{k}(s)|(1+s)^{-3}ds$
$\displaystyle\leq|c_{l,1}^{k}g_{l1}^{k}(r)|+\frac{o(\epsilon_{k}^{\epsilon})}{l^{2}}.$

By $g_{l}(\tau\epsilon_{k}^{-1})=s_{l,k}$ we have $|c_{l,1}^{k}|\leq
Cs_{l,k}\epsilon_{k}^{l}$. Thus the summation of projections on all nodes is
convergent and the summation leads to

$|W_{k}(z)|\leq o(1)\epsilon_{k}^{\epsilon}(1+|z|).$

Using this new bound in the previous machinery we will obtain

$|W_{k}(y)|\leq o(1)\epsilon_{k}^{2\epsilon}(1+|z|).$

A repetitive application of this process eventually makes $E_{1}$ the leading
error term in $E_{k}$ and we have

$|W_{k}(z)|\leq o(1)\epsilon_{k}(1+|z|),\quad|z|\leq\tau\epsilon_{k}^{-1}.$

Lemma 3.2 is established. $\Box$

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

######  Lemma 3.3.

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

Proof of Lemma 3.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
B(0,\tau\epsilon_{k}^{-1}).$

Here we point out that based on (3.5) 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}h_{0,k}(e^{i\beta_{l}}+\epsilon_{k}z)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.$

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 consists of proving $c_{0}=0$ and
$c_{1}=c_{2}=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}=1+O(\epsilon_{k}z),$

and

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

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

(3.38) $\Delta W_{k}(z)+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

(3.39) $\Delta\phi_{0}^{k}+e^{U_{k}}\phi_{0}^{k}=0,\quad\mbox{in
}\quad\mathbb{R}^{2}.$

When $k\to\infty$, $\phi_{0}^{k}\to\phi_{0}$. From multiplying $\phi_{0}^{k}$
to (3.38) and multiplying $W_{k}$ to (3.39) we have

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

Thus from (3.40) we have

(3.41) $\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}\Big{(}\frac{o(\epsilon_{k}^{\epsilon})}{s}+O(s^{-3})\Big{)}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 3.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_{0,k}k$ is

$v_{0,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<\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,\delta_{k}^{-1})}$. A contradiction to $\max|\tilde{w}_{k}|=1$.

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

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

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

(3.42) $D=8(N+1)^{2}.$

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

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

Note that there is no need to define a $\phi_{l,k}$ and
$v_{l,k}=v_{k}-\phi_{l,k}$ because the difference is insignificant.

Since $h_{0,k}(y)=e^{\phi_{1}}$ for some harmonic function
$\phi_{1}(y)=O(\delta_{k}^{*})$ for $|y|\sim 1$, we have

(3.43) $h_{0,k}(Q_{l}^{k})-h_{0,k}(Q_{s}^{k})=O(\delta_{k}^{*})=o(1)M_{k}.$

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

(3.44)
$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_{0,k}-V_{l,k}|$ and we claim that all
these $M_{l,k}$ are comparable:

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

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

$\frac{(v_{0,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 3.2 we have

(3.46)
$v_{0,k}(Q_{s}^{k}+\epsilon_{k}z)-V_{s,k}(Q_{s}^{k}+\epsilon_{k}z)=o(\epsilon_{k})(1+|z|)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_{0,k}(Q_{s}^{k}+\epsilon_{k}z)-V_{l,k}(Q_{s}^{k}+\epsilon_{k}z)}{M_{k}}$
$\displaystyle=$ $\displaystyle
o(\epsilon_{k})(1+|z|)+\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}|y^{N+1}-e_{1}-p_{l}^{k}|^{2}-\frac{e^{\mu_{s}^{k}}}{D}|y^{N+1}-e_{1}-p_{s}^{k}|^{2}}{1+\frac{e^{\mu_{s}^{k}}}{D}|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. The difference between
$h_{l,k}(Q_{l}^{k})$ and $h_{l,k}(Q_{s}^{k})$ is in (3.43). From $A$ we claim
that

(3.47) $\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

(3.48)
$\displaystyle\phi_{1}=(\mu_{s}^{k}-\mu_{l}^{k})(1-\frac{h_{s}}{8}|z+\frac{N}{2}\epsilon_{k}z^{2}e^{-i\beta_{s}}|^{2})/B,$
$\displaystyle\phi_{2}=\frac{h_{s}}{2B}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\phi_{3}=\frac{|p_{s}^{k}-p_{l}^{k}|^{2}}{4(N+1)^{2}\epsilon_{k}^{2}B}\bigg{(}1-\frac{|z|^{2}(1+2\cos(2\theta-2\theta_{st}-2\beta_{s})}{8B}\bigg{)},$
$\displaystyle\phi_{4}=\frac{h_{s}}{4}\frac{h_{l}-h_{s}}{h_{l}}|z+O(\epsilon_{k})|z|^{2}|^{2}/B.$
$\displaystyle
B=1+\frac{h_{s}}{8}|z+\frac{N}{2}z^{2}e^{-i\beta_{s}}\epsilon_{k}|^{2},$

where $h_{s}=h_{0,k}(Q_{s}^{k})$, $h_{l}=h_{0,k}(Q_{l}^{k})$,
$\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_{3}+\phi_{4}$ is the second order term
which will play a leading role later. The derivation of (3.47) is put in the
appendix.

Here $\phi_{1}$, $\phi_{2}$ 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}(1+o(1))$. 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_{2}$ of (3.47)) we have

(3.49) $\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

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

Equation (3.49) 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}$ (3.45) 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}$. (3.45) is established. From now on for convenience we shall just use
$M_{k}$, which has an upper bound in (3.50).

Set $w_{l,k}=(v_{0,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

(3.51)
$\Delta\tilde{w}_{l,k}+|y|^{2N}h_{0,k}(Q_{l}^{k})e^{\xi_{l}^{k}}\tilde{w}_{l,k}=o(\sigma_{k})(y-Q_{l}^{k})|y|^{2N}e^{V_{l,k}}+o(\sigma_{k})|y-Q_{l}^{k}|^{2}e^{V_{l,k}},$

where $\xi_{l}^{k}$ comes from the Mean Value Theorem and satisfies

(3.52)
$e^{\xi_{l}^{k}}=e^{V_{l,k}}(1+\frac{1}{2}w_{l,k}+O(w_{l,k}^{2}))=e^{V_{l,k}}(1+\frac{1}{2}M_{k}\tilde{w}_{l,k}+O(M_{k}^{2})).$

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

(3.53)
$\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}$ (3.46) holds with $M_{s,k}$ replaced by $M_{k}$.

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

(3.54) $\displaystyle\tilde{w}_{l,k}(y)$
$\displaystyle=\int_{\Omega_{k}}H_{y,l}(\eta)|\eta|^{2N}h_{0,k}(\eta)e^{\xi_{l}}\tilde{w}_{l,k}(\eta)d\eta$
$\displaystyle+o(1)\delta_{k}^{N+1}e^{-\mu_{k}}.$

where the last term is based on (3.16) and

$H_{y,l}(\eta):=G_{k}(y,\eta)-G_{k}(Q_{l}^{k},\eta).$

If we only consider $|y|\sim 1$, we use

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

To evaluate the right hand side of (3.54), we only need to concentrate on
$B(Q_{s}^{k},\tau)$ for $s\neq l$ because the integration around $Q_{l}^{k}$
and outside bubbling disks only contribute $o(\epsilon_{k})$. Around each
$Q_{s}^{k}$ the $e^{\xi_{l}}$ can be replaced by $e^{V_{s,k}}$ with
controllable error because from Lemma 3.2, the difference between them only
leads to $o(\epsilon_{k})$ in error. Also, $h_{0,k}(\eta)$ can be replaced by
$h_{0,k}(Q_{s}^{k})$ because

$h_{0,k}(\eta)=h_{0,k}(Q_{s}^{k})+\nabla h_{0,k}(Q_{s}^{k})(\eta-
Q_{s}^{k})+O(\delta_{k}^{*})(\eta-Q_{s}^{k})^{2}.$

Since $\nabla h_{0,k}(Q_{s}^{k})=O(\delta_{k}^{*})$, we see that the
difference between $h_{0,k}(\eta)$ and $h_{0,k}(Q_{s}^{k})$ only leads to
$o(\epsilon_{k})$ in its corresponding integration.

From the decomposition in (3.47) we can now estimate the integral of
$\tilde{w}_{l,k}$ more precisely. Clearly we only need to evaluate integrals
around each $Q_{s}^{k}$. For this we have

(3.56)
$\int_{B(Q_{s}^{k},\tau)}\tilde{w}_{l,k}(\eta)h_{0,k}(Q_{l})|\eta|^{2N}e^{V_{s,k}}d\eta+o(\epsilon_{k})=D_{s,l}^{k}+o(\epsilon_{k}).$

where

$D_{s,l}^{k}=\frac{\pi}{2(N+1)^{2}}\frac{|p_{s}^{k}-p_{l}^{k}|^{2}}{\epsilon_{k}^{2}M_{k}^{2}}M_{k}+2\pi\frac{h_{l}-h_{s}}{M_{k}}$

and for convenience we skip $k$ in $Q_{l}$, and we shall write the two
components of $Q_{l}$ as $Q_{l}=(Q_{l}^{1},Q_{l}^{2})$. The derivation of
(3.56) can be found in Appendix B. Then for $|y|\sim 1$

$\displaystyle\tilde{w}_{l,k}(y)=-\sum_{s\neq
l}(H_{y,l}(Q_{s})+O(\delta_{k}^{2}))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})h_{0,k}(Q_{l})|\eta|^{2N}e^{\xi_{l}}\tilde{w}_{l,k}(\eta)\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}|}{|Q_{l}-Q_{s}|}+O(\delta_{k}^{2}))D_{s,l}^{k}$
$\displaystyle-\sum_{s\neq
l}8\Big{(}\frac{y_{1}-Q_{s}^{1}}{|y-Q_{s}|^{2}}-\frac{Q_{l}^{1}-Q_{s}^{1}}{|Q_{l}-Q_{s}|^{2}}\Big{)}c_{1,s,l}\epsilon_{k}$
$\displaystyle-8\Big{(}\frac{y_{2}-Q_{s}^{2}}{|y-Q_{s}|^{2}}-\frac{Q_{l}^{2}-Q_{s}^{2}}{|Q_{s}-Q_{l}|^{2}}\Big{)}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_{0,k}(y)=V_{l,k}(y)+M_{k}\tilde{w}_{l,k}(y)$

and

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

Thus for $s\neq l$ we have

(3.57)
$\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

$\displaystyle V_{s,k}(y)-V_{l,k}(y)$ $\displaystyle=$
$\displaystyle\mu_{l}^{k}-\mu_{s}^{k}+2\log\frac{h_{s}}{h_{l}}$
$\displaystyle+2\log(1+\frac{\frac{D}{h_{l}}e^{-\mu_{l}^{k}}-\frac{D}{h_{s}}e^{-\mu_{s}^{k}}+|y^{N+1}-1-p_{l}|^{2}-|y^{n+1}-1-p_{s}|^{2}}{\frac{D}{h_{s}}e^{-\mu_{s}^{k}}|y^{N+1}-e_{1}-p_{s}|}$

To evaluate the above, we use

$\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}.$

and set

$\sigma_{k}:=(\mu_{l}^{k}-\mu_{s}^{k}+2\log\frac{h_{s}}{h_{l}})/M_{k}.$

Then we see that $\sigma_{k}=o(1)$ and
$\frac{D}{h_{l}}e^{-\mu_{l}^{k}}-\frac{D}{h_{s}}e^{-\mu_{s}^{k}}=o(M_{k})\epsilon_{k}^{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}+\sigma_{k}+o(\epsilon_{k}).$

We are going to derive a contradiction based on (3.57). For this purpose we
choose $s=0$ in (3.57), which means $\tilde{w}_{s,k}=o(\epsilon_{k})$, then by
taking $|p_{l}^{k}|\sim\epsilon_{k}M_{k}$ we have

(3.58) $\displaystyle
4Re(\frac{y^{N+1}-1}{|y^{N+1}-1|^{2}}\frac{\bar{p}_{l}}{\epsilon_{k}M_{k}})\epsilon_{k}+\sigma_{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}|}+O(\delta_{k}^{2}))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 $|\sigma_{k}|\geq C\epsilon_{k}$, it is easy to see that (3.58) does not
hold for different $|y|\sim 1$. When $\sigma_{k}=o(1)\sigma_{k}$ we observe
$D_{k}$. If $|D_{k}|\geq C\epsilon_{k}$ one can also see the contradiction for
different $y$. Finally for $\sigma_{k}=o(\epsilon_{k}),D_{k}=o(\epsilon_{k})$,
the equality still cannot hold if choose $|y|>>1$ because the last two terms
of (3.58) will majorize the left-hand-side. Lemma 3.3 is established. $\Box$

Proposition 3.1 is an immediate consequence of Lemma 3.3. $\Box$.

##  4\. Proof of the main theorems

First we make a simple observation. If we let

$\hat{w}_{k}=w_{0,k}/\delta_{k}^{*}.$

It is easy to see that $\hat{w}_{k}$ tends to a harmonic function away from
the $N+1$ local maximums of $v_{0,k}$. Since $\hat{w}_{k}$ is bounded, the
global harmonic function $\hat{w}_{k}$ tends to has to be a constant. If we
focus on the neighborhood of $Q_{0}^{k}$, by standard Fourier analysis as
before we see that this constant is zero. Thus over compact subset away from
the $N+1$ local maximum points of $v_{0,k}$ we have $\hat{w}_{k}=o(1)$ and
$\nabla\hat{w}_{k}=o(1)$, which means in this region
$w_{0,k}=o(\delta_{k}^{*})$ and $\nabla w_{0,k}=o(\delta_{k}^{*})$.

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.1)
$\displaystyle\int_{\Omega_{s,k}}\partial_{\xi}(|y|^{2N}h_{0,k}(y))e^{v_{0,k}}-\int_{\partial\Omega_{s,k}}e^{v_{0,k}}|y|^{2N}h_{0,k}(y)(\xi\cdot\nu)$
$\displaystyle=\int_{\partial\Omega_{s,k}}(\partial_{\nu}v_{0,k}\partial_{\xi}v_{0,k}-\frac{1}{2}|\nabla
v_{0,k}|^{2}(\xi\cdot\nu))dS.$

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

(4.2)
$\displaystyle\int_{\Omega_{s,k}}\partial_{\xi}(|y|^{2N}h_{0,k}(Q_{0}^{k}))e^{V_{0,k}}-\int_{\partial\Omega_{s,k}}e^{V_{k}}|y|^{2N}h_{0,k}(Q_{0}^{k})(\xi\cdot\nu)$
$\displaystyle=\int_{\partial\Omega_{s,k}}(\partial_{\nu}V_{k}\partial_{\xi}V_{0,k}-\frac{1}{2}|\nabla
V_{0,k}|^{2}(\xi\cdot\nu))dS.$

Using $w_{0,k}=v_{0,k}-V_{0,k}$ and $|w_{0,k}(y)|\leq C\delta_{k}^{*}$ we have

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

If we just use crude estimate: $\nabla w_{0,k}=o(\delta_{k}^{*})$, we have

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

The difference on the second terms is minor: If we use the expansion of
$v_{0,k}=V_{0,k}+w_{0,k}$ around $e_{1}$, it is easy to obtain

$\int_{\partial\Omega_{s,k}}e^{v_{0,k}}|y|^{2N}h_{0,k}(y)(\xi\cdot\nu)-\int_{\partial\Omega_{s,k}}h_{0,k}(Q_{0}^{k})e^{V_{0,k}}|y|^{2N}(\xi\cdot\nu)=o(\delta_{k}^{*}).$

To evaluate the first term, we use

(4.3)
$\displaystyle\int_{\Omega_{s,k}}\partial_{\xi}(|y|^{2N}h_{0,k}(y))e^{v_{0,k}}$
$\displaystyle=$
$\displaystyle\int_{\Omega_{s,k}}\partial_{\xi}(|y|^{2N}(h_{0,k}(Q_{s}^{k})+\nabla
h_{0,k}(Q_{s}^{k})(y-Q_{s}^{k}))e^{V_{0,k}}(1+w_{0,k})dy+o(\delta_{k}^{*})$
$\displaystyle=$
$\displaystyle\int_{\Omega_{s,k}}(\partial_{\xi}(|y|^{2N})h_{0,k}(Q_{s}^{k})e^{V_{0,k}}+\partial_{\xi}h_{0,k}(Q_{s}^{k})|y|^{2N}e^{V_{0,k}}$
$\displaystyle+\partial_{\xi}(|y|^{2N})h_{0,k}(Q_{s}^{k})e^{V_{0,k}}w_{0,k})dy+o(\delta_{k}^{*}).$

where in the evaluation above the second derivative of $h_{0,k}$ contributes
$o(\delta_{k}^{*})$ because its magnitude is $O(\delta_{k}^{*})$ but it has a
scaling factor $O(|y-Q_{s}^{k}|^{2}$ which makes the whole term
$o(\delta_{k}^{*})$.

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

$\Delta
w_{0,k}+e^{V_{0,k}}|y|^{2N}h_{0,k}(Q_{s}^{k})w_{0,k}=(h_{0,k}(Q_{s}^{k})-h_{0,k}(y))|y|^{2N}e^{V_{0,k}}+D_{2}h$

where $D_{2}h$ represents the term with the second derivatives of $h_{0,k}$,
since it only leads to $o(\delta_{k}^{*})$ in evaluation, we omit its specific
form.

From integration by parts we have

(4.4)
$\displaystyle\int_{\Omega_{s,k}}\partial_{\xi}(|y|^{2N})h_{0,k}(Q_{s}^{k})e^{V_{0,k}}w_{0,k}$
$\displaystyle=$ $\displaystyle
2N\int_{\Omega_{s,k}}|y|^{2N-2}y_{\xi}h_{0,k}(Q_{s}^{k})e^{V_{0,k}}w_{0,k}$
$\displaystyle=$ $\displaystyle
2N\int_{\Omega_{s,k}}\frac{y_{\xi}}{|y|^{2}}(-\Delta
w_{0,k}+(h_{0,k}(y)-h_{0,k}(Q_{s}^{k}))|y|^{2N}e^{V_{k}}dy+o(\delta_{k}^{*})$
$\displaystyle=$ $\displaystyle
2N\int_{\partial\Omega_{s,k}}(\partial_{\nu}(\frac{y_{\xi}}{|y|^{2}})w_{0,k}-\partial_{\nu}w_{0,k}\frac{y_{\xi}}{|y|^{2}})+o(\delta_{k}^{*})=o(\delta_{k}^{*})$

where we have used $\nabla w_{0,k},w_{0,k}=o(\delta_{k}^{*})$ on
$\partial\Omega_{s,k}$. For the second term on the right hand side of (4.3),
By (3.12) we have, for some $0\leq s\leq N$ and $\xi\in\mathbb{S}^{1}$,

(4.5)
$\int_{\Omega_{s,k}}\partial_{\xi}h_{0,k}(Q_{s})|y|^{2N}e^{V_{s,k}}=c\delta_{k}^{*}e^{i\beta_{s}}\cdot\xi+o(\delta_{k}^{*}).$

for some $c\neq 0$. Using (4.4) and (4.5) in the difference between (4.1) and
(4.2), we have

$\delta_{k}^{*}=o(\delta_{k}^{*}).$

This contradiction finishes the proof of Theorem 2.1. $\Box$

Proof of Theorem 1.2: First it is a well known fact that all the blowup points
are in the interior of $\Omega$ (see [33]). Let $p$ be a blowup point of
$\\{\mathfrak{u}_{k}\\}$, if $p$ is a regular point or a non-quantized
singular point, the blowup solutions are known to be simple [3, 44]. So we
only consider the case that $q$ is a quantized singular source. Obviously we
can use a fundamental solution to reduce the study to equation (2.2). Because
there is at least one other blowup point, we can see that around $p$, the
harmonic function that eliminates the oscillation of $u_{k}$ around $p$ does
not tend to $0$, because around the other blowup point, the harmonic function
is tending to infinity. Theorem 1.2 is an immediate consequence of Theorem
2.1. $\Box$

Proof of Theorem 1.1: We consider the blowup solutions of
$u_{k}+\log\lambda_{k}$. It is well known that if $\lambda_{k}$ does not tend
zero, $u_{k}$ is uniformly bounded. So we use $\hat{u}_{k}$ to denote
$u_{k}+\log\lambda_{k}$ and we assume that $\hat{u}_{k}$ is a sequence of
blowup solutions. If there are at least two blowup points, the argument for
Theorem 1.2 can be applied here to say that each blowup point is simple. So we
only consider the case that there is one blowup point. Without loss of
generality we take this blowup point to be the origin. If the origin is a
regular point or a singular point with non-quantized Dirac mass, the blowup
sequence is well known to be simple. Thus we assume that at $0$ there is a
quantized singular source $4\pi N\delta_{0}$. If the domain is a disk centered
at the origin, the proof of Theorem 2.1 also applies in this case that the
only blowup point is simple. So the only case is that $\Omega$ is not a disk
centered around the origin. In this case we assume that
$B_{1}\subset\Omega\subset B_{R}$ for some $R>1$. Now we recall the asymptotic
expansion of

$V_{k}(y)=\mu_{k}-2\log(1+\frac{e^{\mu_{k}}}{D}|y^{N++1}-1-p_{k}|^{2})$

where $|p_{k}|=O(\mu_{k}e^{-\mu_{k}})$. Thenf or $r$ large, say, comparable to
$\delta_{k}^{-1}$,

$\displaystyle V_{k}(y)=-\mu_{k}+2\log D-4(N+1)\log
r+\frac{4}{r^{N+1}}\cos((N+1)\theta)$
$\displaystyle+O(\delta_{k}^{N+1}\mu_{k}e^{-\mu_{k}})+O(\delta_{k}^{2N+2}).$

We define a harmonic function $\phi_{0,k}$ that eliminates all the oscillation
of $V_{k}$ on $\partial\Omega_{k}$. Consider the Fourier expansion of
$\phi_{0,k}$:

$\sum_{n=1}^{\infty}(a_{n}^{k}\cos(n\theta)+b_{n}^{k}\sin(n\theta))r^{n}.$

Note that there is no projection on $1$ because $\phi_{0,k}(0)=0$. Since the
value of the absolute value of $\phi_{0,k}$ on $\partial\Omega_{k}$ (which is
between $|y|=\delta_{k}^{-1}$ and $|y|=R\delta_{k}^{-1}$ ) is comparable to
$\delta_{k}^{N+1}$, we make two observations. First the restriction of
$\phi_{0,k}$ on $\partial B(0,\delta_{k}^{-1})$ is between
$-C\delta_{k}^{N+1}$ and $C\delta_{k}^{N+1}$. This is because of the
asymptotic expansion of $V_{k}$ for $|y|\sim\delta_{k}^{-1}$. This provides an
upper bound of $a_{n}^{k}$ and $b_{n}^{k}$:

$|a_{n}^{k}|+|b_{n}^{k}|\leq C\delta_{k}^{N+1+n}.$

Second, we shall prove there is a point $p_{k}\in\partial
B(0,\delta_{k}^{-1})$ such that

(4.6) $|\phi_{0,k}(p_{k})|\geq
C\delta_{k}^{N+1},\quad|p_{k}|=\delta_{k}^{-1}.$

Without loss of generality we can assume that $\partial\Omega_{k}$ is tangent
to $\partial B(0,\delta_{k})$. It is more convenient to consider this problem
before scaling: Let $\phi$ be a harmonic function on $\Omega$ with $\phi(0)=0$
and the value of $\phi$ on $\partial\Omega$ can be described as

$\phi(x)=f(|x|)\cos((N+1)\theta),\quad x=|x|e^{i\theta}\in\partial\Omega,$

where $f(r)\sim 1$ for $1<r<R$. Let $e^{i\theta_{0}}$ be the tangent point
between $\partial B_{1}$ and $\partial\Omega$, if $\cos((N+1)\theta_{0}\neq
0$, we already have $\phi(e^{i\theta_{0}})\neq 0$, which certainly gives
$|\phi(e^{i\theta_{0}})|\geq C$ for some $C>0$ independent of $k$. Once this
is done, after scaling and normalization with $\delta_{k}^{N+1}$ we have (4.6)
for $p_{k}=\delta_{k}^{-1}e^{i\theta_{0}}$ and
$\phi_{0,k}(y)=\delta_{k}^{N+1}\phi(\delta_{k}y)$.

If $\cos((N+1)\theta_{0})=0$, we just consider a neighborhood of
$e^{i\theta_{0}}$. If $\theta$ is just slightly different from $\theta_{0}$,
the magnitude of $\phi$ on $|x|e^{i\theta}\in\partial\Omega$ is comparable to
$1$ because $f(|x|)$ is comparable to $1$. The curvature of $\partial\Omega$
is bounded, it is easy to use the Green’s representation formula of $\phi$ to
find $x_{1}=e^{i\theta_{1}}\in\partial B_{1}$, close to $e^{i\theta_{0}}$ such
that $|\phi(x_{1})|\geq C$ for some $C>0$. Scaling back to
$\partial\Omega_{k}$ we have $p_{k}=\delta_{k}^{-1}e^{i\theta_{1}}$ and (4.6)
holds.

Based on these two observations we set

$\delta_{k}^{*}=\sum_{n=1}^{N+1}(|a_{n}^{k}|+|b_{n}^{k}|),$

we have

$C\delta_{k}^{2N+2}\leq\delta_{k}^{*}\leq C\delta_{k}^{N+2}.$

By setting $h_{0,k}(y)=e^{\phi_{0,k}}$ as before, we have, using the same
proof of Lemma 3.1, that

$|\nabla h_{0,k}(Q_{s}^{k})|\geq C\delta_{k}^{*},$

for at least one $Q_{s}^{k}$. Then the remaining part of the proof is similar
to that of Theorem 1.2. Theorem 1.1 is established. $\Box$.

##  5\. Appendix A

In this section we derive (3.47) and (3.48). We use simplified notations:

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

where $D_{s}=\frac{8(N+1)^{2}}{h_{s}}$, $h_{s}=h_{0,k}(Q_{s}^{k})$.

$V_{s}-V_{l}=\mu_{s}-\mu_{l}+2\log(1+A)=\mu_{s}-\mu_{l}+2A-A^{2}+O(A^{3}),$

where

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

(5.1) $\displaystyle\frac{e^{\mu_{l}}}{D_{l}}$
$\displaystyle=\frac{e^{\mu_{s}}}{D_{s}}e^{\mu_{l}-\mu_{s}}\frac{D_{s}}{D_{l}}$
$\displaystyle=\frac{e^{\mu_{s}}}{D_{s}}(1+(\mu_{l}-\mu_{s})+O(\mu_{l}-\mu_{s})^{2})(1+\frac{h_{l}-h_{s}}{h_{s}}+O(\delta_{k}^{2*}))$
$\displaystyle=\frac{e^{\mu_{s}}}{D_{s}}(1+\mu_{l}-\mu_{s}+\frac{h_{l}-h_{s}}{h_{s}}+O((\mu_{l}-\mu_{s})^{2}+\delta_{k}^{2*})).$

where we used $(D_{s}-D_{l})/D_{l}=(h_{l}-h_{s})/h_{s}=O(\delta_{k}^{*})$ and
we use $\delta_{k}^{2*}$ to denote $(\delta_{k}^{*})^{2}$.

Using this we now write $A$ as

$\displaystyle A$
$\displaystyle=\bigg{(}\frac{e^{\mu_{s}}}{D_{s}}(|y^{N+1}-1-p_{l}|^{2}-|y^{N+1}-1-p_{s}|^{2})$
$\displaystyle+\frac{e^{\mu_{s}}}{D_{s}}(\mu_{l}-\mu_{s}+\frac{h_{l}-h_{s}}{h_{s}}+E)|y^{N+1}-1-p_{l}|^{2}\bigg{)}/\bigg{(}1+\frac{e^{\mu_{s}}}{D_{s}}|y^{N+1}-1-p_{s}|^{2}\bigg{)}$

where $E=O((\mu_{l}-\mu_{s})^{2}+\delta_{k}^{2*})$.

The derivation of (3.47) is as follows: Recall that the expression of $Q_{s}$
is in (3.44):

$Q_{s}=e^{i\beta_{s}}(1+\frac{p_{s}}{N+1})+O(|p_{s}|^{2}).$

and $y=Q_{s}+\epsilon_{k}z$. Direct computation shows

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

Then

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

With these expressions one can check easily that the expression of $B$ in
(3.47) holds:

$B=1+\frac{h_{s}}{8}|z+\frac{N}{2}\epsilon_{k}z^{2}e^{-i\beta_{s}}|^{2}+O(\epsilon_{k}^{2}\mu_{k})|z|^{2}.$

Then the expression of $A$ is

$\displaystyle
A=\frac{h_{s}}{8}\bigg{(}2Re(z+\frac{N}{2}\epsilon_{k}z^{2}e^{-i\beta_{s}})\frac{\bar{p}_{s}-\bar{p}_{l}}{(N+1)\epsilon_{k}}e^{-i\beta_{s}}+\frac{|p_{s}-p_{l}|^{2}}{(N+1)^{2}\epsilon_{k}^{2}}$
$\displaystyle+(\mu_{l}-\mu_{s}+\frac{h_{l}-h_{s}}{h_{s}}+E_{\mu})|z+\frac{N}{2}\epsilon_{k}z^{2}e^{-i\beta_{s}}+O(\epsilon_{k}^{2}\mu_{k}z)|^{2}\bigg{)}/B.$

where $E_{\mu}=O(|\mu_{l}-\mu_{s}|^{2}+\delta_{k}^{2*})$. Note that
$D_{s}=o(M_{k})$, so $E_{\mu}=o(M_{k}^{2})$.

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

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

$\frac{1}{32(N+1)^{2}}\bigg{(}|z|^{2}(\frac{|p_{s}-p_{l}|}{\epsilon_{k}})^{2}\big{(}1+2\cos(2\theta-2\theta_{sl}-2\beta_{s})\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 (3.47) by direct computation.

##  6\. Appendix B

In this section we prove (3.56) based on (3.48). The terms of $\phi_{1}$ and
$\phi_{2}$ lead to $o(\epsilon_{k})$ The integrations involving $\phi_{3}$ and
$\phi_{4}$ provide 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

(6.1)
$\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

(6.2)
$\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.$

Based on (6.1) and (6.2), the expressions of $\phi_{1}$, $\phi_{2}$ and $B$ in
(3.48) lead to

$\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_{2}}{M_{k}}B^{-2}=o(\epsilon_{k}).$

The integrations involving $\phi_{3}$ and $\phi_{4}$ lead to the expression of
$D_{s,l}^{k}$ in (3.56). Thus (3.56) holds.

##  7\. Appendix C: Location of local maximums of $v_{k}$

In this section we use standard arguments to prove that the local maximum
points of $v_{k}$ is evenly distributed around the origin. Recall that $v_{k}$
satisfies (3.3). Here we abuse the notation by assuming that $Q_{1}^{k}=e_{1}$
is one local maximum of $v_{k}$ and we use $\mu_{k}$ instead of
$\bar{\mu}_{k}$ to denote $v_{k}(e_{1})$ for convenience.

Now we consider $v_{k}$ around $Q_{l}^{k}$. Using the results in [15, 45, 21]
we have, for $v_{k}$ in $B(Q_{l}^{k},\epsilon)$, the following gradient
estimate:

(7.1)
$2N\frac{\tilde{Q}_{l}^{k}}{|\tilde{Q}_{l}^{k}|^{2}}+\nabla\phi_{l}^{k}(\tilde{Q}_{l}^{k})=O(\mu_{k}e^{-\mu_{k}}),$

where $\phi_{l}^{k}$ is the harmonic function that eliminates the oscillation
of $v_{k}$ on $\partial B(Q_{l}^{k},\epsilon)$ and $\tilde{Q}_{l}^{k}$ is the
maximum of $v_{k}-\phi_{l}^{k}$ that satisfies

(7.2) $\tilde{Q}_{l}^{k}-Q_{l}^{k}=O(e^{-\mu_{k}}).$

Using (7.2) in (7.1) we have

(7.3)
$2N\frac{Q_{l}^{k}}{|Q_{l}^{k}|^{2}}+\nabla\phi_{l}^{k}(Q_{l}^{k})=O(\mu_{k}e^{-\mu_{k}}).$

The first estimate of $\nabla\phi_{l}^{k}(Q_{l}^{k})$ is

######  Lemma 7.1.

For $l=0,...,N$,

(7.4) $\nabla\phi_{l}^{k}(Q_{l}^{k})=-4\sum_{m=0,m\neq
l}^{N}\frac{Q_{l}^{k}-Q_{m}^{k}}{|Q_{l}^{k}-Q_{m}^{k}|^{2}}+E$

where

$E=O(\mu_{k}e^{-\mu_{k}}).$

Proof of Lemma 7.1:

From the expression of $v_{k}$ on $\Omega_{k}=B(0,\tau\delta_{k}^{-1})$ we
have, for $y$ away from bubbling disks,

(7.5) $\displaystyle v_{k}(y)$
$\displaystyle=v_{k}|_{\partial\Omega_{k}}+\int_{\Omega_{k}}G(y,\eta)|\eta|^{2N}e^{v_{k}(\eta)}d\eta$
$\displaystyle=v_{k}|_{\partial\Omega_{k}}+\sum_{l=0}^{N}G(y,Q_{l}^{k})\int_{B(Q_{l},\epsilon)}|\eta|^{2N}e^{v_{k}}d\eta$
$\displaystyle+\sum_{l}\int_{B(Q_{l},\epsilon)}(G(y,\eta)-G(y,Q_{l}^{k}))|\eta|^{2N}e^{v_{k}}d\eta+O(\mu_{k}e^{-\mu_{k}}).$

We have, by standard estimates

$v_{k}(y)=v_{k}|_{\partial\Omega_{k}}-4\sum_{l=0}^{N}\log|y-Q_{l}^{k}|+8\pi\sum_{l=0}^{N}H(y,Q_{l}^{k})+O(\mu_{k}e^{-\mu_{k}}).$

The harmonic function that kills the oscillation of $v_{k}$ around $Q_{m}^{k}$
is

$\displaystyle\phi_{m}^{k}=-4\sum_{l=0,l\neq
m}^{N}(\log|y-Q_{l}^{k}|-\log|Q_{m}^{k}-Q_{l}^{k}|)$
$\displaystyle+8\pi\sum_{l=0}^{N}(H(y,Q_{l}^{k})-H(Q_{m}^{k},Q_{l}^{k}))+O(\mu_{k}e^{-\mu_{k}}).$

The corresponding estimate for $\nabla\phi_{m}^{k}$ is

$\nabla\phi_{m}^{k}(Q_{m}^{k})=-4\sum_{l=0,l\neq
m}^{N}\frac{Q_{m}^{k}-Q_{l}^{k}}{|Q_{m}^{k}-Q_{l}^{k}|^{2}}+8\pi\sum_{l=0}^{N}\nabla_{1}H(Q_{m}^{k},Q_{l}^{k})+O(\mu_{k}e^{-\mu_{k}}).$

where $\nabla_{1}$ stands for the differentiation with respect to the first
component. From the expression of $H$, we have

(7.6) $\displaystyle\nabla_{1}H(Q_{m}^{k},Q_{l}^{k})$
$\displaystyle=\frac{1}{2\pi}\frac{Q_{m}^{k}-\tau^{2}\delta_{k}^{-2}Q_{l}^{k}/|Q_{l}^{k}|^{2}}{|Q_{m}^{k}-\tau^{2}\delta_{k}^{-2}Q_{l}^{k}/|Q_{l}^{k}|^{2}|^{2}}$
$\displaystyle=\frac{1}{2\pi}\tau^{-2}\delta_{k}^{2}\frac{\tau^{-2}\delta_{k}^{2}Q_{m}^{k}-Q_{l}^{k}/|Q_{l}^{k}|^{2}}{|Q_{l}^{k}/|Q_{l}^{k}|^{2}-\tau^{-2}\delta_{k}^{2}Q_{m}^{k}|^{2}}$
$\displaystyle=-\frac{1}{2\pi}\tau^{-2}\delta_{k}^{2}e^{\frac{2\pi
il}{N+1}}+O(\sigma_{k}\delta_{k}^{2}).$

where $\sigma_{k}=\max_{l}|Q_{l}^{k}-e^{\frac{2\pi il}{N+1}}|$. Later we shall
obtain more specific estimate of $\sigma_{k}$.

Thus

(7.7) $\displaystyle\nabla\phi_{m}^{k}(Q_{m}^{k})$ $\displaystyle=$
$\displaystyle-4\sum_{l=0,l\neq
m}^{N}\frac{Q_{m}^{k}-Q_{l}^{k}}{|Q_{m}^{k}-Q_{l}^{k}|^{2}}-4\tau^{-2}\delta_{k}^{2}\sum_{l=0}^{N}e^{\frac{2\pi
il}{N+1}}+O(\sigma_{k}\delta_{k}^{2})+O(\mu_{k}e^{-\mu_{k}})$ $\displaystyle=$
$\displaystyle-4\sum_{l=0,l\neq
m}^{N}\frac{Q_{m}^{k}-Q_{l}^{k}}{|Q_{m}^{k}-Q_{l}^{k}|^{2}}+O(\sigma_{k}\delta_{k}^{2})+O(\mu_{k}e^{-\mu_{k}})$

where we have used $\sum_{l=0}^{N}e^{2\pi li/(N+1)}=0$. Since we don’t have
the estimate of $\sigma_{k}$ now we have

$\nabla\phi_{m}^{k}(Q_{m}^{k})=-4\sum_{l=0,l\neq
m}^{N}\frac{Q_{m}^{k}-Q_{l}^{k}}{|Q_{m}^{k}-Q_{l}^{k}|^{2}}+O(\mu_{k}e^{-\mu_{k}})+O(\sigma_{k}\delta_{k}^{2}).$

Lemma 7.1 is established. $\Box$

Using

$E_{1}=O(\mu_{k}e^{-\mu_{k}})+O(\sigma_{k}\delta_{k}^{2}).$

The Pohozaev identity around $Q_{l}^{k}$ now reads

$-4\sum_{j=0,j\neq
l}^{N}\frac{Q_{l}^{k}-Q_{j}^{k}}{|Q_{l}^{k}-Q_{j}^{k}|^{2}}+2N\frac{Q_{l}^{k}}{|Q_{l}^{k}|^{2}}=E_{1}.$

We have, treating every term as a complex number,

$N\frac{1}{\bar{Q}_{l}^{k}}=2\sum_{j=0,j\neq
l}^{N}\frac{1}{\bar{Q}_{l}^{k}-\bar{Q}_{j}^{k}}+E_{1},$

where $\bar{Q}_{l}^{k}$ is the conjugate of $Q_{l}^{k}$. Thus

(7.8) $N=2\sum_{j=0,j\neq l}^{N}\frac{Q_{l}^{k}}{Q_{l}^{k}-Q_{j}^{k}}+E_{1}.$

Let $\beta_{l}=2\pi l/(N+1)$, we write $Q_{l}^{k}=e^{i\beta_{l}}+p_{l}^{k}$
for $p_{l}^{k}\to 0$. Then we write the first term on the right hand side of
(7.8) as

$\displaystyle\frac{Q_{l}^{k}}{Q_{l}^{k}-Q_{j}^{k}}=\frac{e^{i\beta_{l}}+p_{l}^{k}}{e^{i\beta_{l}}-e^{i\beta_{j}}+p_{l}^{k}-p_{j}^{k}}$
$\displaystyle=$
$\displaystyle\frac{e^{i\beta_{l}}+p_{l}^{k}}{(e^{i\beta_{l}}-e^{i\beta_{j}})(1+(p_{l}^{k}-p_{j}^{k})/(e^{i\beta_{l}}-e^{i\beta_{j}}))}$
$\displaystyle=$
$\displaystyle\frac{e^{i\beta_{l}}}{e^{i\beta_{l}}-e^{i\beta_{j}}}+\frac{p_{l}^{k}}{e^{i\beta_{l}}-e^{i\beta_{j}}}-\frac{e^{i\beta_{l}}}{(e^{i\beta_{l}}-e^{i\beta_{j}})^{2}}(p_{l}^{k}-p_{j}^{k})+O(\sigma_{k}^{2})$
$\displaystyle=$
$\displaystyle\frac{e^{i\beta_{l}}}{e^{i\beta_{l}}-e^{i\beta_{j}}}+\frac{e^{i\beta_{l}}p_{j}^{k}-e^{i\beta_{j}}p_{l}^{k}}{(e^{i\beta_{l}}-e^{i\beta_{j}})^{2}}+O(\sigma_{k}^{2}).$

Using

(7.9) $N=2\sum_{j=0,j\neq
l}^{N}\frac{e^{i\beta_{l}}}{e^{i\beta_{l}}-e^{i\beta_{j}}},$

we write (7.8) as

(7.10) $\sum_{j=0,j\neq
l}^{N}\frac{e^{i\beta_{l}}p_{j}^{k}-e^{i\beta_{j}}p_{l}^{k}}{(e^{i\beta_{l}}-e^{i\beta_{j}})^{2}}=E_{1}+O(\sigma_{k}^{2})$

for $l=0,1,2,....,N$. For convenience we set

$p_{l}^{k}=e^{i\beta_{l}}m_{l}^{k}\quad\mbox{ and
}\quad\beta_{jl}=\beta_{j}-\beta_{l}$

to reduce (7.10) to

(7.11) $\displaystyle\sum_{j=0,j\neq
l}^{N}\frac{e^{i\beta_{jl}}m_{j}^{k}}{(1-e^{i\beta_{jl}})^{2}}-\bigg{(}\sum_{j=0,j\neq
l}^{N}\frac{e^{i\beta_{jl}}}{(1-e^{i\beta_{jl}})^{2}}\bigg{)}m_{l}^{k}$
$\displaystyle=E+O(\sigma_{k}^{2})+O(\delta_{k}\sigma_{k})$

for $l=0,1.....,N$. It is easy to verify that

(7.12)
$\frac{e^{i\theta}}{(1-e^{i\theta})^{2}}=\frac{1}{2(\cos\theta-1)}=(-\frac{1}{4})\frac{1}{\sin^{2}(\theta/2)}.$

To deal with coefficients of $m_{j}^{k}$ in (7.11) we set

$d_{j}=\frac{1}{\sin^{2}(\frac{j\pi}{N+1})},\quad j=1,...,N$

and

$D=\sum_{j=0,j\neq l}^{N}d_{|j-l|}.$

Since $d_{l}=d_{N+1-l}$ it is easy to check that $D$ does not depend on $l$:

(7.13)
$D=\sum_{k=1}^{N}d_{k}=\sum_{k=1}^{N}\frac{1}{\sin^{2}(\frac{k\pi}{N+1})}=\frac{N^{2}+2N}{3}.$

Now (7.11) can be written as

(7.14) $-\sum_{j\neq
l,j=0}^{N}d_{|j-l|}m_{j}^{k}+Dm_{l}^{k}=E_{1}+O(\sigma_{k}^{2}),\quad
l=0,....,N.$

For $l=0$, we have $\beta_{0}=0$ and $m_{0}^{k}=0$. Thus from (7.14) we have

(7.15) $-\sum_{j=1}^{N}d_{j}m_{j}^{k}=E_{1}+O(\sigma_{k}^{2}).$

If we take $(m^{k}_{1},...,m^{k}_{n})$ as unknowns in (7.14), the last $N$
equations of (7.14) ( for $l=1,...,N$) can be written as

(7.16) $A\left(\begin{array}[]{c}m^{k}_{1}\\\ m^{k}_{2}\\\ \vdots\\\
m^{k}_{N}\end{array}\right)=E_{1}+O(\sigma_{k}^{2}).$

where

$A=\left(\begin{array}[]{cccc}D&-d_{1}&...&-d_{N-1}\\\
-d_{1}&D&...&-d_{N-2}\\\ \vdots&\vdots&...&\vdots\\\
-d_{N-1}&-d_{N-2}&...&D\end{array}\right)$

Since $D=|d_{1}|+...+|d_{N}|$ and each $d_{i}>0$, we see that the matrix is
invertible, thus $|m^{k}_{i}|=E$ for all $i$.

One consequence about $v_{0,k}$ is that the local maximum points of $v_{0,k}$,
being perturbed by the amount of $O(\delta_{k}^{*})$, are
$O(e^{-\mu_{k}}\delta_{k}^{*})$ away from the corresponding locations of the
local maximum points of $v_{k}$. The reason is $v_{k}$ is non-degenerate near
each blowup point.

## References

  * [1] M. Ahmedou; L. Wu; L. Zhang. Classification and a priori estimates for the singular prescribing Q-curvature equation on 4-manifold, J. Funct. Anal. 283 (2022), no. 9, 109649.
  * [2] J. Ambjorn; 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] D. Bartolucci; C.C. Chen; C.S. Lin; G. Tarantello. Profile of blow-up solutions to mean field equations with singular data, Comm. Partial Differential Equations 29 (2004), no. 7-8, 1241-1265.
  * [4] D. Bartolucci; A. Jevnikar; Y. Lee; W. Yang. Uniqueness of bubbling solutions of mean field equations, J. Math. Pures Appl. (9) 123 (2019), 78–126.
  * [5] D. Bartolucci; G. Tarantello, Asymptotic Blow-up Analysis for Singular Liouville type Equations with Applications, J. Differential Equations 262, no. 7 (2017), 3887–3931
  * [6] D. Bartolucci; G. Tarantello. The Liouville equation with singular data: a concentration-compactness principle via a local representation formula, J. Differential Equations 185 (2002), no. 1, 161–180.
  * [7] 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.
  * [8] Luca Battaglia, Uniform bounds for solutions to elliptic problems on simply connected planar domains. Proc. Amer. Math. Soc. 147 (2019), no. 10, 4289–4299.
  * [9] L.A. Caffarelli; B. Gidas; J. Spruck. Asymptotic symmetry and local behavior of semilinear elliptic equations with critical Sobolev growth, Comm. Pure Appl. Math. 42 (1989), no. 3, 271–297.
  * [10] 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.
  * [11] 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.
  * [12] S. Chanillo; M. Kiessling. Rotational symmetry of solutions of some nonlinear problems in statistical mechanics and in geometry, Comm. Math. Phys.  160 (1994), no. 2, 217–238.
  * [13] 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.
  * [14] C.C. Chen; C.S. Lin. Mean field equation of Liouville type with singular data: topological degree, Comm. Pure Appl. Math. 68 (2015), no. 6, 887–947
  * [15] 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.
  * [16] C.C. Chen; C.S. Lin. Mean field equation of Liouville type with singular data: topological degree, Comm. Pure Appl. Math. 68 (2015), no. 6, 887–947.
  * [17] T. D’Aprile; J. Wei. Bubbling solutions for the Liouville equation with a singular source: non-simple blow-up, J. Funct. Anal. 279 (2020), no. 6, 108605, 38 pp.
  * [18] 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.
  * [19] M. del Pino, P. Esposito, M. Musso, Two-dimensional Euler flows with concentrated vorticities, Trans. Amer. Math. Soc. 362 (2010), no. 12, 6381–6395.
  * [20] A. Eremenko; A. Gabrielov; V. Tarasov. Metrics with conic singularities and spherical polygons, Illinois J. Math 58 (2014), no. 3, 739–755.
  * [21] M. Gluck. Asymptotic behavior of blow up solutions to a class of prescribing Gauss curvature equations. Nonlinear Anal. 75 (2012), no. 15, 5787–5796.
  * [22] Y. Gu; L. Zhang. Degree counting theorems for singular Liouville systems, Ann. Sc. Norm. Super. Pisa Cl. Sci. (5) Vol. XXI (2020), 1103–1135.
  * [23] Y. Gu; L. Zhang. Structure of bubbling solutions of Liouville systems with negative singular sources, preprint 2022, https://arxiv.org/abs/2112.10031
  * [24] R. Jackiw; E.J. Weinberg. Self-Dual Chern Simons vortices, Phys. Rev. Lett. 64 (1990), 2234–2237.
  * [25] T.J. Kuo; C.S. Lin. Estimates of the mean field equations with integer singular sources: non-simple blowup, J. Differential Geom. 103 (2016), no. 3, 377–424.
  * [26] Y.Y. Li. Harnack type inequality: the method of moving planes, Comm. Math. Phys. 200 (1999), no. 2, 421–444.
  * [27] Y.Y. Li, I. Shafrir, Blow-up analysis for solutions of $-\Delta u=V(x){e^{u}}$ in dimension two, Indiana Univ. Math. J. 43 (4) (1994), 1255–1270.
  * [28] C.S. Lin; J. Wei; D. Ye. Classifcation and nondegeneracy of $SU(n+1)$ Toda system, Invent. Math.  190 (2012), no.1, 169–207.
  * [29] C.S. Lin; J. Wei; L. Zhang. Local profile of fully bubbling solutions to SU(n+1) Toda systems. J. Eur. Math. Soc. (JEMS)  18 (2016), no. 8, 1707–1728.
  * [30] C.S. Lin; J. Wei; L. Zhang. Convergence rate, location and $\partial_{z}^{2}$ condition for fully bubbling solutions to SU(n+1) Toda systems, Adv. Math. 285 (2015), 1189–1224.
  * [31] C. S. Lin; S. Yan. On the mean field type bubbling solutions for Chern–Simons–Higgs equation, Adv. Math. 338 (2018), 1141–1188.
  * [32] A. Malchiodi and D. Ruiz, New improved Moser–Trudinger inequalities and singular Liouville equations on compact surfaces, Geom. Funct. Anal. 21 (2011), no. 5, 1196–1217.
  * [33] Ma, Li; Wei, Juncheng C. Convergence for a Liouville equation. Comment. Math. Helv. 76 (2001), no. 3, 506–514.
  * [34] M. Nolasco; G. Tarantello. Double vortex condensates in the Chern-Simons-Higgs theory, Calc. Var. Partial Differential Equations 9 (1999), no. 1, 31–94.
  * [35] J. Prajapat; G. Tarantello. On a class of elliptic problems in $\mathbb{R}^{2}$: symmetry and uniqueness results. Proc. Roy. Soc. Edinburgh Sect. A 131 (2001), no. 4, 967–985.
  * [36] M. Troyanov. Prescribing curvature on compact surfaces with conical singularities, Trans. Amer. Math. Soc. 324 (1991) 793–821.
  * [37] M. Troyanov. Metrics of constant curvature on a sphere with two conical singularities, Lect. Notes Math., 1410, Springer, NY, 1989, 296-308.
  * [38] J. Wei; L. Wu; L. Zhang. Estimates of bubbling solutions of $SU(3)$ Toda systems at critical parameters - Part 2, preprint, https://arxiv.org/abs/2206.08114
  * [39] J. Wei; L. Zhang. Estimates for Liouville equation with quantized singularities. Adv. Math.  380 (2021), 107606.
  * [40] J. Wei; L. Zhang. Vanishing estimates for Liouville equation with quantized singularities, Proc. London Math. Soc. 3 (2022), 1–26.
  * [41] J. Wei; L. Zhang. Laplacian vanishing theorem for quantized singular Liouville equation, preprint https://arxiv.org/abs/2202.10825
  * [42] L. Wu; L. Zhang. Uniqueness of bubbling solutions of mean field equations with non-quantized singularities, Commun. Contemp. Math. 23 (2021), no. 4, 2050038, 39 pp.
  * [43] Y. Yang. Solitons in Field Theory and Nonlinear Analysis, Springer-Verlag, 2001.
  * [44] 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.
  * [45] L. Zhang. Blowup solutions of some nonlinear elliptic equations involving exponential nonlinearities, Comm. Math. Phys. 268 (2006), no. 1, 105–133.
  * [46] L. Zhang. Vanishing estimates for fully bubbling solutions of $SU(n+1)$ Toda Systems at a singular source, Int. Math. Res. Not.  18 (2020), 5774–5795.