111Lei Zhang is partially supported by a Simons Foundation Collaboration Grant

# On Liouville systems at critical parameters, Part 2: Multiple bubbles

Hsin-Yuan Huang  Department of Applied Mathematics
National Chiao-Tung University
Hsinchu, Taiwan<EMAIL_ADDRESS>and  Lei Zhang  Department of
Mathematics
University of Florida
358 Little Hall P.O.Box 118105
Gainesville FL 32611-8105<EMAIL_ADDRESS>
###### Abstract.

In this paper, we continue to consider the generalized Liouville system:

$\Delta_{g}u_{i}+\sum_{j=1}^{n}a_{ij}\rho_{j}\left(\frac{h_{j}e^{u_{j}}}{\int
h_{j}e^{u_{j}}}-{1}\right)=0\quad\text{in \,}M,\quad i\in I=\\{1,\cdots,n\\},$

where $(M,g)$ is a Riemann surface $M$ with volume $1$, $h_{1},..,h_{n}$ are
positive smooth functions and $\rho_{j}\in\mathbb{R}^{+}$($j\in I$). In
previous works Lin-Zhang identified a family of hyper-surfaces $\Gamma_{N}$
and proved a priori estimates for $\rho=(\rho_{1},..,\rho_{n})$ in areas
separated by $\Gamma_{N}$. Later Lin-Zhang also calculated the leading term of
$\rho^{k}-\rho$ where $\rho\in\Gamma_{1}$ is the limit of $\rho^{k}$ on
$\Gamma_{1}$ and $\rho^{k}$ is the parameter of a bubbling sequence. This
leading term is particularly important for applications but it is very hard to
be identified if $\rho^{k}$ tends to a higher order hypersurface $\Gamma_{N}$
($N>1$). Over the years numerous attempts have failed but in this article we
overcome all the stumbling blocks and completely solve the problem under the
most general context: We not only capture the leading terms of
$\rho^{k}-\rho\in\Gamma_{N}$, but also reveal new robustness relations of
coefficient functions at different blowup points.

###### Key words and phrases:

Liouville system, asymptotic analysis, a priori estimate, classification of
solutions, method of unique continuation, Pohozaev identity, blowup phenomenon

###### 1991 Mathematics Subject Classification:

35J60, 35J55

##  1\. Introduction

Let $(M,g)$ be a compact Riemann surface with volume $1$, $h_{1},...,h_{n}$ be
positive $C^{3}$ functions on $M$, $\rho_{1},..,\rho_{n}$ be nonnegative
constants. In this article we continue our study of the following Liouville
system defined on $(M,g)$:

(1.1)
$\Delta_{g}u_{i}+\sum_{j=1}^{n}\rho_{j}a_{ij}\left(\frac{h_{j}e^{u_{j}}}{\int_{M}h_{j}e^{u_{j}}dV_{g}}-1\right)=0,\quad
i\in I:=\\{1,..,n\\}$

where $dV_{g}$ is the volume form, $\Delta_{g}$ is the Laplace-Beltrami
operator $\Delta_{g}\leq 0$. When $n=1$, equation (1.1) is the mean field
equation of the Liouville type:

(1.2)
$\Delta_{g}u+\rho\left(\frac{he^{u}}{\int_{M}he^{u}dV_{g}}-{1}\right)=0\quad\text{in
\,}M$

when $a_{11}=1$. Therefore, the Liouville system (1.1) is a natural extension
of the classical Liouville equation, which has been extensively studied for
decades because of its profound connections with various fields in geometry
and physics. Since the general form of Liouville systems includes many models
from Biology, Physics and other disciplines of sciences, it is very desirable
to study generical Liouville systems and derive common features. Recently, the
Liouville system has drawn a lot of attention because it also arises from the
stationery solutions of multi-species Patlak-Keller-Segel system[41] and self-
dual condensate solutions of Ablian Chern-Simons model with $N$ Higgs
particles[38, 31] when certain parameter tends to zero. In particular, these
two examples exhibit the bubbling phenomenon. The study of bubbling solutions
represents an essential difficulty of Liouville system and it not only impacts
the immediately related fields but also depends on the development of them.
The readers may look into the following references for closely related
discussions [1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 14, 17, 18, 19, 20, 23, 25, 26,
27, 28, 30, 35, 36, 37, 39, 40, 42, 43, 44].

For system (1.1) the Sobolev spaces for solutions are

${}^{\text{\r{}}}H^{1,n}=\,^{\text{\r{}}}H^{1}(M)\times\cdots\times\,^{\text{\r{}}}H^{1}(M)$

where

${}^{\text{\r{}}}H^{1}(M)=\left\\{u\in L^{2}(M)\,\Big{|}\,|\nabla_{g}u|\in
L^{2}(M)\text{ and }\int_{M}u\,dV_{g}=0\right\\}.$

For any $\rho=(\rho_{1},\cdots,\rho_{n})$, $\rho_{i}>0$, let $\varPhi_{\rho}$
be a nonlinear functional defined in ${}^{\text{\r{}}}H^{1,n}$ by

$\varPhi_{\rho}(u)=\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}\int_{M}\nabla_{g}u_{i}\cdot\nabla_{g}u_{j}dV_{g}-\sum_{j=1}^{n}\rho_{j}\log\int_{M}h_{j}e^{u_{j}}dV_{g}$

where $(a^{ij})_{n\times n}$ is the inverse of $A=(a_{ij})_{n\times n}$, $I$
is the set of indexes: $I=\\{1,...,n\\}$. It is easy to see that equation
(1.1) is the Euler-Lagrangian equation of $\varPhi_{\rho}$.

In [32, 33], Lin and the second author completed a degree counting program for
(1.2) under the following two assumptions on the matrix A:

$\displaystyle(H1):\quad A\mbox{ is symmetric, nonnegative, irreducible and
invertible.}$ $\displaystyle(H2):\quad a^{ii}\leq 0,\,\,\forall i\in I,\quad
a^{ij}\geq 0\,\,\forall i\neq j\in I,\quad\sum_{j=1}^{n}a^{ij}\geq
0\,\,\forall i\in I.$

Roughly speaking $(H1)$ is a rather standard assumption for Liouville systems,
$(H2)$ says the interaction between equations has to be strong. In [33] Lin-
Zhang identified a family of hypersurfaces

$\Gamma_{N}=\left\\{\rho\,\big{|}\,\rho_{i}>0,i\in
I;\,\,\Lambda_{I}(\rho)=0,\quad\Lambda_{J}(\rho)>0,\,\,\forall\,\,\emptyset\neq
J\subsetneq I,\,\,\right\\}.$

where

$\Lambda_{I}(\rho)=4\sum_{i=1}^{n}\frac{\rho_{i}}{2\pi
N}-\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}\frac{\rho_{i}}{2\pi
N}\frac{\rho_{j}}{2\pi N}$

($\Lambda_{J}(\rho)$ is understood similarly). Furthermore they proved that if
$\rho=(\rho_{1},...,\rho_{n})$ stays in a regions bounded by these critical
hypersurfaces, a priori estimate holds. Based on the a priori estimate, Lin-
Zhang proved the following degree counting formula which implies existence of
solution if the degree is not zero:

(1.3) $d_{\rho}=\left\\{\begin{array}[]{ll}1\quad\text{if
}\rho\in\mathcal{O}_{0}\\\
\frac{1}{N!}\bigg{(}(-\chi_{M}+1)...(-\chi_{M}+N)\bigg{)}\quad\text{if
}\rho\in\mathcal{O}_{N}.\end{array}\right.$

where $\chi_{M}$ is the Euler characteristic of $M$, $\mathcal{O}_{N}$ is the
region between $\Gamma_{N}$ and $\Gamma_{N+1}$. The jump of the degree from
$\mathcal{O}_{N-1}$ to $\mathcal{O}_{N}$ is contributed by the blowup
solutions to (1.1).

In this article we consider the case $\rho\to\Gamma_{N}$ from
$\mathcal{O}_{N-1}$ or $\mathcal{O}_{N}$: suppose
$\rho=(\rho_{1},..,\rho_{n})\in\Gamma_{N}$ is a limit point,
$\rho^{k}=(\rho_{1}^{k},..,\rho_{n}^{k})$ is a sequence of parameters
corresponding to bubbling solutions, the aim of this article is to identity
the leading term of $\rho^{k}-\rho$ as $k\to\infty$. Since the normal vector
at $\rho$ is proportional to

$(\sum_{j=1}^{n}a_{1j}\frac{\rho_{j}}{2\pi
N}-2,...,\sum_{j=1}^{n}a_{nj}\frac{\rho_{j}}{2\pi N}-2),$

which has all components positive (implied by $\Lambda_{J}>0$ in the
definition of $\Gamma_{N}$, see [33]), we assume

(1.4) $\frac{\rho_{i}^{k}-\rho_{i}}{\rho_{j}^{k}-\rho_{j}}\sim
1\quad\forall\quad i\neq j\in I.$

Note that $A_{k}\sim B_{k}$ means $CA_{k}\leq B_{k}\leq C_{1}A_{k}$ for some
$C,C_{1}>0$ independent of $k$. It is established in [33] that when
$\rho^{k}=(\rho_{1}^{k},..,\rho_{n}^{k})$ tends to $\rho\in\Gamma_{N}$, blowup
solutions have exactly $N$ disjoint blowup points: $p_{1},...,p_{N}$. In [34]
Lin-Zhang derived the leading term when $N=1$. It is interesting to observe in
[34] that there is one particular point $Q\in\Gamma_{1}$ such that if $\rho\to
Q$ the leading term contains local curvature at $Q$ only, but if $\rho$ tends
to any other point, the leading term is involved with global integration of
the whole manifold.

The main purpose of this article is to extend the result in [34] to
$\rho^{k}\to\Gamma_{N}$ when $N>1$. Among other things, we obtain the leading
term of $\Lambda_{I}(\rho^{k})$ as $\rho^{k}\to\Gamma_{N}$ which gives us the
sufficient conditions of a uniform bound of solutions as
$\rho^{k}\to\Gamma_{N}$. Let $p_{1}$,…,$p_{N}$ be $N$ disjoint blowup points,
which means for each $p_{t}$, there exist $p_{t}^{k}\to p_{t}$ such that
$\max_{i\in I}u_{i}^{k}(p_{t}^{k})\to\infty$, $t\in\\{1,\cdots,N\\}$. For
$h_{i}^{k}$ we assume that they are uniformly bounded by positive constants:

(1.5) $\frac{1}{C}\leq h_{i}^{k}(x)\leq C,\quad\|h_{i}^{k}\|_{C^{3}(M)}\leq C$

for all $i$ and a $C>0$ independent of $k$. Throughout the paper we use
$u^{k}=(u_{1}^{k},...,u_{n}^{k})$ to denote blowup solutions and $M_{k,t}$ to
denote the magnitude of $u^{k}$ near the blowup point $p_{t}$ and use
$\epsilon_{k,t}=e^{-\frac{1}{2}M_{k,t}}$ to measure the errors:

(1.6) $M_{k,t}=\max_{i\in I}\max_{x\in
B(p_{t},\delta_{0})}\\{u_{i}^{k}(x)-\log\int_{M}h_{i}^{k}e^{u_{i}^{k}}dV_{g},\\},\quad\epsilon_{k,t}=e^{-\frac{1}{2}M_{k,t}},$

Here we require $\delta_{0}$ to be small enough so that
$B(p_{t_{1}},\delta_{0})\cap B(p_{t_{2}},\delta_{0})=\emptyset$ for all
$t_{1}\neq t_{2}$. Let $p_{t}^{k}$ be the point where

$\max_{i\in I}\max_{x\in
B(p_{t},\delta_{0})}u_{i}^{k}(x)-\log\int_{M}h_{i}^{k}e^{u_{i}^{k}}dV_{g}$

is taken in $B(p_{t},\delta_{0})$. Under the assumptions $(H1)$ and $(H2)$, we
have a full blow-up picture in all balls (see [33]).

To understand more precise information for the blowup phenomenon to (1.1), we
shall study the convergence rate of $\Lambda_{I}(\rho^{k})$ as
$\rho^{k}\to\rho$ in terms of the magnitude of $u^{k}$. The following two
fundamental questions will be answered in this paper:

  * (1)

Are the magnitude of $u^{k}$ at different blowup points comparable to each
other?

  * (2)

What is the convergence rate of the difference of the local masses
$\sigma_{i,s}^{k}-\sigma_{i,t}^{k}$ where
$\sigma_{i,s}^{k}=\int_{B(p_{s},\delta_{0})}e^{u_{i}^{k}}dV_{g}$ and
$\sigma_{i,t}^{k}=\int_{B(p_{t},\delta_{0})}e^{u_{i}^{k}}dV_{g}$?

The first main result is to answer the first question.

######  Theorem 1.1.

Let $u^{k}\in\,\,^{\text{\r{}}}H^{1,n}(M)$ be a sequence of blowup solutions
of (1.1). Suppose $(H1),(H2)$ holds for $A$, (1.4) holds for $\rho^{k}$ and
(1.5) holds for $h_{i}^{k}$. Then

(1.7) $|M_{k,s}-M_{k,t}|=O(1),\quad\text{ for
}\,\,s\not=t,\,\,s,t\in\\{1,\cdots,N\\}.$

Here, $O(1)$ is independent of $k$.

Note that without knowing $|M_{k,s}-M_{k,t}|=O(1)$ we do not even know if
$O(\epsilon_{k,s})=O(\epsilon_{k,t})$. We thus can use
$\epsilon_{k}=e^{-\frac{1}{2}M_{k}}$ to measure the errors, where

$M_{k}=\max_{i\in I}\max_{x\in
M}\\{u_{i}^{k}(x)-\log\int_{M}h_{i}^{k}e^{u_{i}^{k}}dV_{g},\\}.$

Let

(1.8) $m_{i}=\sum_{j=1}^{n}a_{ij}\frac{\rho_{j}}{2\pi N},\quad i=1,..,n,\quad
m=\min_{i\in I}m_{i}.$

Here we note that either $2<m<4$ or all $m_{i}=4$ for all $m_{i}$ (see
(2.15)). Now we define a special point $Q_{N}=(q_{1},...,q_{n})$ on
$\Gamma_{N}$, which satisfies

$\sum_{j=1}^{n}a_{ij}q_{j}=8\pi N\quad\forall i\in I.$

The second result is showing the tightness of the local masses.

######  Theorem 1.2.

Under the same assumptions in Theorem 1.1.

  * (1)

If $\rho^{k}\to\rho\in\Gamma_{N}(\rho\neq Q_{N})$ from $\mathcal{O}_{N-1}$ or
$\mathcal{O}_{N}$, then

(1.9) $|\sigma_{i,s}^{k}-\sigma_{i,t}^{k}|=O(\epsilon_{k}^{m-2}),\quad\text{
for }s\not=t,\,\,s,t\in\\{1,\cdots,N\\},\,\,i\in I.$

  * (2)

If $\rho^{k}\to Q_{N}$ from $\mathcal{O}_{N-1}$ or $\mathcal{O}_{N}$, then

(1.10)
$|\sigma_{i,s}^{k}-\sigma_{i,t}^{k}|=O(\epsilon_{k}^{2}\log\frac{1}{\epsilon_{k}}),\quad\text{
for }s\not=t,\,\,s,t\in\\{1,\cdots,N\\},\,\,i\in I.$

We remark here that the techniques developed in the proof of Theorems 1.1 and
1.2 also play the key roles in the study of the local uniqueness of the
bubbling solutions in Liouville systems[24]. For the sake of contradiction,
one also need to compare two sequence of the bubbling solutions at the blowup
points with the same limit $\rho$.

The leading terms of $\Lambda_{I}(\rho^{k})$ are different in two cases.
Before we state the results, we give some notations:

We define $N$ open sets $\Omega_{t,\delta_{0}}$ such that they are mutually
disjoint, each of them contains a bubbling disk and their union is $M$:

(1.11)
$B(p_{t}^{k},\delta_{0})\subset\Omega_{t,\delta_{0}},\quad\cup_{t=1}^{N}\overline{\Omega_{t,\delta_{0}}}=M,\quad\Omega_{t,\delta_{0}}\cap\Omega_{s,\delta_{0}}=\emptyset,\,\,\forall
t\neq s.$

Let

$I_{1}=\\{i\in I;\quad\lim_{k\to\infty}m_{i}^{k}=m.\quad\\}.$

and $G$ be the Green’s function defined by

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

and $\gamma$ is the regular part of the Green’s function. Note that in local
coordinates of a point, say $\eta$, $G$ is of the form

$G(x,\eta)=-\frac{1}{2\pi}\log|x-\eta|+\gamma(x,\eta).$

We also define

(1.12)
$G^{*}(p_{t}^{k},p_{s}^{k})=\left\\{\begin{array}[]{ll}\gamma(p_{t}^{k},p_{t}^{k}),\quad
s=t,\\\ G(p_{t}^{k},p_{s}^{k}),\quad s\neq t.\end{array}\right.\qquad
s,t=1,...,N.$

The third result is the leading terms of $\Lambda_{I}(\rho^{k})$ of the first
case.

######  Theorem 1.3.

Under the same assumptions in Theorem 1.1. If
$\rho^{k}\to\rho\in\Gamma_{N}(\rho\neq Q_{N})$ from $\mathcal{O}_{N-1}$ or
$\mathcal{O}_{N}$, then

(1.13) $\Lambda_{I}(\rho^{k})=(D+o(1))\frac{\epsilon_{k}^{m-2}}{N}.$

Here, The quantity $D$ is defined as follows

(1.14) $D=\sum_{i\in
I_{1}}e^{D_{i}-\alpha_{i}}\sum_{t=1}^{N}c_{t}\lim_{\delta_{0}\to
0}\bigg{(}\delta_{0}^{2-m}-\frac{(m-2)}{2\pi}\int_{\hat{\Omega}_{t,\delta_{0}}}(\frac{h_{i}^{k}(x)}{h_{i}^{k}(p_{t}^{k})}e^{2\pi
m\sum_{l=1}^{N}(G(x,p_{l}^{k})-G^{*}(p_{t}^{k},p_{l}^{k}))}dV_{g}\bigg{)}.$

where $\hat{\Omega}_{t,\delta_{0}}=\Omega_{t,\delta_{0}}\setminus
B(p_{t}^{k},\delta_{0})$, $c_{t}=\frac{h_{i}^{k}(p_{t}^{k})e^{2\pi
m\sum_{l=1}^{N}G^{*}(p_{t}^{k},p_{l}^{k})}}{h_{i}^{k}(p_{1}^{k})e^{2\pi
m\sum_{l=1}^{N}G^{*}(p_{1}^{k},p_{l}^{k})}}$ and $D_{i}$, $\alpha_{i}$ are
constants defined in (3.14).

We remark here that $\Lambda_{I}(\rho^{k})>0$ for
$\rho^{k}\in\mathcal{O}_{N-1}$ and $\Lambda_{I}(\rho^{k})<0$ for
$\rho^{k}\in\mathcal{O}_{N}$. Thus, if $D\not=0$, then the blowup solutions
with $\rho^{k}$ occur as $\rho^{k}\to\rho\in\Gamma_{N}(\rho\not=Q_{N})$ only
from one side of $\Gamma_{N}$. Furthermore, it yields a uniform bound of
solutions as $\rho^{k}$ converges to $\Gamma_{N}(\rho\not=Q_{N})$ from
$\mathcal{O}_{N-1}$ provided that $D<0$. It is easy to see that for fixed $k$
the following limit exists:

(1.15) $\lim_{\delta_{0}\to
0}\bigg{(}\delta_{0}^{2-m}-\frac{m-2}{2\pi}\int_{\hat{\Omega}_{t,\delta_{0}}}\frac{h_{i}^{k}(x)}{h_{i}^{k}(p_{t}^{k})}e^{2\pi
m(\sum_{l=1}^{N}G(x,p_{l}^{k})-G^{*}(p_{t}^{k},p_{l}^{k}))}dV_{g}\bigg{)}$

because the leading terms from the Green’s function is
$-\frac{1}{2\pi}\log|x-p_{t}^{k}|$. The study of the sign of $D$ is another
interesting fundamental question. However, it is out of the scope of of the
present article. We will come back to this issue in the further study.

The fourth major result is concerned with the leading term of
$\Lambda_{I}(\rho^{k})$ when $\rho^{k}\to Q_{N}$ and the major difference is
that the leading term only depends on the curvature and coefficient functions
at blowup points: $p_{1}^{k}$, …, $p_{N}^{k}$:

######  Theorem 1.4.

Under the same assumptions in Theorem 1.1. If $\rho^{k}\to Q_{N}$ from
$\mathcal{O}_{N-1}$ or $\mathcal{O}_{N}$, then

(1.16)
$\Lambda_{I}(\rho^{k})=-4\sum_{i=1}^{n}\sum_{t=1}^{N}b_{it}^{k}\epsilon_{k}^{2}\log\epsilon_{k}^{-1}+O(\epsilon_{k}^{2})$

where

(1.17) $\displaystyle b_{it}^{k}$
$\displaystyle=e^{D_{i}-\alpha_{i}}\bigg{(}\frac{1}{4}\frac{\Delta
h_{i}^{k}(p_{t}^{k})}{h_{i}^{k}(p_{t}^{k})}-K(p_{t}^{k})+4\pi N$
$\displaystyle+4\pi\frac{\nabla
h_{i}^{k}(p_{t}^{k})}{h_{i}^{k}(p_{t}^{k})}\cdot\sum_{l=1}^{N}\nabla_{1}G^{*}(p_{t}^{k},p_{l}^{k})+16\pi^{2}|\sum_{l=1}^{N}\nabla_{1}G^{*}(p_{t}^{k},p_{l}^{k})|^{2}\bigg{)}.$

and $K$ is the Gaussian curvature.

As an application of the formula (1.16), we obtain the uniform bound of
solutions as $\rho^{k}$ converges to $Q_{N}$ from $\mathcal{O}_{N-1}$ provided
$b_{it}^{k}>0$ for $t=1,\cdots,N.$

The fifth result is about the locations of the blowup points and the mutual
relation of coefficient functions.

######  Theorem 1.5.

If $\rho^{k}\to\rho\neq Q$ as in (1.4), then for $t=1,...,N$

(1.18) $\sum_{i=1}^{n}\big{(}\nabla(\log h_{i}^{k})(p_{t}^{k})+2\pi
m_{i}\sum_{s=1}^{N}\nabla_{1}G^{*}(p_{t}^{k},p_{s}^{k})\big{)}=O(\epsilon_{k}^{m-2}),\,\,t=1,..,N$

where $\nabla_{1}$ means the differentiation with respect to the first
component. If $\rho^{k}\to Q$ as in (1.4),

(1.19) $\sum_{i=1}^{n}\big{(}\nabla(\log
h_{i}^{k})(p_{t}^{k})+8\pi\sum_{s=1}^{N}\nabla_{1}G^{*}(p_{t}^{k},p_{s}^{k})\big{)}=O(\epsilon_{k}^{2}\log\epsilon_{k}^{-1}),\,\,t=1,..,N$

Besides the location of blowup points, we also reveal new information about
coefficient functions in the next theorem. For convenience we set

(1.20) $H_{i,t}^{k}:=\log h_{i}^{k}(p_{t}^{k})+2\pi
m_{i}\sum_{l=1}^{N}G^{*}(p_{t}^{k},p_{l}^{k}).$

######  Theorem 1.6.

Let $H_{i,t}^{k}$ be defined as in (1.20), then we have

(1.21) $H_{i,t}^{k}=H_{i,s}^{k}+O(\epsilon_{k}^{m-2}),\quad\forall i\in
I,\quad\forall t\neq s,\quad\mbox{if}\quad m<4,$ (1.22)
$H_{i,t}^{k}=H_{i,s}^{k}+O(\epsilon_{k}^{2}\log\frac{1}{\epsilon_{k}}),\quad\forall
i\in I,\quad t\neq s,\quad\mbox{if}\quad m=4.$

The tightness of the coefficient functions is not seen in the case of $N=1$.
This show that the occurrence of multi bubbles forces the coefficient
functions at the blowup points to be the same. In construction of bubbling
solutions one needs to know the precise information about bubbling
interactions, exact location of blowup points, accurate vanishing rate of
coefficient functions and specific leading terms in asymptotic expansions. All
these have been covered in the main results of this article. Until now the
construction of bubbling solutions for Liouville systems is still in the early
stage of development, as the constructions so far are still restricted to
single blowup point situations [22].

In all these main results the readers can see that sharp estimates are
obtained for all the error terms. This is why we think these results will play
as a central role in applications. We expect the theorems of this article to
serve as a benchmark for more sophisticated discussions in the near future.
Also Gu-Zhang [21] completed a degree counting program for singular Liouville
systems, the corresponding discussion for leading terms of approximating
critical hyper-surfaces of singular Liouville systems is another exciting
unconquered land to explore. Besides these immediate impacts to closely
related fields, the idea of the proof in this article, the way to overcome
major difficulties in bubble interaction could lead to major advance in Chern-
Simons type equations (see [23]) as well.

The main difficulty in the proof of the main theorems is on the interaction of
bubbling solutions, which has a large to do with the nature of Liouville
systems. For global solutions defined in $\mathbb{R}^{2}$, it is established
in [15, 32] that total integrations of all components form a $n-1$ dimensional
hyper-surface similar to $\Gamma_{1}$. This continuum of energy brings great
difficulty for bubbling interaction: If we only have two bubbling disks, the
energy in each disk is very close to $\Gamma_{1}$, but to identify the leading
term of $\rho^{k}-\rho\in\Gamma_{2}$ one has to prove that they are both
tending the _same point_ on $\Gamma_{1}$. It does not help to use the fact
that they have the same limit, because the energy sequence may tend to its
limit position very slowly. This difficulty does not exist for Toda systems,
because the energy set for Toda systems is discrete. When we have only one
bubble, this bubble interaction situation can be avoided ( see [34]). So the
main contribution in this article is to prove that the bubbling solutions have
almost the same energy in each bubbling disk. The key idea is as follows:
Around each blowup point, we first use an approximation theorem of Lin-Zhang
[34] to have an initial expansion of the bubbling solution. The first term of
this expansion is a sequence of global solutions. In this article we identify
what the global sequence is around each bubbling point and compare them using
a key idea of Lin-Zhang [32] in their proof of the classification theorem for
Liouville systems. It turns out that after scaling, the global sequences are
extremely close to one another. Then we further prove that the energy of the
global sequence is not too much different from the bubbling solutions in each
bubbling disk. All the error terms must be carefully identified in order to
single out the leading terms in the main theorems.

The organization of this article is as follows: In section two we deploy the
basic setting for all the topics in this article and invoke a few
approximation results in previous works of Lin-Zhang. In section three we
prove the closeness of bubbling solutions around different blowup points, in
which Theorems 1.1, 1.2 and 1.6 will be proved. The proof is set-up in two
stages as the approximation becomes better in the second stage. Finally in
section four all other main theorems will be proved based on the precise
estimates in section three.

##  2\. Approximation around a blowup point

First we claim that we can assume $u^{k}=(u_{1}^{k},...,u_{n}^{k})$ to satisfy

(2.1) $\int_{M}h_{i}^{k}e^{u_{i}^{k}}dV_{g}=1,\quad Vol(M)=1.$

because otherwise we just consider

(2.2) $\Theta_{i}^{k}=u_{i}^{k}-\log\int_{M}h_{i}^{k}e^{u_{i}^{k}}dV_{g},\quad
i\in I.$

Then we have

(2.3)
$-\Delta_{g}\Theta_{i}^{k}=\sum_{j=1}^{n}a_{ij}\rho_{j}^{k}(h_{j}^{k}e^{\Theta_{j}^{k}}-1).$

where $\Theta^{k}$ satisfies (2.1). Here we first set up preliminary
discussions about the profile of $u^{k}$ near a blowup point. Suppose $p$ is a
blowup point and in $B(p,\delta)$ there is only one blowup point of $u^{k}$.
Let

$\tilde{M}_{k}=\max_{i\in I,x\in
B(p,\delta)}u_{i}^{k}(x)+\log(\rho_{i}^{k}h_{i}^{k}(\tilde{p}_{k}))\quad\mbox{and
}\quad\tilde{\epsilon}_{k}=e^{-\frac{1}{2}\tilde{M}_{k}},$

and $\tilde{p}_{k}$ be where $\tilde{M}_{k}$ is attained ( $\tilde{p}_{k}\to
p$). Then the functions

$\tilde{v}_{i}^{k}(y)=u_{i}^{k}(\tilde{p}_{k}+\tilde{\epsilon}_{k}y)+2\log\tilde{\epsilon}_{k}$

converge in $C^{2}_{loc}(\mathbb{R}^{2})$ to the limit function
$v=(v_{1},...,v_{n})$ which is a global solution of

(2.4) $\left\\{\begin{array}[]{ll}-\Delta
v_{i}=\sum_{j}a_{ij}e^{v_{j}},\quad\mathbb{R}^{2},\quad i=1,..,n\\\ \\\
\int_{\mathbb{R}^{2}}e^{v_{i}}<\infty,\quad i=1,..,n,\quad\max_{i\in
I}v_{i}(0)=0.\end{array}\right.$

Here we note that it is established in [33] that with assumptions (H1), (H2)
all the bubbling solutions are fully bubbling: the limit must have $n$
equations and no component is lost in the limiting taking process. The
classification of all global solutions of (2.4) was completed in the work of
Chipot-Shafrir-Wolansky [15] and Lin-Zhang [32]. All components of
$v=(v_{1},...,v_{n})$ have one common point of symmetric symmetry. In this
context, this common point is the origin.

To state more precise approximation results we write the equation in local
coordinates around $\tilde{p}_{k}$. In this coordinate $ds^{2}$ has the form
$e^{\phi(y_{\tilde{p}_{k}})}(dy_{1}^{2}+dy_{2}^{2})$ where

(2.5)
$|\nabla\phi(0)|=0,\quad\phi(0)=0,\quad\Delta_{y_{\tilde{p}_{k}}}\phi=-2Ke^{\phi}$

where $K$ is the Gauss curvature. In local coordinates, (2.2) becomes

(2.6) $-\Delta
u_{i}^{k}=\sum_{j=1}^{n}a_{ij}\rho_{j}^{k}e^{\phi}(h_{j}^{k}e^{u_{j}^{k}}-1),\quad\mbox{in}\quad
B(0,\delta).$

Let $f_{i}^{k}$ be defined by

(2.7) $\left\\{\begin{array}[]{ll}\Delta
f_{i}^{k}=\sum_{j=1}^{n}a_{ij}\rho_{j}^{k}e^{\phi},\quad\mbox{in}\quad
B(0,\delta).\\\ \\\ f_{i}^{k}(0)=0,\quad\nabla
f_{i}^{k}(0)=0.\end{array}\right.$

Then we have

(2.8)
$-\Delta(u_{i}^{k}-f_{i}^{k})=\sum_{j=1}^{n}a_{ij}\rho_{j}^{k}h_{j}^{k}e^{u_{j}^{k}-f_{j}^{k}}e^{f_{j}^{k}}e^{\phi},\quad\mbox{in}\quad
B(0,\delta).$

If we set

$\tilde{h}_{i}^{k}(x)=\frac{h_{i}^{k}(x)}{h_{i}^{k}(\tilde{p}_{k})}e^{\phi+f_{i}^{k}},$

we have $\tilde{h}_{i}^{k}(0)=1$,

(2.9) $\nabla(\log\tilde{h}_{i}^{k})(0)=\nabla(\log
h_{i}^{k})(\tilde{p}_{k}).$ (2.10)
$\Delta(\log\tilde{h}_{i}^{k})(0)=\Delta(\log
h_{i}^{k})(\tilde{p}_{k})-2K(\tilde{p}_{k})+\sum_{j=1}^{n}a_{ij}\rho_{j}^{k}.$

Thus we set

(2.11) $\tilde{u}_{i}^{k}=u_{i}^{k}+\log\rho_{i}^{k}+\log
h_{i}^{k}(\tilde{p}_{k})-f_{i}^{k},$

and write the equation of $\tilde{u}_{i}^{k}$ as

(2.12)
$\Delta\tilde{u}_{i}^{k}+\sum_{j=1}^{n}a_{ij}\tilde{h}_{j}^{k}e^{\tilde{u}_{j}^{k}}=0,\quad\mbox{
in}\quad B(0,\delta).$

Now we introduce $\phi_{i}^{k}$ to be a harmonic function defined by the
oscillation of $\tilde{u}_{i}^{k}$ on $B(\tilde{p}_{k},\delta)$:

(2.13) $\left\\{\begin{array}[]{ll}-\Delta\phi_{i}^{k}=0,\quad\mbox{in}\quad
B(0,\delta),\\\ \\\
\phi_{i}^{k}=\tilde{u}_{i}^{k}-\frac{1}{2\pi\delta}\int_{\partial
B(0,\delta)}\tilde{u}_{i}^{k},\quad\mbox{on}\quad\partial
B(0,\delta).\end{array}\right.$

Obviously $\phi_{i}^{k}(0)=0$ by the mean value theorem and $\phi_{i}^{k}$ is
uniformly bounded on $B(0,\delta/2)$ because $u_{i}^{k}$ has finite
oscillation away from blowup points. It is a standard fact (see [32, 33]) that
the location of $\max_{i\in I}\max_{x\in
B(0,\delta)}\tilde{u}_{i}^{k}(x)-\phi_{i}^{k}(x)$ is $O(\epsilon_{k}^{2})$
(roughly speaking, the reason is $0$ is a non-degenerate maximum of
$u_{i}^{k}$ and $\phi_{i}^{k}(0)=0$).

Going back to the original coordinate system, we call the maximum point after
perturbation $p_{k}$. Now we set

$M_{k}=\max_{i\in I}\max_{x\in
B(0,\delta)}u_{i}^{k}(x)+\log(\rho_{i}^{k}h_{i}^{k}(p_{k}))-\phi_{i}^{k}(x),\quad\epsilon_{k}=e^{-\frac{1}{2}M_{k}},$

and we let $V^{k}=(V_{1}^{k},..,V_{n}^{k})$ be the radial solutions of

(2.14) $\left\\{\begin{array}[]{ll}-\Delta
V_{i}^{k}=\sum_{j=1}^{n}a_{ij}e^{V_{j}^{k}}\quad\mbox{in}\quad\mathbb{R}^{2},\quad
i\in I,\\\ \\\
V_{i}^{k}(0)=u_{i}^{k}(p_{k})+\log(\rho_{i}^{k}h_{i}^{k}(p_{k}))-\phi_{i}^{k}(p_{k}),\quad
i\in I.\end{array}\right.$

Note that since $p_{k}=\tilde{p}_{k}+O(\epsilon_{k}^{2})$, it is easy to
obtain that the oscillation of $V_{i}^{k}$ on $\partial B(p_{k},\delta)$ is
$O(\epsilon_{k}^{2})$. The sequence of function
$V^{k}=(V_{1}^{k},...,V_{n}^{k})$, which agrees with
$u_{i}^{k}(x)+\log(\rho_{i}^{k}h_{i}^{k}(p_{k}))-\phi_{i}^{k}(x)$ at $p_{k}$,
gives the first term in the approximation of $u_{i}^{k}$ near $p$. To state
more precise approximation terms, we use

$v_{i}^{k}(y)=u_{i}^{k}(p_{k}+\epsilon_{k}y)+\log(\rho_{i}^{k}h_{i}^{k}(p_{k}))-\phi_{i}^{k}(\epsilon_{k}y)+2\log\epsilon_{k},\quad|y|<\frac{\delta}{2}\epsilon_{k}^{-1}$

and the following rough approximation theorem is established in [34]:

######  Remark 2.1.

The notations $M_{k}$, $\epsilon_{k}$ are the same as those in the
introduction. It is confusing at this moment, later in the multiple bubbling
situation we will use $M_{k,t}$ and $\epsilon_{k,t}$ to denote the maximum of
bubbling solutions and decay rate in each bubbling disk
$B(p_{t}^{k},\delta_{0})$. $M_{k}$ is the maximum of $M_{k,t}$. But our
analysis will show that we can replace $M_{k}$ by $M_{k,t}$, $\epsilon_{k}$ by
$\epsilon_{k,t}$ for any $t$ and the nature of the proof does not change.

Before citing the approximation theorems in [34] we mention one simple fact
implied by the Pohozaev identity. Let
$\sigma_{i}^{k}=\frac{1}{2\pi}\int_{B(p_{k},\delta)}h_{i}^{k}e^{u_{i}^{k}}$
and $m_{i}^{k}=\sum_{j=1}^{n}a_{ij}\sigma_{j}^{k}$. Let
$\sigma_{i}=\lim_{k\to\infty}\sigma_{i}^{k}$ and
$m_{i}=\lim_{k\to\infty}m_{i}^{k}$. As usual we set
$m=\min\\{m_{1},...,m_{n}\\}$. Then it is established in [32] that each
$m_{i}>2$ and

$\sum_{i=1}^{n}\sigma_{i}(m_{i}-4)=0.$

Since each $\sigma_{i}>0$ it is easy to see that either

(2.15) $m<4,\quad\mbox{or}\quad m_{i}=4\,\,\forall i\in I.$

The first approximation theorem established in [34] is a rough one that does
not distinguish $m<4$ or $m=4$.

######  Theorem 2.1.

Given $\delta>0$, there exist $C(\delta)>0$, $k_{0}(\delta)>1$ such that for
$|y|\leq\frac{\delta}{2}\epsilon_{k}^{-1}$ and $|\alpha|=0,1$, the following
holds for all $k\geq k_{0}$

(2.16)
$|D^{\alpha}(v_{i}^{k}(y)-V_{i}^{k}(\epsilon_{k}y)-2\log\epsilon_{k}-\Phi_{i}^{k}(y))|\leq
C\epsilon_{k}^{2}(1+|y|)^{4-m-|\alpha|+\delta}.$

where

$\Phi_{i}^{k}(y)=\epsilon_{k}(G_{1,i}^{k}(r)\cos\theta+G_{2,i}^{k}(r)\sin\theta)$

with

(2.17) $|G_{t,i}^{k}(r)|\leq Cr(1+r)^{2-m+\delta}\quad t=1,2.$

Note that $\Phi^{k}=(\Phi_{1}^{k},..,\Phi_{n}^{k})$ denotes the projection of
$v_{i}^{k}$ onto $span\\{\sin\theta,\cos\theta\\}$. i.e.

(2.18)
$\Phi^{k}_{i}(r\cos\theta,r\sin\theta)=\epsilon_{k}(G_{1,i}^{k}(r)\cos\theta+G_{2,i}^{k}(r)\sin\theta),\quad
i\in I$

with $G_{t,i}^{k}(r)$ ($t=1,2$) satisfying some ordinary differential
equations to be specified later. The estimate for $|\alpha|=1$ follows from
standard gradient estimate for elliptic equations.

Theorem 2.1 does not distinguish $m<4$ or $m=4$. But using Theorem 2.1 in more
careful computation for $m<4$ and $m=4$ gives rise to more accurate results as
follows: Here it is important to observe that $m-2<2$ if $m<4$. If we use

$m_{i}^{k}=\frac{1}{2\pi}\int_{B(p,\delta)}\rho_{i}^{k}h_{i}^{k}(x)e^{u_{i}^{k}}dV_{g},\quad
m_{k}=\min\\{m_{1}^{k},...,m_{n}^{k}\\}.$

Clearly $m_{k}\to m\in(2,4)$.

######  Theorem 2.2.

Suppose $m<4$, then for $|y|\leq\frac{\delta_{0}}{2}\epsilon_{k}^{-1}$ and
$i\in I$

$\displaystyle|D^{\alpha}\big{(}v_{i}^{k}(y)-(V_{i}^{k}(\epsilon_{k}y)+2\log\epsilon_{k})-\Phi^{k}_{i}(y)\big{)}|$
$\displaystyle\leq$ $\displaystyle
C\epsilon_{k}^{2}(1+|y|)^{4-m_{k}-|\alpha|}\log(2+|y|).\quad|\alpha|=0,1$

Moreover $G_{t,i}^{k}$ ($t=1,2,i\in I$) satisfy

(2.20) $|G_{t,i}^{k}(r)|\leq Cr(1+r)^{2-m_{k}}\quad 0<r<\epsilon_{k}^{-1}.$

Note that Theorem 2.2 is slightly stronger than Theorem 4.2 of [34] because
the latter has a logarithmic term. The reason is in the context of Theorem
2.2, the function $v_{i}^{k}$ agrees with its approximation at the origin.
Theorem 2.2 can be proved just like Theorem 4.2 in [34].

######  Theorem 2.3.

If $m=4$ and $|m_{i}^{k}-4|\leq C/\log\epsilon_{k}^{-1}$ for all $i\in I$,
then we have, for $|y|\leq\frac{\delta_{0}}{2}\epsilon_{k}^{-1}$ and $i\in I$

$\displaystyle|D^{\alpha}(v_{i}^{k}(y)-(V_{i}^{k}(\epsilon_{k}y)+2\log\epsilon_{k})-\Phi_{i}^{k}(y)\big{)}|$
$\displaystyle\leq$ $\displaystyle
C\epsilon_{k}^{2}(1+|y|)^{-|\alpha|}(\log(2+|y|))^{2}.\quad|\alpha|=0,1,$

where $\Phi^{k}$ is of the form stated in (2.18) with $G_{t,i}^{k}$ ($t=1,2$)
satisfying

(2.22) $|G_{t,i}^{k}(r)|\leq Cr(1+r)^{-2},\quad 0<r<\epsilon_{k}^{-1},\quad
i\in I.$

##  3\. Rough estimates about bubbling magnitudes

In this section, we will prove Theorems 1.1, 1.2 and 1.6. First in this
section for simplicity we assume there are only two blowup points $p$ and $q$.
The nature of analysis does not change if we have more blowup points.

Now we use Green’s representation to describe the neighborhood of $p_{k}$. The
expression of $u_{i}^{k}$ is

$\displaystyle u_{i}^{k}(x)$
$\displaystyle=\bar{u}_{i}^{k}+\int_{M}G(x,\eta)\sum_{j=1}^{n}a_{ij}\rho_{j}^{k}h_{j}^{k}e^{u_{j}^{k}}dV_{g}$
$\displaystyle=\bar{u}_{i}^{k}+\left(\int_{B(p_{k},\delta)}+\int_{B(q_{k},\delta)}+\int_{M\setminus(B(p_{k},\delta)\cup
B(q_{k},\delta))}\right)G(x,\eta)\sum_{j=1}^{n}a_{ij}\rho_{j}^{k}h_{j}^{k}e^{u_{j}^{k}}dV_{g}$
$\displaystyle=\bar{u}_{i}^{k}+I+II+III.$

Here $\bar{u}_{i}^{k}=\int_{M}u_{i}^{k}dV_{g}$ and we use
$E=O(\epsilon_{k}^{m-2+\delta})$ to denote a very rough order of the error. If
we use Theorem 2.1 in the evaluation, after cancellation we have

$\displaystyle I$
$\displaystyle=\int_{B_{\delta}}(-\frac{1}{2\pi}\log|\eta|+\gamma(x,\eta))\sum_{j=1}^{n}a_{ij}\tilde{h}_{j}^{k}e^{\tilde{u}_{j}^{k}}d\eta$
$\displaystyle=-m_{i}^{k}\log|x|+2\pi m_{i}^{k}\gamma(x,p_{k})+E.$

where $m_{i}^{k}=\frac{1}{2\pi}\sum_{j=1}^{n}a_{ij}\sigma_{j}^{k}$ and
$\sigma_{i}^{k}=\int_{B(p,\delta)}h_{i}^{k}e^{u_{i}^{k}}dV_{g}$, we use
$\bar{\sigma}_{i}^{k}$, $\bar{m}_{i}^{k}$ to denote the integrations in
$B(q_{k},\delta)$. Here we note that in the evaluation of the integrals, the
$\Phi_{i}^{k}$ terms disappear because of the symmetry of the domain.
Similarly the integration around $q_{k}$ gives

$II=2\pi\bar{m}_{i}^{k}G(x,q_{k})+E,\quad III=E.$

Thus in the neighborhood of $p_{k}$ we have

(3.1) $u_{i}^{k}(x)=\bar{u}_{i}^{k}-m_{i}^{k}\log|x|+2\pi
m_{i}^{k}\gamma(x,p_{k})+2\pi\bar{m}_{i}^{k}G(x,q_{k})+E$

in, say $B(p_{k},2\delta)\setminus B(p_{k},\delta/2)$. If we use the
approximation theorems to evaluate $u_{i}^{k}$ at $p_{t}^{k}$, it is easy to
obtain

(3.2) $\bar{u}_{i}^{k}=(1-\frac{m_{i}^{k}}{2})M_{k}+O(1).$

Before more advanced estimates we first establish an elementary one:

Proof of Theorem 1.1: From (3.2) we have

(3.3)
$(1-\frac{m_{i}^{k}}{2})M_{k}=(1-\frac{\bar{m}_{i}^{k}}{2})\bar{M}_{k}+O(1).$

Let

$\lambda_{k}=M_{k}/\bar{M}_{k},\quad\delta_{i}^{k}=O(1)/\bar{M}_{k},$

we have

(3.4)
$(\frac{m_{i}^{k}}{2}-1)\lambda_{k}+\delta_{k}=(\frac{\bar{m}_{i}^{k}}{2}-1)+o(\epsilon_{k}^{\delta}),\quad\mbox{for
some }\delta>0.$

It is established in [32, 34] that $m^{k}=(m_{1}^{k},...,m_{n}^{k})$ satisfies

(3.5)
$\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}(\frac{m_{i}^{k}-2}{2})(\frac{m_{j}^{k}-2}{2})=\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}+E.$

and (3.5) also holds for $\bar{m}^{k}=(\bar{m}_{1}^{k},...,\bar{m}_{n}^{k})$.
Thus using (3.4) in (3.5) for $\bar{m}_{k}$, we have

$\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}(\frac{m_{i}^{k}-2}{2}\lambda_{k}+\delta_{i}^{k})(\frac{m_{j}^{k}-2}{2}\lambda_{k}+\delta_{j}^{k})=\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}+E,$

which can be written as a quadratic expression of $\lambda_{k}$:

$\displaystyle\lambda_{k}^{2}\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}(\frac{m_{i}^{k}-2}{2})(\frac{m_{j}^{k}-2}{2})+2\lambda_{k}\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}(\frac{m_{i}^{k}-2}{2})\delta_{j}^{k}+\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}\delta_{i}^{k}\delta_{j}^{k}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}+E/\bar{M}_{k}.$

Let

$B_{k}=\frac{\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}(m_{i}^{k}-2)\delta_{j}^{k}}{\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}},\quad
C_{k}=\frac{\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}\delta_{i}^{k}\delta_{j}^{k}}{\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}}.$

Here we note that $\sum_{i=1}^{n}\sum_{j=1}^{n}a^{ij}>0$ because $(H2)$
requires $A$ to be invertible and $\sum_{j=1}^{n}a^{ij}\geq 0$ for all $i$.
Then

$\lambda_{k}^{2}+B_{k}\lambda_{k}+C_{k}=1+E/\bar{M}_{k}.$

First it is obvious to observe that $\lim_{k\to\infty}\lambda_{k}=1$. Thus by
$C_{k}=O(\bar{M}_{k}^{-2})$ and $B_{k}=O(\bar{M}_{k}^{-1})$

$\lambda_{k}=\sqrt{1-C_{k}+\frac{B_{k}^{2}}{4}}-\frac{B_{k}}{2}=1-\frac{B_{k}}{2}-\frac{1}{2}(C_{k}-\frac{B_{k}^{2}}{4})+O(\bar{M}_{k}^{-4}).$

This verifies that $M_{k}-\bar{M}_{k}=O(1)$ which justifies
$O(\epsilon_{k}^{m_{k}-2})=O(\bar{\epsilon}_{k}^{\bar{m}_{k}-2})$. Theorem 1.1
is established. $\Box$

By Theorem 1.1, all the error terms above can be improved to
$E=O(\epsilon_{k}^{m_{k}-2})$. Note that it is not $O(\epsilon_{k}^{m-2})$
yet, because the closeness of $m_{k}$ and $m$ is not derived yet. Another
consequence of Theorem 1.1 is that

(3.6)
$\frac{\sum_{j=1}^{n}a_{ij}\rho_{i}^{k}}{2\pi}-m_{i}^{k}-\bar{m}_{i}^{k}=O(\epsilon_{k}^{m_{k}-2}).$

Indeed, integrating the equation for $u_{i}^{k}$ (which is (1.1)), we have

$\sum_{j=1}^{n}\int_{M}a_{ij}\rho_{j}^{k}h_{j}^{k}e^{u_{j}^{k}}=\sum_{j=1}^{n}a_{ij}\rho_{j}^{k}.$

The integration of the left in $B(p,\delta)$ and $B(q,\delta)$ gives

$m_{i}^{k}+\bar{m}_{i}^{k}+O(\epsilon_{k}^{m_{k}-2})=\sum_{j=1}^{n}a_{ij}\rho_{j}^{k}/(2\pi).$

Thus (3.6) is verified.

Here we recall a theorem in [34] about location of blowup points:

Let $p_{t}^{k}$ be blowup points described as before. Then at each blowup
point $p_{t}^{k}$, let $\phi_{it}^{k}$ be the harmonic function that
eliminates the oscillation of $\tilde{u}_{i}^{k}$ on $\partial
B(\tilde{p}_{t}^{k},\delta)$ for $\delta>0$ small. Then it is proved in [34]
that

######  Theorem 3.1.

If $m<4$

(3.7) $|\sum_{i=1}^{n}\bigg{(}\partial_{l}(\log
h_{i}^{k})(p_{t}^{k})+\partial_{l}\phi_{it}^{k}(p_{t}^{k})\bigg{)}\rho_{it}^{k}|\leq
C\epsilon_{k}^{m_{k}-2},\quad l=1,2,$

where $C$ is independent of $k$. On the other hand, if $m=4$

(3.8) $|\sum_{i=1}^{n}\bigg{(}\partial_{l}(\log
h_{i}^{k})(p_{t}^{k})+\partial_{l}\phi_{i}^{k}(p_{t}^{k})\bigg{)}\rho_{it}^{k}|\leq
C\epsilon_{k}^{2}\log\epsilon_{k}^{-1},\quad l=1,2,$

where
$\rho_{it}^{k}=\int_{B(p_{t}^{k},\delta)}\rho_{i}^{k}h_{i}^{k}e^{u_{i}^{k}}dV_{g}$.

In $B(p_{k},\delta)$, by the definition of $\tilde{u}_{i}^{k}$ in (2.11) and
the estimate of $u_{i}^{k}$ in (3.1) we now have

(3.9) $\displaystyle\tilde{u}_{i}^{k}(x)=\bar{u}_{i}^{k}-m_{i}^{k}\log|x|+2\pi
m_{i}^{k}\gamma(x,p_{k})+2\pi\bar{m}_{i}^{k}G(x,q_{k})$
$\displaystyle+\log\rho_{i}^{k}+\log
h_{i}^{k}(p_{k})-f_{i}^{k}(x)+O(\epsilon_{k}^{m_{k}-2}).$

In this neighborhood, $\tilde{u}_{i}^{k}$ is of the form

$\tilde{u}_{i}^{k}(x)=V_{i}^{k}(x)+\phi_{i}^{k}(x)+O(\epsilon_{k}^{m_{k}-2}),\quad
x\in B(0,\delta/2)\setminus B(0,\delta/8)$

where $\phi_{i}^{k}$ be the harmonic function on $B(p_{k},\delta)$ that
eliminates the oscillation of $\tilde{u}_{i}^{k}$:

$\Delta\phi_{i}^{k}(x)=0,\quad\mbox{in}\quad
B(p_{k},\delta),\quad\phi_{i}^{k}(x)=\tilde{u}_{i}^{k}(x)-\frac{1}{2\pi\delta}\int_{\partial
B(p_{k},\delta)}\tilde{u}_{i}^{k}.$

Note that we have used the fact that the first order terms
$\Phi_{i}^{k}(x)=O(\epsilon_{k}^{m_{k}-2})$ when $x\sim 1$. Notice that by
(3.6) and the fact that $q_{k}$ is not in $B(p_{k},\delta)$,

$\displaystyle\Delta_{g}(2\pi
m_{i}^{k}\gamma(x,p_{k})+2\pi\bar{m}_{i}^{k}G(x,q_{k})-f_{i}^{k})$
$\displaystyle=$ $\displaystyle 2\pi
m_{i}^{k}+2\pi\bar{m}_{i}^{k}-\sum_{j=1}^{n}a_{ij}\rho_{j}^{k}=O(\epsilon_{k}^{m_{k}-2})$

Thus

(3.10) $\Delta(2\pi
m_{i}^{k}\gamma(x,p_{k})+2\pi\bar{m}_{i}^{k}G(x,q_{k})-f_{i}^{k})=O(\epsilon_{k}^{m_{k}-2}).$

By the definition of $\tilde{u}_{i}^{k}$ in (2.11) and the comparison of (3.1)
and (3.10) we have

(3.11) $\displaystyle\phi_{i}^{k}(x)=2\pi
m_{i}^{k}(\gamma(x,p_{k})-\gamma(p_{k},p_{k}))+2\pi\bar{m}_{i}^{k}(G(x,q_{k})-G(p_{k},q_{k}))$
$\displaystyle-f_{i}^{k}+O(\epsilon_{k}^{m_{k}-2}).$

So if we rewrite the expression of $\tilde{u}_{i}^{k}$ as

$\displaystyle\tilde{u}_{i}^{k}(x)$
$\displaystyle=\bar{u}_{i}^{k}-m_{i}^{k}\log|x|+\phi_{i}^{k}(x)$
$\displaystyle+2\pi
m_{i}^{k}\gamma(p_{k},p_{k})+2\pi\bar{m}_{i}^{k}(G(p_{k},q_{k}))+\log(\rho_{i}^{k}h_{i}^{k}(p_{k}))+O(\epsilon_{k}^{m_{k}-2}).$

then we see that for $x\in B(0,\delta)\setminus B(0,\delta/8)$,

(3.12) $\displaystyle V_{i}^{k}(x)=-m_{i}^{k}\log|x|+\bar{u}_{i}^{k}$
$\displaystyle+2\pi
m_{i}^{k}\gamma(p_{k},p_{k})+2\pi\bar{m}_{i}^{k}G(p_{k},q_{k})+\log(\rho_{i}^{k}h_{i}^{k}(p_{k}))+O(\epsilon_{k}^{m_{k}-2}).$

Similarly around $q_{k}$ we have

(3.13)
$\displaystyle\tilde{V}_{i}^{k}(x)=-\bar{m}_{i}^{k}\log|x|+\bar{u}_{i}^{k}$
$\displaystyle+2\pi\bar{m}_{i}^{k}\gamma(q_{k},q_{k})+2\pi
m_{i}^{k}G(q_{k},p_{k})+\log(\rho_{i}^{k}h_{i}^{k}(q_{k}))+O(\epsilon_{k}^{m_{k}-2}).$

Let

$M_{k}=\max_{i\in I}\max_{x}\tilde{u}_{i}^{k}(x),\quad\mbox{in}\quad
B(p_{k},\delta),$

and $\bar{M}_{k}$ be the corresponding maximum in $B(q_{k},\delta)$. Then it
is proved in [33] that

$M_{k}-\tilde{u}_{i}^{k}(p_{k})=O(1).$

We shall use the following notations:

$D_{i}^{k}=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}\sum_{j=1}^{n}a_{ij}e^{V_{j}^{k}},\quad\alpha_{i}^{k}=M_{k}-\tilde{u}_{i}^{k}(p_{k}).$

$\bar{D}_{i}^{k}$, $\bar{\alpha}_{i}^{k}$ can be understood similarly.

In order to obtain accurate estimate of $|m_{i}^{k}-\bar{m}_{i}^{k}|$, we
first derive a simple fact about global solutions of Liouville system.

######  Lemma 3.1.

Let $U=(U_{1},...,U_{n})$ be global solution of

$\Delta
U_{i}+\sum_{j=1}^{n}a_{ij}e^{U_{j}}=0,\quad\mbox{in}\quad\mathbb{R}^{2},\quad\int_{\mathbb{R}^{2}}e^{U_{i}}<\infty,\quad
U_{i}\,\,\,\mbox{is radial}$

and suppose $\max_{i\in I}U_{i}(0)=0$. Then

(3.14) $U_{i}(r)=-m_{i}\log
r+D_{i}-\alpha_{i}-\sum_{j=1}^{n}\frac{a_{ij}}{(m_{j}-2)^{2}}e^{D_{j}-\alpha_{j}}r^{2-m_{j}}+O(r^{2-m-\delta}).$

where
$m_{i}=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}\sum_{j=1}^{n}a_{ij}e^{U_{j}(x)}dx$,
$\alpha_{i}=-U_{i}(0)$ and $D_{i}=\int_{0}^{\infty}\log
r\sum_{j=1}^{n}a_{ij}e^{U_{j}(r)}rdr.$.

Proof of Lemma 3.1:

It is easy to see that

$U_{i}(x)=-\frac{1}{2\pi}\int_{\mathbb{R}^{2}}\log|x-y|\sum_{j=1}^{n}a_{ij}e^{U_{j}}dy+c_{i}.$

Recall that $U_{i}(0)=-\alpha_{i}$, then

$-\alpha_{i}=-\int_{0}^{\infty}\log
r\sum_{j=1}^{n}a_{ij}e^{U_{j}(r)}rdr+c_{i}.$

and

(3.15) $\displaystyle U_{i}(x)$
$\displaystyle=-\frac{1}{2\pi}\int_{\mathbb{R}^{2}}(\log|x-y|-\log|x|)\sum_{j=1}^{n}a_{ij}e^{U_{j}(y)}dy+D_{i}-\alpha_{i}-m_{i}\log|x|$
$\displaystyle=-m_{i}\log|x|+D_{i}-\alpha_{i}+O(|x|^{-\delta}).$

for some $\delta>0$. This expression gives

$e^{U_{i}(r)}=r^{-m_{i}}e^{D_{i}-\alpha_{i}}+O(r^{-m-\delta}).$

Then we use this in the ode that $U_{i}$ satisfies:

$U_{i}^{\prime\prime}(r)+\frac{1}{r}U_{i}^{\prime}(r)=-\sum_{j=1}^{n}a_{ij}e^{U_{j}(r)},\quad
0<r<\infty.$

Here we use the fact that

$\lim_{r\to\infty}rU_{i}^{\prime}(r)=-m_{i}.$

Thus

$\displaystyle-
m_{i}-rU_{i}^{\prime}(r)=-\sum_{j=1}^{n}a_{ij}\int_{r}^{\infty}se^{U_{j}(s)}ds$
$\displaystyle=-\sum_{j=1}^{n}a_{ij}\frac{e^{D_{j}-\alpha_{j}}}{m_{j}-2}r^{2-m_{j}}+O(r^{2-m-\delta}).$

Then we have

$U_{i}^{\prime}(r)=-\frac{m_{i}}{r}+\sum_{j=1}^{n}\frac{a_{ij}}{m_{j}-2}e^{D_{j}-\alpha_{j}}r^{1-m_{j}}+O(r^{1-m-\delta}).$

After integration and using (3.15) we see that (3.14) holds. Lemma 3.1 is
established. $\Box$

The main result in this section is:

######  Proposition 3.1.

If $m<4$,

$|m_{i}^{k}-\bar{m}_{i}^{k}|\leq C\delta_{0}^{4-m}\epsilon_{k}^{m_{k}-2},\quad
i\in I.$

Proof of Proposition 3.1: Let $V_{i}^{k}$ be the sequence of global solutions
that approximate $\tilde{u}_{i}^{k}$ around $p_{k}$, and in a neighborhood
centered at $p_{k}$, $V_{i}^{k}(0)$ agrees with $\tilde{u}_{i}^{k}$ at
$p_{k}$. $\bar{V}_{i}^{k}$ is understood as the first approximation around
$q_{k}$ and we have We use $m_{iv}^{k}$ to denote

$m_{iv}^{k}=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}\sum_{j=1}^{n}a_{ij}e^{V_{j}^{k}}dx,\quad
i=1,..,n,$

and $\bar{m}_{iv}^{k}$ is for $\bar{V}_{i}^{k}$. Note that a rough estimate of
$m_{i}^{k}$ based on Theorem 2.2 gives

$m_{i}^{k}-m_{iv}^{k}=O(\epsilon_{k}^{m_{k}-2}).$

If we use

$U_{i}(y)=V_{i}^{k}(x)+2\log\epsilon_{k},\quad\epsilon_{k}=e^{-\frac{1}{2}M_{k}},\quad
x=\epsilon_{k}y.$

Then the expansion of $V_{i}^{k}$ for $|x|$ is

(3.16) $\displaystyle V_{i}^{k}(x)$
$\displaystyle=-m_{iv}^{k}\log|x|-\frac{m_{iv}^{k}-2}{2}M_{k}+D_{i}-\alpha_{i}$
$\displaystyle-\sum_{j=1}^{n}\frac{a_{ij}}{(m_{jv}^{k}-2)^{2}}e^{D_{j}-\alpha_{j}}\epsilon_{k}^{m_{jv}^{k}-2}|x|^{2-m_{jv}^{k}}+O(\epsilon_{k}^{m-2+\delta}).$

$V^{k}=(V_{i}^{k},....,V_{n}^{k})$ is the sequence of global functions that
serves as the first term in the expansion of $\tilde{u}^{k}$ around $p_{k}$.
Similarly

$\bar{V}_{i}^{k}(x)=-\bar{m}_{iv}^{k}\log|x|-\frac{\bar{m}_{iv}^{k}-2}{2}\bar{M}_{k}+\bar{D}_{i}-\bar{\alpha}_{i}+O(\epsilon_{k}^{m_{k}-2})|x|^{2-\bar{m}_{iv}^{k}}$

Since both $V^{k}=(V_{1}^{k},...,V_{n}^{k})$ and
$\bar{V}^{k}=(\bar{V}_{1}^{k},...,\bar{V}_{n}^{k})$ are radial and satisfy the
same Liouville system. The dependence on initial condition gives

(3.17) $|V_{i}^{k}(x)-(\bar{V}_{i}^{k}(\eta x)+2\log\eta)|\leq
C\sum_{i=1}^{n}|\alpha_{i}-\bar{\alpha}_{i}|\quad\mbox{in}\quad B(0,R).$

where $R>0$ is a constant, $\eta$ is chosen to make the heights equal, in this
case $2\log\eta=M_{k}-\bar{M}_{k}$. We also note that one of
$\alpha_{i}=\bar{\alpha}_{i}$. Here we invoke an important result in [32].
Suppose $\alpha_{1}=0$, the mapping from $(\alpha_{2},...,\alpha_{n})$ to
$(\sigma_{2},...,\sigma_{n})$ is invertible. In fact the following matrix

$\mathbb{M}=\left(\begin{array}[]{ccc}\partial_{\alpha_{2}}\sigma_{2}&...&\partial_{\alpha_{n}}\sigma_{2}\\\
\vdots&\vdots&\vdots\\\
\partial_{\alpha_{2}}\sigma_{n}&...&\partial_{\alpha_{n}}\sigma_{n}\end{array}\right)$

is invertible. Note that in this proposition we use

$\sigma_{i}^{k}=\sum_{j=1}^{n}a^{ij}m_{jv}^{k},\quad\bar{\sigma}_{i}^{k}=\sum_{j=1}^{n}a^{ij}\bar{m}_{jv}^{k}.$

Thus (3.17) can be written as

(3.18) $|V_{i}^{k}(x)-(\bar{V}_{i}^{k}(\eta x)+2\log\eta)|\leq
C\sum_{i=2}^{n}|\sigma_{i}^{k}-\bar{\sigma}_{i}^{k}|.$

By the expansion of $\bar{V}_{i}^{k}$, we find that

(3.19) $\bar{V}_{i}^{k}(\eta
x)+2\log\eta=-\frac{\bar{m}_{iv}^{k}-2}{2}M_{k}-\bar{m}_{iv}^{k}\log|x|+\bar{D}_{i}-\bar{\alpha}_{i}+O(\epsilon_{k}^{m_{k}-2}).$

The difference between $V_{i}^{k}$ and $\bar{V}_{i}^{k}(\eta x)+2\log\eta$
gives

$\displaystyle V_{i}^{k}(x)-(\bar{V}_{i}^{k}(\eta x)+2\log\eta)$
$\displaystyle=$
$\displaystyle(\bar{m}_{iv}^{k}-m_{iv}^{k})\log|x|+\frac{\bar{m}_{iv}^{k}-m_{iv}^{k}}{2}M_{k}+D_{i}-\bar{D}_{i}+\bar{\alpha}_{i}-\alpha_{i}+O(\epsilon_{k}^{m_{k}-2}).$

By the dependence of initial condition and fast decay of $V^{k}$ and
$\bar{V}^{k}$, we have

$|D_{i}-\bar{D}_{i}|\leq C\sum_{i=1}^{n}|\alpha_{i}-\bar{\alpha}_{i}|\leq
C\sum_{i=2}^{n}|\sigma_{i}^{k}-\bar{\sigma}_{i}^{k}|.$

And by the one to one correspondence between $(\alpha_{1},...,\alpha_{n})$ to
$(\sigma_{2},...,\sigma_{n})$ and the smoothness of the mapping we also obtain

$\sum_{i=1}^{n}|\alpha_{i}-\bar{\alpha}_{i}|\leq
C\sum_{i=2}^{n}|\sigma_{i}^{k}-\bar{\sigma}_{i}^{k}|.$

Combining these estimates we have

(3.20) $|\frac{\bar{m}_{iv}^{k}-m_{iv}^{k}}{2}M_{k}|\\\ \leq
C\sum_{j=2}^{n}|\sigma_{j}^{k}-\bar{\sigma}_{j}^{k}|+O(\epsilon_{k}^{m_{k}-2})\quad
i=1,...,n.$

After that we multiply $a^{ij}$ with summation on $i$ and take summation on
$j$, then we have

(3.21)
$\sum_{j=1}^{n}|\bar{\sigma}_{j}^{k}-\sigma_{j}^{k}|\leq\frac{C}{M_{k}}\sum_{l=2}^{n}|\sigma_{l}^{k}-\bar{\sigma}_{l}^{k}|+O(\epsilon_{k}^{m_{k}-2}).$

We thus obtain the following important closeness result:

(3.22)
$\sigma_{i}^{k}-\bar{\sigma}_{i}^{k}=O(\epsilon_{k}^{m_{k}-2})/M_{k},\quad
i=1,...,n.$

Thus we have proved that
$m_{iv}^{k}-\bar{m}_{iv}^{k}=O(\epsilon_{k}^{m_{k}-2}/M_{k})$, using the
expansion of $\tilde{u}_{i}^{k}$ in the calculation of
$\int_{B(p_{k},\delta_{0})}\rho_{i}^{k}h_{i}^{k}e^{u_{i}^{k}}dV_{g}$. It is
easy to see that

$|m_{i}^{k}-m_{iv}^{k}|\leq
C\delta_{0}^{m_{k}-4}\epsilon_{k}^{m_{k}-2},\quad|\bar{m}_{i}^{k}-\bar{m}_{iv}^{k}|\leq
C\delta_{0}^{4-m}\epsilon_{k}^{m_{k}-2}.$

Thus

$|m_{i}^{k}-\bar{m}_{i}^{k}|\leq C\delta_{0}^{4-m}\epsilon_{k}^{m_{k}-2}.$

Proposition 3.1 is established. $\Box$

######  Remark 3.1.

By the same argument for the case $m=4$, we also have

(3.23) $m_{iv}^{k}-\bar{m}_{iv}^{k}=O(\epsilon_{k}^{m_{k}-2}/M_{k}).$

Here $m_{k}\to 4$. Later we shall see $m_{k}$ can be replaced by $4$.
Correspondingly

(3.24)
$m_{i}^{k}-\bar{m}_{i}^{k}=O(\epsilon_{k}^{2}\log\frac{1}{\epsilon_{k}}),$

but the leading term (of the order
$O(\epsilon_{k}^{2}\log\frac{1}{\epsilon_{k}})$ ) can be identified as a local
term that involves curvature at the blowup point.

In the more general situation of $N$ bubbles, if we use $m_{it}^{k}$ to denote
the energy in $B(p_{t}^{k},\delta_{0})$, $m_{itv}^{k}$ to denote the energy of
the global sequence as the first term in the approximation, we have, for
$m<4$:

(3.25) $m_{it}^{k}-m_{is}^{k}=O(\epsilon_{k}^{m_{k}-2}),\quad s\neq t,\quad
m<4$ (3.26)
$m_{itv}^{k}-m_{isv}^{k}=O(\epsilon_{k}^{m_{k}-2}/\log\frac{1}{\epsilon_{k}}),\quad
s\neq t,\quad m<4.$

For $m=4$, (3.26) also holds, and
$m_{it}^{k}-m_{is}^{k}=O(\epsilon_{k}^{m_{k}-2}\log\epsilon_{k}^{-1})$. The
following lemma gives an estimate of $\rho_{i}^{k}-\rho$, which determines as
a consequence that

$O(\epsilon_{k}^{m_{k}-2})=O(\epsilon_{k}^{m-2})\quad\mbox{if}\quad m\in(2,4)$
$O(\epsilon_{k}^{m_{k}-2}\log\frac{1}{\epsilon_{k}})=O(\epsilon_{k}^{2}\log\frac{1}{\epsilon_{k}}),\quad\mbox{if}\quad
m=4.$

######  Lemma 3.2.

If $m<4$

$\rho_{i}^{k}-\rho_{i}=O(\epsilon_{k}^{m-2}),\quad i\in I.$

Proof of Lemma 3.2: Recall that $\rho\in\Gamma_{N}$. Let

$\rho_{it}^{k}=\int_{B(p_{t},\delta)}\rho_{i}h_{i}^{k}e^{u_{i}^{k}}dV_{g},\quad
t=1,...,N,\quad E_{i}^{k}=\rho_{i}^{k}-\sum_{t=1}^{N}\rho_{it}^{k}.$

Here $\delta>0$ is small so that bubbling disks are disjoint. Clearly from
$\int_{M}h_{i}^{k}e^{u_{i}^{k}}=1$ and (3.28) we have

$E_{i}^{k}=O(\epsilon_{k}^{m_{k}-2}),\quad\rho_{i}^{k}=\sum_{i=1}^{N}\rho_{it}^{k}+E_{i}^{k}.$

Let $\sigma_{i}=\frac{\rho_{i}}{2\pi N}$ and
$\sigma_{it}^{k}=\frac{\rho_{it}^{k}}{2\pi}$. Then we write
$\sigma_{i}=\sigma_{it}+s_{it}^{k}$. We have known that $s_{it}^{k}=o(1)$ as
$k\to\infty$. Since $\rho\in\Gamma_{N}$ we know

$\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}\sigma_{i}\sigma_{j}-4\sum_{i=1}^{n}\sigma_{i}=0.$

On the other hand the Pohozaev identity around $p_{t}$ gives

$\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}\sigma_{it}^{k}\sigma_{jt}^{k}-4\sum_{i=1}^{n}\sigma_{it}^{k}=O(\epsilon_{k}^{m_{k}-2}).$

The difference between these two equations gives

(3.27)
$\sum_{i=1}^{n}2(m_{i}-2)s_{it}^{k}+\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}s_{it}^{k}s_{jt}^{k}=O(\epsilon_{k}^{m_{k}-2}),\quad
t=1,...,N.$

Taking the sum of $N$ equations we have

$\sum_{i=1}^{n}2(m_{i}-2)(\sum_{t=1}^{N}s_{it}^{k})+\sum_{t=1}^{N}a_{ij}s_{it}^{k}s_{jt}^{k}=O(\epsilon_{k}^{m_{k}-2}).$

From Proposition 3.1 we know that the difference between any two $s_{it}^{k}$
is $O(\epsilon_{k}^{m_{k}-2})$,
$\sum_{t=1}^{N}s_{it}^{N}=Ns_{i1}^{k}+O(\epsilon_{k}^{m_{k}-2})$. Thus we have

$N\sum_{i=1}^{n}s_{i1}^{k}+\sum_{t=1}^{N}\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}s_{it}^{k}s_{jt}^{k}=O(\epsilon_{k}^{m_{k}-2}).$

By the assumption $\frac{\rho_{i}^{k}-\rho_{i}}{\rho_{j}^{k}-\rho_{j}}\sim 1$
we have $s_{i1}^{k}/s_{j1}^{k}\sim 1$ for all $i,j$ because
$\sum_{t=1}^{N}s_{it}^{k}=\frac{1}{2\pi}(\rho_{i}^{k}-\rho_{i})$. Thus

$s_{it}^{k}=O(\epsilon_{k}^{m_{k}-2}),\quad t=1,...,N.$

Thus Lemma 3.2 is established. $\Box$

Now we are ready to prove Theorem 1.2.

Proof of Theorem 1.2: Due to the closeness between $m_{i}^{k}$ and $m$ (which
is of the order $O(\epsilon_{k}^{m_{k}-2})$ we shall use
$O(\epsilon_{k}^{m-2})$ instead of $O(\epsilon_{k}^{m_{k}-2})$ in Proposition
3.1. By a similar argument, based on (3.24) we clearly have

$\mbox{ For
}m=4,\quad\rho_{i}^{k}-\rho_{i}=O(\epsilon_{k}^{2}\log\frac{1}{\epsilon_{k}}),$

which gives $O(\epsilon_{k}^{m_{k}-2})=O(\epsilon_{k}^{2})$. We thus have

$m_{i,t}^{k}-m_{i,s}^{k}=O(\varepsilon_{k}^{2}\log\frac{1}{\varepsilon_{k}}),\,\,t\not=s,\,\,s,t\in\\{1,\cdots,N\\},\,\,i\in
I.$

Recall that $\sigma_{i,t}^{k}=2\pi\sum_{j\in I}a^{ij}m_{i,t}^{k}$ and
$A=(a_{ij})$ is invertible. This theorem is proved. $\Box$

Because of (3.25) the evaluation of $u_{i}^{k}(x)$ away from bubbling disks
now becomes:

$u_{i}^{k}(x)=\bar{u}_{i}^{k}+2\pi
m_{i}^{k}\sum_{l=1}^{N}G^{*}(x,p_{l}^{k})+O(\epsilon_{k}^{m-2}),\quad x\in
M\setminus(\cup B(p_{t}^{k},\delta_{0})),\quad m<4,$
$u_{i}^{k}(x)=\bar{u}_{i}^{k}+2\pi
m_{i}^{k}\sum_{l=1}^{N}G^{*}(x,p_{l}^{k})+O(\epsilon_{k}^{2}\log\frac{1}{\epsilon_{k}}),\quad
x\in M\setminus(\cup B(p_{t}^{k},\delta_{0})),\quad m=4,$

for $\delta_{0}>0$ small. We recall the notation

$G^{*}(p_{t}^{k},p_{l}^{k})=\left\\{\begin{array}[]{ll}\gamma(p_{t}^{k},p_{t}^{k}),\quad\mbox{if}\quad
t=l,\\\ G(p_{t}^{k},p_{l}^{k}),\quad\mbox{if}\quad t\neq l.\end{array}\right.$

With the information available we are in the position to prove Theorem 1.6.

Proof of Theorem 1.6: In $B(p_{t}^{k},\delta_{0})$ we use
$V_{t}^{k}=(V_{1t}^{k},...,V_{nt}^{k})$ to denote the sequence of global
solutions as the first term in the approximation. In the context of multiple
bubbles, $V_{it}^{k}$ has two expressions:

$V_{it}^{k}=-m_{it}^{k}\log|x|+\bar{u}_{i}^{k}+2\pi
m_{i}^{k}\sum_{l=1}^{N}G^{*}(p_{t}^{k},p_{l}^{k})+\log(\rho_{i}^{k}h_{i}^{k}(p_{t}^{k}))+O(\epsilon_{k}^{m-2}),$

and

$V_{it}^{k}=-m_{itv}^{k}\log|x|-\frac{m_{itv}^{k}-2}{2}M_{k,t}+D_{it}^{k}-\alpha_{it}^{k}+O(\epsilon_{k}^{m-2}),$

where $m_{itv}^{k}=\sum_{j=1}^{n}a_{ij}\sigma_{itv}^{k}$,
$\sigma_{itv}^{k}=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}e^{V_{it}^{k}}$.

By comparing the two expressions of $V_{it}^{k}$ we have

$\displaystyle-m_{it}^{k}\log|x|+\bar{u}_{i}^{k}+2\pi
m_{it}^{k}\sum_{l=1}^{N}G^{*}(p_{t}^{k},p_{l}^{k})+\log(h_{i}^{k}(p_{t}^{k})\rho_{i}^{k})$
$\displaystyle=-m_{itv}^{k}\log|x|-\frac{m_{itv}^{k}-2}{2}M_{k,t}+D_{it}^{k}-\alpha_{it}^{k}+O(\epsilon_{k}^{m-2}).$

where all the notations are clearly understandable under this context. Thus,
for $t\neq s$, based on the two different expression of $\bar{u}_{i}^{k}$, we
have

$\sum_{l=1}^{N}2\pi m_{i}^{k}G^{*}(p_{t}^{k},p_{l}^{k})+\log
h_{i}^{k}(p_{t}^{k})=\sum_{l=1}^{N}2\pi
m_{i}^{k}G^{*}(p_{s}^{k},p_{l}^{k})+\log
h_{i}^{k}(p_{s}^{k})+O(\epsilon_{k}^{m-2}),$

where we have used
$m_{ivt}^{k}-m_{ivs}^{k}=O(\epsilon_{k}^{m-2}/\log\frac{1}{\epsilon_{k}})$ and
used $m_{i}^{k}$ to replace all $m_{it}^{k}$. (1.21) is verified. (1.22) can
be verified similarly. Theorem 1.6 is established. $\Box$

The following expression of $\bar{u}_{i}^{k}$ will be used

(3.28) $\displaystyle\bar{u}_{i}^{k}$
$\displaystyle=(1-\frac{m_{i}^{k}}{2})M_{k,t}-\log(\rho_{i}^{k}h_{i}^{k}(p_{t}^{k}))-2\pi
m_{i}^{k}\sum_{l=1}^{N}G^{*}(p_{t}^{k},p_{l}^{k})$
$\displaystyle+D_{i}^{k}-\alpha_{i}^{k}+O(\epsilon_{k}^{m_{k}-2}),\quad
t=1,...,N,$

where we used $\alpha_{i}^{k}$ and $D_{i}^{k}$ to replace any $D_{it}^{k}$ and
$\alpha_{it}^{k}$ for obvious reasons. One of $M_{k,t}$ is $M_{k}$ and the
difference between any two of these is bounded. In fact, for $t\neq s$, the
difference on equations (3.28) for $t$ and $s$ gives

$2\pi
m_{i}(\sum_{l=1}^{N}(G^{*}(p_{t}^{k},p_{l}^{k})-G^{*}(p_{s}^{k},p_{l}^{k}))+\log\frac{h_{i}^{k}(p_{t}^{k})}{h_{i}^{k}(p_{s}^{k})}=-\frac{m_{i}-2}{2}(M_{k,t}-M_{k,s})+O(\epsilon_{k}^{m-2}).$

Equivalent form is

(3.29) $exp(2\pi
m_{i}\sum_{l=1}^{N}(G^{*}(p_{t}^{k},p_{l}^{k})-G^{*}(p_{s}^{k},p_{l}^{k})))\,\frac{h_{i}^{k}(p_{t}^{k})}{h_{i}^{k}(p_{s}^{k})}=\frac{\epsilon_{k,t}^{m_{i}-2}}{\epsilon_{k,s}^{m_{i}-2}}+O(\epsilon_{k}^{m-2}),$

where $\epsilon_{k,t}=e^{-\frac{1}{2}M_{k,t}}$. Also, (3.28) gives

(3.30)
$e^{\bar{u}_{i}^{k}}=\epsilon_{k,t}^{m_{i}-2}\frac{e^{D_{i}^{k}-\alpha_{i}^{k}}}{\rho_{i}^{k}h_{i}^{k}(p_{t}^{k})}e^{-2\pi
m_{i}\sum_{l=1}^{N}G^{*}(p_{t}^{k},p_{l}^{k})}+O(\epsilon_{k}^{m-2+\delta}).$

##  4\. Proof of leading terms

Now we are in the position to prove Theorem 1.3.

Proof of Theorem 1.3:

Since $\int_{M}h_{i}^{k}e^{u_{i}^{k}}dV_{g}=1$ we write

$\rho_{i}^{k}=\sum_{t=1}^{N}\int_{B(p_{t}^{k},\delta_{0})}\rho_{i}^{k}h_{i}^{k}e^{u_{i}^{k}}dV_{g}+\int_{M\setminus\cup_{t}B(p_{t}^{k},\delta_{0})}\rho_{i}^{k}h_{i}^{k}e^{u_{i}^{k}}dV_{g}=\sum_{t=1}^{N}\rho_{it}^{k}+\rho_{ib}^{k}$

where in local coordinates

$\rho_{it}^{k}=\int_{B(0,\delta_{0})}\tilde{h}_{i}^{k}e^{\tilde{u}_{i}^{k}}d\eta.$

Let

$I_{1}=\\{i\in I;\quad\lim_{k\to\infty}m_{i}^{k}=m.\quad\\}.$

Based on Proposition 3.1 for each $i\in I_{1}$,
$m_{i}^{k}-m=O(\epsilon_{k}^{m-2})$. If we use $V_{i}^{k}$ to be the leading
term in the approximation of $\tilde{u}_{i}^{k}$ and $U_{i}^{k}$ be the scaled
version of $V_{i}^{k}$, by (3.16) we have (since $m<4$)

(4.1) $\displaystyle\frac{1}{2\pi}\rho_{it}^{k}$
$\displaystyle=\frac{1}{2\pi}\int_{B(0,\delta_{0}\epsilon_{k}^{-1})}\tilde{h}_{i}^{k}(0)e^{V_{i}^{k}(y)}dy+o(\delta_{0})\epsilon_{k}^{m-2},$
$\displaystyle=\sigma_{iv}^{k}-\frac{e^{D_{i}-\alpha_{i}}}{m-2}\epsilon_{k,t}^{m-2}\delta_{0}^{2-m}+E_{\delta_{0}},\quad
i\in I_{1}.$

where we use
$\sigma_{iv}^{k}=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}e^{V_{it}^{k}}$ and
$E_{\delta_{0}}$ to denote $o(\delta_{0})\epsilon_{k}^{m-2}$,
$\epsilon_{k,t}=e^{-\frac{M_{k,t}}{2}}$. We did not use $t$ in
$\sigma_{iv}^{k}$ because the difference between any two of them is
$O(\epsilon_{k}^{m-2}/M_{k})$. Now for $i\not\in I_{1}$ we have

(4.2) $\frac{1}{2\pi}\rho_{it}^{k}=\sigma_{iv}^{k}+E_{\delta_{0}},\quad
i\not\in I_{1}$

and

(4.3) $|\rho_{ib}^{k}|=E_{\delta_{0}},\quad i\not\in I_{1}.$

Combining (4.1), (4.2), (4.3) we have

$\displaystyle\sum_{i=1}^{n}\frac{4}{2\pi}\rho_{it}^{k}-\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}\frac{\rho_{it}^{k}}{2\pi}\frac{\rho_{jt}^{k}}{2\pi}$
$\displaystyle=$ $\displaystyle
4\sum_{i=1}^{n}(\sigma_{iv}-\frac{e^{D_{i}-\alpha_{i}}}{m_{i}-2}\delta^{2-m_{i}}\epsilon_{k,t}^{m_{i}-2})$
$\displaystyle-\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}(\sigma_{iv}^{k}-\frac{e^{D_{i}-\alpha_{i}}}{m_{i}-2}\delta^{2-m_{i}}\epsilon_{k,t}^{m_{i}-2})(\sigma_{jv}^{k}-\frac{e^{D_{j}-\alpha_{j}}}{m_{j}-2}\delta^{2-m_{j}}\epsilon_{k,t}^{m_{j}-2})+E_{\delta_{0}}$
$\displaystyle=$ $\displaystyle-\frac{4}{m-2}\sum_{i\in
I_{1}}e^{D_{i}-\alpha_{i}}\delta_{0}^{2-m}\epsilon_{k,t}^{m-2}+2\sum_{i\in
I_{1}}\sum_{j=1}^{n}a_{ij}\sigma_{jv}^{k}\frac{e^{D_{j}-\alpha_{j}}}{m_{j}-2}\delta_{0}^{2-m_{j}}\epsilon_{k,t}^{m_{j}-2}+E_{\delta_{0}}$
$\displaystyle=$ $\displaystyle
2\delta_{0}^{2-m}\epsilon_{k,t}^{m-2}\sum_{i\in
I_{1}}e^{D_{i}-\alpha_{i}}+E_{\delta_{0}}.$

Note that
$m_{i}=\sum_{j=1}^{n}a_{ij}\sigma_{jv}^{k}+O(\epsilon_{k}^{m-2}/M_{k})$. For
$i\in I_{1}$, using (3.28) we have

$\displaystyle\rho_{ib}^{k}$
$\displaystyle=\int_{M\setminus(\cup_{t}B(p_{t}^{k},\delta_{0}))}\rho_{i}^{k}h_{i}^{k}e^{u_{i}^{k}}dV_{g}$
$\displaystyle=$
$\displaystyle\rho_{i}^{k}e^{\bar{u}_{i}^{k}}\int_{M\setminus\cup_{t}B(p_{t}^{k},\delta_{0})}h_{i}^{k}e^{2\pi
m\sum_{t=1}^{N}G(x,p_{t}^{k})}+E_{\delta_{0}},\quad i\in I_{1}.$

Now we define $N$ open sets $\Omega_{t,\delta_{0}}$ such that they are
mutually disjoint, each of them contains a bubbling disk and their union is
$M$:

$B(p_{t}^{k},\delta_{0})\subset\Omega_{t,\delta_{0}},\quad\cup_{t=1}^{N}\overline{\Omega_{t,\delta_{0}}}=M,\quad\Omega_{t,\delta_{0}}\cap\Omega_{s,\delta_{0}}=\emptyset,\,\,\forall
t\neq s.$

In each $\Omega_{t,\delta_{0}}$ we use (3.30)) to write $\rho_{ib}^{k}$ as
(for $i\in I_{1}$)

(4.4)
$\displaystyle\rho_{ib}^{k}=\rho_{i}^{k}e^{\bar{u}_{i}^{k}}\sum_{t=1}^{N}\int_{\hat{\Omega}_{t,\delta_{0}}}h_{i}^{k}e^{2\pi
m\sum_{l=1}^{N}G(x,p_{l}^{k})}$
$\displaystyle=\sum_{t=1}^{N}\int_{\hat{\Omega}_{t,\delta_{0}}}\epsilon_{k,t}^{m-2}\frac{h_{i}^{k}(x)}{h_{i}^{k}(p_{t}^{k})}e^{D_{i}-\alpha_{i}}e^{2\pi
m\sum_{l=1}^{N}(G(x,p_{l}^{k})-G^{*}(p_{t}^{k},p_{l}^{k}))}dV_{g}$
$\displaystyle+E_{\delta_{0}},\quad i\in I_{1},$

where $\hat{\Omega}_{t,\delta_{0}}=\Omega_{t,\delta_{0}}\setminus
B(p_{t}^{k},\delta_{0})$. Now we put estimates together to have

$\displaystyle 4\sum_{i=1}^{n}\frac{\rho_{i}^{k}}{2\pi
N}-\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}\frac{\rho_{i}^{k}}{2\pi
N}\frac{\rho_{j}^{k}}{2\pi N}$ $\displaystyle=$ $\displaystyle
4\sum_{i=1}^{n}(\sum_{t=1}^{N}\frac{\rho_{it}^{k}}{2\pi
N}+\frac{\rho_{ib}}{2\pi
N})-\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}(\sum_{t=1}^{N}\frac{\rho_{it}^{k}}{2\pi
N}+\frac{\rho_{ib}^{k}}{2\pi N})(\sum_{s=1}^{N}\frac{\rho_{js}^{k}}{2\pi
N}+\frac{\rho_{jb}^{k}}{2\pi N})$ $\displaystyle=$
$\displaystyle\frac{2}{N}\sum_{t=1}^{N}\sum_{i\in
I_{1}}\delta_{0}^{2-m}\epsilon_{k,t}^{m-2}e^{D_{i}-\alpha_{i}}-2(m-2)\sum_{i\in
I_{1}}\frac{\rho_{ib}^{k}}{2\pi N}+E_{\delta_{0}}.$

Using (4.4) in the expression above we have

(4.5) $\displaystyle\qquad 4\sum_{i=1}^{n}\frac{\rho_{i}^{k}}{2\pi
N}-\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}\frac{\rho_{i}^{k}}{2\pi
N}\frac{\rho_{j}^{k}}{2\pi N}$ $\displaystyle=\frac{2}{N}\sum_{i\in
I_{1}}e^{D_{i}-\alpha_{i}}\sum_{t=1}^{N}\epsilon_{k,t}^{m-2}\bigg{(}\delta_{0}^{2-m}-\frac{(m-2)}{2\pi}\int_{\hat{\Omega}_{t,\delta_{0}}}(\frac{h_{i}^{k}(x)}{h_{i}^{k}(p_{t}^{k})}e^{2\pi
m\sum_{l=1}^{N}(G(x,p_{l}^{k})-G^{*}(p_{t}^{k},p_{l}^{k}))}dV_{g}\bigg{)}$
$\displaystyle+E_{\delta_{0}}$

Since one of $\epsilon_{k,t}$ is $\epsilon_{k}$, say
$\epsilon_{1,t}=\epsilon_{k}$, based on (3.29) we use

(4.6)
$c_{t}=\epsilon_{k,t}^{m-2}/\epsilon_{k}^{m-2}=\frac{h_{i}^{k}(p_{t}^{k})e^{2\pi
m\sum_{l=1}^{N}G^{*}(p_{t}^{k},p_{l}^{k})}}{h_{i}^{k}(p_{1}^{k})e^{2\pi
m\sum_{l=1}^{N}G^{*}(p_{1}^{k},p_{l}^{k})}},$

Here we note that (3.29) implies that for $i\in I_{1}$, $c_{t}$ is independent
of $i\in I_{1}$. Then we use $D+o(1)$ to represent

(4.7) $\displaystyle D+o(1):=$ $\displaystyle\sum_{i\in
I_{1}}e^{D_{i}-\alpha_{i}}\sum_{t=1}^{N}c_{t}\bigg{(}\delta_{0}^{2-m}-\frac{(m-2)}{2\pi}\int_{\hat{\Omega}_{t,\delta_{0}}}(\frac{h_{i}^{k}(x)}{h_{i}^{k}(p_{t}^{k})}e^{2\pi
m\sum_{l=1}^{N}(G(x,p_{l}^{k})-G^{*}(p_{t}^{k},p_{l}^{k}))}dV_{g}\bigg{)}.$

and the leading term of $4\sum_{i=1}^{n}\frac{\rho_{i}^{k}}{2\pi
N}-\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{\rho_{i}^{k}}{2\pi
N}\frac{\rho_{j}^{k}}{2\pi N}$ is written as

$4\sum_{i=1}^{n}\frac{\rho_{i}^{k}}{2\pi
N}-\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{\rho_{i}^{k}}{2\pi
N}\frac{\rho_{j}^{k}}{2\pi N}=\frac{2}{N}\epsilon_{k}^{m-2}(D+o(1)),$

where $o(1)$ stands for the infinitesimal quantity when $\delta_{0}\to 0$.
Here it is important to observe that $D$ is involved with integration on the
whole manifold. Theorem 1.3 is established. $\Box$

Proof of Theorem 1.5:

Around each $p_{t}^{k}$, an extension of (3.11) can be easily determined to be

(4.8) $\displaystyle\phi_{i}^{k}(x)$ $\displaystyle=2\pi
m_{i}(\gamma(x,p_{t}^{k})-\gamma(p_{t}^{k},p_{t}^{k}))$
$\displaystyle+\sum_{l\neq t}2\pi
m_{i}(G^{*}(x,p_{l}^{k})-G^{*}(p_{t}^{k},p_{l}^{k}))-f_{i}^{k}(x)+E_{1}$

where $E_{1}=O(\epsilon_{k}^{m-2})$ if $m<4$ and is
$O(\epsilon_{k}^{2}\log\frac{1}{\epsilon_{k}})$ if $m=4$. Correspondingly,

(4.9) $\nabla\phi_{i}^{k}(p_{t}^{k})=2\pi
m_{i}\sum_{l=1}^{N}\nabla_{1}G^{*}(p_{t}^{k},p_{l}^{k})+E_{1}.$

With these notations, (1.18) and (1.19) follow immediately. Theorem 1.5 is
established. $\Box$

Finally we prove Theorem 1.4.

Proof of Theorem 1.4:

$\rho_{i}^{k}=\sum_{t=1}^{N}\int_{B(p_{t}^{k},\delta_{0})}\rho_{i}^{k}h_{i}^{k}e^{u_{i}^{k}}dV_{g}+\int_{M\setminus\cup_{t}B(p_{t}^{k},\delta_{0})}\rho_{i}^{k}h_{i}^{k}e^{u_{i}^{k}}dV_{g}.$

We continue to use the notation $\rho_{it}^{k}$ and $\rho_{ib}^{k}$. By (3.2)
and Theorem 2.3 the second integral is $O(\epsilon_{k}^{2})$, this is the same
as the computation for the single equation [12].

Now we use the expansion of bubbles to compute each $\rho_{it}^{k}$. By the
expansion of $\tilde{u}_{i}^{k}$ around $p_{t}^{k}$, we have

$\displaystyle\rho_{it}^{k}=\int_{B(p_{t}^{k},\delta_{0})}\rho_{i}^{k}h_{i}^{k}e^{u_{i}^{k}}dV_{g}=\int_{B(0,\delta_{0})}\tilde{h}_{i}^{k}e^{\phi_{i}^{k}}e^{\tilde{u}_{i}^{k}-\phi_{i}^{k}}d\eta$
(4.10)
$\displaystyle=\int_{B(0,\delta_{0}\epsilon_{k,t}^{-1})}\rho_{i}^{k}h_{i}^{k}(p_{t}^{k})e^{U_{i}^{k}(\eta)}d\eta+O(\epsilon_{k}^{2})$
$\displaystyle+\int_{B(0,\delta_{0}\epsilon_{k,t}^{-1})}\epsilon_{k,t}^{2}(\frac{1}{4}\Delta\tilde{h}_{i}^{k}(0)+\frac{1}{2}\nabla\tilde{h}_{i}^{k}(0)\cdot\nabla\phi_{i}^{k}(0)+\frac{1}{4}|\nabla\phi_{i}^{k}(0)|^{2})|y|^{2}e^{U_{i}^{k}}d\eta.$

The first integral on the right hand side of the above is
$O(\epsilon_{k}^{2})$ different from the global solution in the approximation
of $\tilde{u}_{i}^{k}$ around $p_{t}^{k}$. So we use $\sigma_{ivt}^{k}$ to
denote it. For $t\neq s$, from (3.23) we see that

$\sigma_{ivt}^{k}-\sigma_{ivs}^{k}=O(\epsilon_{k}^{2}/\log\frac{1}{\epsilon_{k}}).$

To evaluate the last term, we first use the definition of the
$\tilde{h}_{i}^{k}$ to have

$\nabla\tilde{h}_{i}^{k}(0)\cdot\nabla\phi_{i}^{k}(0)=2\pi m_{i}\frac{\nabla
h_{i}^{k}(p_{t}^{k})}{h_{i}^{k}(p_{t}^{k})}\cdot\sum_{l=1}^{N}\nabla_{1}G^{*}(p_{t}^{k},p_{l}^{k})+O(\epsilon_{k}^{2}\log\frac{1}{\epsilon_{k}}),$

and

$\displaystyle\Delta\tilde{h}_{i}^{k}(0)=\frac{\Delta
h_{i}^{k}(p_{t}^{k})}{h_{i}^{k}(p_{t}^{k})}-2K(p_{t}^{k})+\sum_{j=1}^{n}a_{ij}\rho_{j}+O(\epsilon_{k}^{2}\log\frac{1}{\epsilon_{k}})$
(4.11) $\displaystyle=\frac{\Delta
h_{i}^{k}(p_{t}^{k})}{h_{i}^{k}(p_{t}^{k})}-2K(p_{t}^{k})+8\pi
N+O(\epsilon_{k}^{2}\log\frac{1}{\epsilon_{k}}).$

Then we define $b_{it}^{k}$ as

(4.12) $\displaystyle b_{it}^{k}$
$\displaystyle=e^{D_{i}-\alpha_{i}}\bigg{(}\frac{1}{4}\frac{\Delta
h_{i}^{k}(p_{t}^{k})}{h_{i}^{k}(p_{t}^{k})}-K(p_{t}^{k})+4\pi N$
$\displaystyle+4\pi\frac{\nabla
h_{i}^{k}(p_{t}^{k})}{h_{i}^{k}(p_{t}^{k})}\cdot\sum_{l=1}^{N}\nabla_{1}G^{*}(p_{t}^{k},p_{l}^{k})+16\pi^{2}|\sum_{l=1}^{N}\nabla_{1}G^{*}(p_{t}^{k},p_{l}^{k})|^{2}\bigg{)}.$

With this $b_{it}^{k}$ we have

$\frac{\rho_{it}^{k}}{2\pi}=\sigma_{iv}^{k}+b_{it}^{k}\epsilon_{k,t}^{2}\log\frac{1}{\epsilon_{k}}+O(\epsilon_{k}^{2}).$

From (3.29) we see that $\epsilon_{k,t}$ can be replaced by $\epsilon_{k}$.
Consequently,

(4.13) $\displaystyle 4\sum_{i=1}^{n}\frac{\rho_{i}^{k}}{2\pi
N}-\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}\frac{\rho_{i}^{k}}{2\pi
N}\frac{\rho_{j}^{k}}{2\pi N}$ $\displaystyle=$ $\displaystyle
4\sum_{i=1}^{n}\sum_{t=1}^{N}\frac{\rho_{it}^{k}}{2N\pi}-\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}(\sum_{t=1}^{N}\frac{\rho_{it}}{2\pi
N})(\sum_{s=1}^{N}\frac{\rho_{js}}{2\pi N})$ $\displaystyle=$ $\displaystyle
4\sum_{i=1}^{n}\sum_{t=1}^{N}(\frac{\sigma_{ivt}}{N}+\tilde{\epsilon}_{k}\frac{b_{it}^{k}}{N})-\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{t=1}^{N}\sum_{s=1}^{N}a_{ij}(\frac{\sigma_{ivt}}{N}+\tilde{\epsilon}_{k}\frac{b^{k}_{it}}{N})(\frac{\sigma_{jv}}{N}+\tilde{\epsilon}_{k}\frac{b^{k}_{js}}{N})$
$\displaystyle=$ $\displaystyle
4\sum_{i=1}^{n}\frac{\sigma_{iv}}{N}-\sum_{i=1}^{n}\sum_{j=1}^{n}a_{ij}\frac{\sigma_{iv}\sigma_{jv}}{N}+4\sum_{i=1}^{n}\sum_{t=1}^{N}\tilde{\epsilon}_{k}\frac{b^{k}_{it}}{N}-2\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{t=1}^{N}a_{ij}\sigma_{iv}\tilde{\epsilon}_{k}\frac{b^{k}_{jt}}{N}+O(\epsilon_{k}^{2})$
$\displaystyle=$
$\displaystyle-4\epsilon_{k}^{2}\log\epsilon_{k}^{-1}\sum_{i=1}^{n}\sum_{t=1}^{N}b_{it}^{k}+O(\epsilon_{k}^{2}).$

where $\tilde{\epsilon}_{k}$ stands for
$\epsilon_{k}^{2}\log\frac{1}{\epsilon_{k}}$. Theorem 1.4 is established.
$\Box$

## Acknowledgment

The authors are grateful to Professor C.-S. Lin for his helpful discussions
and comments on the article.

## References

  * [1] J. J. Aly, Thermodynamics of a two-dimensional self-gravitating system, Phy. Rev. A 49, No. 5, Part A (1994), 3771–3783.
  * [2] D. Bartolucci, Compactness result for periodic multivortices in the electroweak theory, Nonlinear Anal., 53 (2003), 277–297.
  * [3] D. Bartolucci and G. Tarantello, Liouville type equations with singular data and their applications to periodic multivortices for the electroweak theory, Comm. Math. Phys., 229 (2002), 3–47.
  * [4] W. H. Bennet, Magnetically self-focusing streams, Phys. Rev. 45 (1934), 890–897.
  * [5] L. A. Caffarelli, Y. Yang, Vortex condensation in the Chern-Simons Higgs model: an existence theorem, Comm. Math. Phys. 168 (1995), no. 2, 321–336.
  * [6] E. Caglioti, P. L. Lions, C. Marchioro, C, M. Pulvirenti, A special class of stationary flows for two-dimensional Euler equations: a statistical mechanics description, Comm. Math. Phys. 143 (1992), no. 3, 501–525.
  * [7] E. Caglioti, P. L. Lions, C. Marchioro, C, M. Pulvirenti, A special class of stationary flows for two-dimensional Euler equations: a statistical mechanics description. II, Comm. Math. Phys. 174 (1995), no. 2, 229–260.
  * [8] S. A. Chang, C. C. Chen, C. S. Lin, Extremal functions for a mean field equation in two dimension. (English summary) Lectures on partial differential equations, 61–93, New Stud. Adv. Math., 2, Int. Press, Somerville, MA, 2003.
  * [9] S. Chanillo, M. K-H Kiessling, Rotational symmetry of solutions of some nonlinear problems in statistical mechanics and in geometry, Comm. Math. Phys. 160 (1994), no. 2, 217–238.
  * [10] S. Chanillo, M. K-H Kiessling, Conformally invariant systems of nonlinear PDE of Liouville type, Geom. Funct. Anal. 5 (1995), no. 6, 924–947.
  * [11] J. L. Chern, Z. Y. Chen and C. S. Lin, Uniqueness of topological solutions and the structure of solutions for the Chern-Simons with two Higgs particles, Comm. Math. Phys. 296 (2010), 323–351.
  * [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] C. C. Chen, C. S. Lin, Topological degree for a mean field equation on Riemann surfaces. Comm. Pure Appl. Math. 56 (2003), no. 12, 1667–1727.
  * [14] S. Childress and J. K. Percus, Nonlinear aspects of Chemotaxis, Math. Biosci. 56 (1981), 217–237.
  * [15] M. Chipot, I. Shafrir, G. Wolansky, On the solutions of Liouville systems. J. Differential Equations 140 (1997), no. 1, 59–105.
  * [16] M. Chipot, I. Shafrir, G. Wolansky, Erratum: “On the solutions of Liouville systems” [J. Differential Equations 140 (1997), no. 1, 59–105; MR1473855 (98j:35053)]. J. Differential Equations 178 (2002), no. 2, 630.
  * [17] P. Debye and E. Huckel, Zur Theorie der Electrolyte, Phys. Zft 24 (1923), 305–325.
  * [18] G. Dunne, Self-dual Chern-Simons Theories, Lecture Notes in Physics, vol. m36, Berlin: Springer-Verlag, 1995.
  * [19] J. Dziarmaga, Low energy dynamics of $[U(1)]^{N}$ Chern-Simons solitons and two dimensional nonlinear equations, Phys. Rev. D 49 (1994), 5469–5479.
  * [20] L. Ferretti, S.B. Gudnason, K. Konishi, Non-Abelian vortices and monopoles in SO(N) theories, Nuclear Physics B, Volume 789, Issues 1-2, 21 January 2008, Pages 84-110.
  * [21] Y. Gu, L. Zhang, Degree counting theorems for singular Liouville systems, with Yi Gu. to appear on Annali della Scuola Normale Superiore di Pisa, Classe di Scienze, DOI Number: $10.2422/2036-2145.201812_{0}07$
  * [22] H. Huang Existence of bubbling solutions for the Liouville system in a torus. Calc. Var. Partial Differential Equations 58 (2019), no. 3, Paper No. 99, 26 pp.
  * [23] H. Huang ; L. Zhang, The domain geometry and the bubbling phenomenon of rank two gauge theory. Comm. Math. Phys. 349 (2017), no. 1, 393–424.
  * [24] H. Huang ; L. Zhang, Local uniqueness of bubbling solutions to the Liouville system, in preparation
  * [25] J. Hong, Y. Kim and P. Y. Pac, Multivortex solutions of the Abelian Chern-Simons-Higgs theory, Phys. Rev. Letter 64 (1990), 2230–2233.
  * [26] R. Jackiw and E. J. Weinberg, Selfdual Chern Simons vortices, Phys. Rev. Lett. 64 (1990), 2234–2237.
  * [27] E. F. Keller and L. A. Segel, Traveling bands of Chemotactic Bacteria: A theoretical analysis, J. Theor. Biol. 30 (1971), 235–248.
  * [28] M. K.-H. Kiessling, Statistical mechanics of classical particles with logarithmic interactions, Comm. Pure Appl. Math. 46 (1993), no. 1, 27–56.
  * [29] M. K.-H. Kiessling, Symmetry results for finite-temperature, relativistic Thomas-Fermi equations, Commun. Math. Phys. 226 (2002), no. 3, 607–626.
  * [30] M. K.-H. Kiessling and J. L. Lebowitz, Dissipative stationary Plasmas: Kinetic Modeling Bennet Pinch, and generalizations, Phys. Plasmas 1 (1994), 1841–1849.
  * [31] C. Kim, C. Lee and B.-H. Lee, Schrödinger fields on the plane with $[U(1)]^{N}$ Chern-Simons interactions and generalized self-dual solitons, Phys. Rev. D (3) 48 (1993), 1821–1840.
  * [32] C. S. Lin, L. Zhang, Profile of bubbling solutions to a Liouville system. Ann. Inst. H. Poincare Anal. Non Lineaire 27 (2010), no. 1, 117–143.
  * [33] C. S. Lin, L. Zhang, A topological degree counting for some Liouville systems of mean field equations, Comm. Pure Appl. Math. volume 64, Issue 4, pages 556–590, April 2011.
  * [34] C. S. Lin, L. Zhang, On Liouville systems at critical parameters, Part 1: One bubble. J. Funct. Anal. 264 (2013), no. 11, 2584–2636.
  * [35] M. Nolasco, G. Tarantello, Vortex condensates for the SU(3) Chern-Simons theory, Comm. Math. Phys. 213 (2000), no. 3, 599–639.
  * [36] I. Rubinstein, Electro diffusion of Ions, SIAM, Stud. Appl. Math. 11 (1990).
  * [37] J. Spruck, Y. Yang, Topological solutions in the self-dual Chern-Simons theory: existence and approximation, Ann. Inst. H. Poincare Anal. Non Lineaire 12 (1995), no. 1, 75–97.
  * [38] F. Wilczek, Disassembling anyons. Physical review letters 69.1 (1992): 132.
  * [39] G. Wolansky, On steady distributions of self-attracting clusters under friction and fluctuations, Arch. Rational Mech. Anal. 119 (1992), 355–391.
  * [40] G. Wolansky, On the evolution of self-interacting clusters and applications to semi-linear equations with exponential nonlinearity, J. Anal. Math. 59 (1992), 251–272.
  * [41] G. Wolansky. Multi-components chemotactic system in the absence of conflicts. European Journal of Applied Mathematics, Volume 13, Issue 6, 2002.
  * [42] Y. Yang, Solitons in field theory and nonlinear analysis, Springer-Verlag, 2001.
  * [43] L. Zhang, Blowup solutions of some nonlinear elliptic equations involving exponential nonlinearities. Comm. Math. Phys. 268 (2006), no. 1, 105–133.
  * [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.


2021

These authors contributed equally to this work.

These authors contributed equally to this work.

[1]J. Y. Süngü These authors contributed equally to this work.

These authors contributed equally to this work.

[1]Physics, Kocaeli University, Umuttepe Yerleşkesi, Izmit, 41001, Turkey

# Possible Molecular Explanation for the Resonance $Y(4500)$

E. Güngör<EMAIL_ADDRESS>H. Sundu<EMAIL_ADDRESS><EMAIL_ADDRESS>E. V.Veliev<EMAIL_ADDRESS>*

###### Abstract

The BESIII collaboration has discovered a new state with hidden charm-strange.
Its mass is intriguingly close to the $D_{s}\bar{D}_{s1}$ threshold and does
not have the properties of the charmonium states. Working with the QCD sum
rules (QCDSR) approach, we test if the charmonium-like structure $Y(4500)$,
detected in the invariant mass spectrum $K^{+}K^{-}J/\psi$ may be interpreted
as an exotic $D_{s}\bar{D}_{s1}$ molecular structure with $J^{PC}=1^{--}$.
Considering the contributions of QCD condensates up to operator dimension ten,
we estimate the mass and decay constant of $Y(4500)$ resonance. We get
$m_{Y}=(4488.35\pm 11.54)$ MeV in excellent agreement with the meson mass
reported by BESIII and $f_{Y}=(4.04\pm 0.36)\times 10^{-3}$$\mathrm{GeV^{4}}$.
Our findings indicate that a pseudoscalar-axialvector molecule current can
well describe this state.

###### keywords:

QCD Sum rules, Exotic molecule meson, Hidden charm-strange mesons, Mass and
decay constant of mesons

##  1 Introduction

Observing a plethora of the surprising charmonium-like states named $XYZ$
states, which do not fit the quarkonia interpretation at the $B$ factories and
the Tevatron, has been a new improvement in high energy physics and
highlighted the need for a more comprehensive theoretical understanding of the
hadron spectrum. All $XYZ$ states are above or at least in the vicinity of
open-charm thresholds. So, we expect them to decay into open-charm states
because of the OZI rule dominantly. But, many have only decayed to charmonium
and light mesons/photons. Considering these unusual properties, the $Y$-states
are regarded as promising candidates for unconventional hadron states. Their
production mechanism, spin-parity assignments, masses, decay constants,
widths, and modes have been widely discussed in review papers Chen:2022asf ;
Brambilla:2019esw ; Liu:2013waa ; Chen:2016qju ; Agaev:2020zad .

Over the past several decades, many charmonium-like exotic candidates have
been announced, such as the scalar fully open-flavor structure $X_{0}(2900)$
Albuquerque:2020ugi ; Agaev:2022eeh , hidden charm-strange scalar meson
$X(3960)$ Agaev:2022pis ; Mutuk:2022ckn , scalar $X(5568)$ with four different
quark flavors Agaev:2016urs ; Sungu:2019ybf , axial-vector hidden charm-
strange state $Z_{cs}(3985)$ Maiani:2021tri ; Wang:2020rcx ; Sungu:2020zvk ,
axial-vector hidden charm $X(3872)$ Mutuk:2018zxs ; Sundu:2016oda , axial-
vector charmed-strange state $X(4140)$ Agaev:2017foq ; Turkan:2017pil , vector
charmed-strange $X(4630)$ Agaev:2022iha , vector hidden bottom-strange
$Y_{b}(10890)$ Ali:2013xba ; Sungu:2018iew , pseudo-vector meson $Z_{c}(3900)$
and its excited state $Z(4430)$ Azizi:2020yhs ; Agaev:2017tzv , etc.

Possible interpretations include meson molecules, tetraquarks, glueballs,
hybrid mesons and pentaquarks; but the status of these exotic particles is not
yet fully clarified. Thus, searches of all production and decay patterns plus
the measurement of resonance parameters from other experiments, and also
calculations from various theoretical studies are worthwhile in giving
information that will help reveal the nature of $XYZ$ states Chen:2022asf ;
Brambilla:2019esw ; Liu:2013waa ; Chen:2016qju ; Agaev:2020zad ; Wang:2021gml
; Tiwari:2022azj ; Liu:2021xje .

Specifically interesting among them are the vector mesons detected in
$e^{+}e^{-}$ annihilation. Recently, a series of charmonium-like resonances
with $J^{PC}=1^{--}$, so-called $Y$ states, have been reported by countless
experiments. Still, the inner structure of the $Y$ states remains uncertain.
For instance, the state $Y(4230)$, earlier named $Y(4260)$, is the first
detected vector charmonium-like state in the invariant-mass spectrum of
$\pi^{+}\pi^{-}J/\psi$ in 2005 by BaBar detector BaBar:2005hhc . Also, it
decays into $\pi^{0}\pi^{0}J/\psi$, and $K^{+}K^{-}J/\psi$ were reported in a
study of $12.6pb^{-1}$ data collected at $4.26$ GeV by the CLEO-c experiment
CLEO:2006ike  and confirmed by Belle Belle:2007dxy .

Recently, $Y(4230)$ and $Y(4500)$ have been observed by BESIII Collaboration
in $e^{+}e^{-}\rightarrow K^{+}K^{-}J/\psi$ annihilation BESIII:2022joj . The
second new structure is detected for the first time with a statistical
significance greater than $8\sigma$, called $Y(4500)$. The mass and width of
the newly observed resonance are measured to be:

$\displaystyle M_{Y(4500)}$ $\displaystyle=$ $\displaystyle 4484.7\pm 13.3\pm
24.1~{}\mathrm{MeV},$ $\displaystyle\Gamma_{Y(4500)}$ $\displaystyle=$
$\displaystyle 111.1\pm 30.1\pm 15.2~{}\mathrm{MeV}.$ (1)

and also in the $e^{+}e^{-}\rightarrow D^{\ast 0}D^{\ast-}\pi^{+}+c.c$ process
the following results are obtained BESIII:2023unv :

$\displaystyle M_{Y(4500)}$ $\displaystyle=$ $\displaystyle 4469.1\pm 26.2\pm
3.6~{}\mathrm{MeV},$ $\displaystyle\Gamma_{Y(4500)}$ $\displaystyle=$
$\displaystyle 246.3\pm 36.7\pm 9.4~{}\mathrm{MeV}.$ (2)

Hadronic molecules attract special attention since many uncommon resonances
are located close to the threshold of a pair of hadrons. For the
$D_{s}\bar{D}_{s1}$ system with $J^{PC}=1^{--}$, Peng et al. Peng:2022nrj
found a bound state at $4.494$ GeV, which is consistent with the mass and
quantum numbers of the $Y(4500)$ observed recently by BESIII BESIII:2022joj .
They have considered the $Y(4230),Y(4360),Y(4500)$ and $Y(4620)$ resonances
from a molecular perspective by formulating a contact-range effective field
theory for the $D\bar{D}_{1}$ and $D_{s}\bar{D}_{s1}$ family of $S$\- and
$P$\- wave charmed meson-antimeson systems. The $Y(4500)$ and $Y(4620)$ are
estimated as hidden-strange partners of the $Y(4230)$ and $Y(4360)$,
respectively.

Additionally, in Ref. Dong:2021juy , charm-strange meson pairs in the energy
range relevant for the $B^{+}\rightarrow J/\psi\phi K^{+}$ is studied, and
$Y(4500)$ is predicted in the threshold $D_{s}(1968)\bar{D}_{s1}(2536)$. This
issue will become clear due to other theoretical analyses and more precise
experiments.

This paper is organized as follows: Detailed information on the QCDSR is given
in Section 2. In Section 3, the results are presented and compared with other
works. Finally, we discuss our findings in Section 4.

##  2 Molecule Assignment for the $Y(4500)$ via QCDSR

QCDSR method is an entirely analytical tool and one of the most active
nonperturbative models in Quantum Chromodynamics (QCD) for obtaining valuable
information on hadrons Shifman:1978bx . The main idea of the QCDSR approach is
to equalize a theoretical expression of a relevant correlator with a
phenomenological one. The underlying concept of this technique is duality,
which builds a relation between hadronic (physical) and quark-gluon (QCD)
degrees of freedom. The sum rules obtained allow us to compute observable
properties of the hadronic state.

The starting point of the QCDSR method is the two-point correlation function,
and it is determined as the vacuum expectation value of the time-ordered
product of interpolation currents as stated below:

$\Pi_{\mu\nu}(q)=i\int d^{4}x~{}e^{iq\cdot
x}\bra{0}T[j_{\mu}(x)j_{\nu}^{\dagger}(0)]\ket{0},$ (3)

where $T$ denotes the time-ordered product operator, $j_{\mu}$ represents the
colourless local current of the $Y$ state with the quantum number
$J^{PC}=1^{--}$ and it is chosen to be Chen:2021kol

$\displaystyle
j_{\mu}(x)=\frac{i}{\sqrt{2}}\bigg{(}\bar{c}_{a}(x)\gamma_{\mu}{\gamma_{5}}{s_{a}(x)}\bar{s}_{b}(x){\gamma_{5}}{c_{b}(x)}-\bar{s}_{a}(x)\gamma_{\mu}{\gamma_{5}}{c_{a}(x)}\bar{c}_{b}(x){\gamma_{5}}s_{b}(x)\bigg{)},$
(4)

where $a,b$ are colour indices. Chosen current $j_{\mu}(x)$ must contain all
the information, like quark contents and quantum numbers of the examined
meson.

###  2.1 Physical (Phenomenological) Side

To extract the mass and decay constant QCDSR expressions, we first compute the
correlation function in terms of the hadronic degrees of freedom. Inserting
the complete set of hadronic states with the same quantum numbers and
performing integral over $x$ in Eq. (3) we get:

$\displaystyle\Pi^{Phys.}_{\mu\nu}(q)=\frac{\bra{0}j_{\mu}(0)\ket{Y(q)}\bra{Y(q)}{j}_{\nu}^{\dagger}(0)\ket{0}}{m^{2}_{Y}-q^{2}}+subtracted~{}terms,$
(5)

where $m_{Y}$ denotes the mass of $Y$ meson. The coupling of the state, Y, to
the current, $j_{\mu}$, can be parametrized in terms of the decay constant
$f_{Y}$ as:

$\bra{0}j_{\mu}(0)Y(q)\rangle=f_{Y}m_{Y}{\varepsilon_{\mu}},$ (6)

where $\varepsilon_{\mu}$ is the polarization vector of the resonance $Y$
satisfying

$\displaystyle\varepsilon_{\mu}\varepsilon_{\nu}^{*}=-g_{\mu\nu}+\frac{q_{\mu}q_{\nu}}{m_{Y}^{2}}.$
(7)

After employing polarization relations, the correlation function is written in
terms of Lorentz structures in the form

$\Pi_{\mu\nu}^{Phys.}(q)=\frac{m_{Y}^{2}f_{Y}^{2}}{m_{Y}^{2}-q^{2}}\left(-g_{\mu\nu}+\frac{q_{\mu}q_{\nu}}{m_{Y}^{2}}\right)+\ldots,$
(8)

here dots represent the contributions originating from the continuum and
higher states. Any chosen Lorentz structure’s coefficient can be employed to
derive the required sum rules. Here, the structure $g_{\mu\nu}$ is selected to
find the sum rules and the typical Borel transformation in terms of $q^{2}$ is
performed to suppress the continuum. In the end, we get the physical side as

$\displaystyle\mathcal{B}\Pi^{Phys.}=m_{Y}^{2}f_{Y}^{2}~{}e^{-m_{Y}^{2}/M^{2}},$
(9)

where $M^{2}$ shows the auxiliary Borel parameter and choosing the scalar
function for the physical side is:

$\Pi^{Phys.}(q^{2})=\frac{m_{Y}^{2}f_{Y}^{2}}{m_{Y}^{2}-q^{2}}.$ (10)

###  2.2 QCD (Theoretical) Side

We represent hadrons by their interpolating quark currents instead of model-
dependent parameters in this part. A correlation function accompanies these
currents and is treated with the operator product expansion (OPE), separating
the short and long-distance quark-gluon interactions. The former is handled by
QCD perturbation theory, while the latter are parametrized according to
universal vacuum condensates.

To evaluate the theoretical part, the current expression in Eq. (4) is
replaced into the correlation function given in Eq. (3). After the heavy and
light quark fields are contracted, we have the correlation function
$\Pi_{\mu\nu}^{\mathrm{QCD}}(q)$ in the molecular picture for the $Y$ state:

$\displaystyle\Pi_{\mu\nu}^{\mathrm{QCD}}(q)$ $\displaystyle=$
$\displaystyle\frac{i}{2}\int d^{4}xe^{iq\cdot
x}\big{\\{}\mathrm{Tr}[\gamma_{5}S_{c}^{bb^{\prime}}(x)\gamma_{5}S_{s}^{b^{\prime}b}(-x)]\mathrm{Tr}[\gamma_{\mu}\gamma_{5}S_{s}^{aa^{\prime}}(x)\gamma_{5}\gamma_{\nu}S_{c}^{a^{\prime}a}(-x)]$
(11) $\displaystyle-$
$\displaystyle\mathrm{Tr}[\gamma_{5}S_{c}^{ba^{\prime}}(x)\gamma_{5}\gamma_{\nu}S_{s}^{a^{\prime}b}(-x)]\mathrm{Tr}[\gamma_{\mu}\gamma_{5}S_{s}^{ab^{\prime}}(x)\gamma_{5}S_{c}^{b^{\prime}a}(-x)]$
$\displaystyle-$
$\displaystyle\mathrm{Tr}[\gamma_{5}S_{s}^{ba^{\prime}}(x)\gamma_{5}\gamma_{\nu}S_{c}^{a^{\prime}b}(-x)]\mathrm{Tr}[\gamma_{\mu}\gamma_{5}S_{c}^{ab^{\prime}}(x)\gamma_{5}S_{s}^{b^{\prime}a}(-x)]$
$\displaystyle+$
$\displaystyle\mathrm{Tr}[\gamma_{5}S_{s}^{bb^{\prime}}(x)\gamma_{5}S_{c}^{b^{\prime}b}(-x)]\mathrm{Tr}[\gamma_{\mu}\gamma_{5}S_{c}^{aa^{\prime}}(x)\gamma_{5}\gamma_{\nu}S_{s}^{a^{\prime}a}(-x)]\\},$

where $S_{c,s}^{ij}(x)$ are the charm and strange quark propagators. These
quark propagators are described in connection with the quark and gluon
condensates Reinders:1984sr  and can be written down in coordinate space using
the light and heavy-quark propagators in Eqs. (12) and (13):

$\displaystyle S_{s}^{ij}(x)$ $\displaystyle=$ $\displaystyle
i\delta_{ij}\frac{\not{x}}{2\pi^{2}x^{4}}-\delta_{ij}\frac{m_{s}}{4\pi^{2}x^{2}}-\delta_{ij}\frac{\langle\overline{s}s\rangle}{12}+i\delta_{ij}\frac{\not{x}m_{s}\langle\overline{s}s\rangle}{48}-\delta_{ij}\frac{x^{2}}{192}\langle\overline{s}g\sigma
Gs\rangle$ (12) $\displaystyle+$ $\displaystyle
i\delta_{ij}\frac{x^{2}\not{x}m_{s}}{1152}\langle\overline{s}g\sigma
Gs\rangle-i\frac{gG_{ij}^{\alpha\beta}}{32\pi^{2}x^{2}}\left[\not{x}{\sigma_{\alpha\beta}+\sigma_{\alpha\beta}}\not{x}\right]-i\delta_{ij}\frac{x^{2}\not{x}g^{2}\langle\overline{s}s\rangle^{2}}{7776}$
$\displaystyle-$
$\displaystyle\delta_{ij}\frac{x^{4}\langle\overline{s}s\rangle\langle
g^{2}G^{2}\rangle}{27648}+\ldots,$ $\displaystyle S_{c}^{ij}(x)$
$\displaystyle=$ $\displaystyle i\int\frac{d^{4}k}{(2\pi)^{4}}e^{-ik\cdot
x}\Bigg{\\{}\frac{\delta_{ij}({\not{k}}+m_{c})}{k^{2}-m_{c}^{2}}-\frac{gG_{ij}^{\alpha\beta}}{4}\frac{\sigma_{\alpha\beta}({\not{k}}+m_{c})+({\not{k}}+m_{c})\sigma_{\alpha\beta}}{(k^{2}-m_{c}^{2})^{2}}$
(13) $\displaystyle+$
$\displaystyle\frac{g^{2}G^{2}}{12}\delta_{ij}m_{c}\frac{k^{2}+m_{c}{\not{k}}}{(k^{2}-m_{c}^{2})^{4}}+\frac{g^{3}G^{3}}{48}\delta_{ij}\frac{({\not{k}}+m_{c})}{(k^{2}-m_{c}^{2})^{6}}[{\not{k}}(k^{2}-3m_{c}^{2})$
$\displaystyle+$ $\displaystyle
2m_{c}(2k^{2}-m_{c}^{2})]({\not{k}}+m_{c})+\ldots\Bigg{\\}},$

where $G_{ij}^{\alpha\beta}\equiv G_{A}^{\alpha\beta}t_{ij}^{A}$ is the
external gluon field strength tensors, $t_{ij}^{A}=\lambda_{ij}^{A}/2$ with
$\lambda_{ij}^{A}$ Gell-Mann matrices, $A=1,2,...8$ shows color indices,
$m_{s}$ and $m_{c}$ symbolizes the strange and charm quark masses
respectively, $\langle\bar{s}s\rangle$ is the strange quark condensate. To get
rid of the contributions emerging from higher states, the standard Borel
transformation is applied to the invariant amplitude $\Pi^{QCD}(q^{2})$,
selecting the same structure in $\Pi^{Phys.}(q^{2})$, one gets the QCD side:

$\Pi^{\mathrm{QCD}}(M^{2},s_{0})=\int_{s_{min}}^{s_{0}}\rho^{\mathrm{QCD}}(s)e^{-s/M^{2}}ds+\Pi(M^{2}),$
(14)

where $s_{min}=4(m_{c}+m_{s})^{2}$ and effective threshold $s_{0}$ shows the
mass of the first excited state with the same quantum number as the selection
of interpolating currents for the considered particle. Based on the
analyticity, spectral density can be described as the imaginary part of the
correlator by the dispersion relation as $\rho^{QCD}(s)=(1/\pi)Im[\Pi^{QCD}]$,
and we evaluate up to operator dimension ten. The spectral density can be
separated according to operator dimensions as

$\displaystyle\rho^{\mathrm{QCD}}(s)$ $\displaystyle=$
$\displaystyle\rho^{\mathrm{pert.}}(s)+\sum_{N=3}^{8}\rho^{DimN}(s),$ (15)

where $\rho^{\mathrm{pert.}}$ denotes the spectral densities of the
perturbative part and the second term implies contributions through the
operator dimensions from three to eight, respectively. We do not give the
explicit forms of spectral density for shortness here. Also, $\Pi(M^{2})$ in
Eq. (14) expresses some non-perturbative contributions:

$\displaystyle\Pi(M^{2})$ $\displaystyle=$
$\displaystyle\sum_{N=6}^{10}\Pi^{DimN}(M^{2}).$ (16)

By matching the coefficients of chosen structure, $g_{\mu\nu}$, in both
physical and QCD sides, and employing the quark hadron duality ansatz, the
final decay constant sum rule for $Y$ meson is extracted as

$f_{Y}^{2}=\frac{e^{m_{Y}^{2}/M^{2}}}{m_{Y}^{2}}\Pi^{QCD}(M^{2},s_{0}).$ (17)

To obtain the final form of the mass and sum rule, the hadronic decay constant
should be removed from the sum rule expression in Eq. (17). Leaving
$\Pi^{QCD}(M^{2},s_{0})$ alone in Eq. (17) and taking the derivative in terms
of $d/d(-1/M^{2})$, the mass sum rule is obtained:

$m_{Y}^{2}=\frac{\Pi^{QCD^{\prime}}(M^{2},s_{0})}{\Pi^{QCD}(M^{2},s_{0})},$
(18)

where $s_{0}$ is the continuum threshold and

$\Pi^{QCD^{\prime}}(M^{2},s_{0})=\frac{d\Pi^{QCD}(M^{2},s_{0})}{d(-1/M^{2})}.$
(19)

##  3 Numerical Analysis

Now, it is time to perform the numerical analysis of the sum rules in Eqs.
(17) and (18). To this aim, the input parameters and condensate values are
presented in Table 1.

Table 1: Input parameters Parameters | Values
---|---
$m_{c}$ |  $1.27\pm 0.02~{}\mathrm{GeV}$ ParticleDataGroup:2022pth
$m_{s}$ |  $93_{3.4}^{+8.6}~{}\mathrm{MeV}$ ParticleDataGroup:2022pth
$\langle\bar{s}s\rangle(1\mbox{GeV})$ |  $0.8\times(-0.24\pm 0.01)^{3}~{}\mathrm{GeV}^{3}$ Shifman:1978bx
$m_{0}^{2}$ |  $(0.8\pm 0.1)$ $\mathrm{GeV}^{2}$ Dosch:1988vv
$\langle\frac{\alpha_{s}}{\pi}G^{2}\rangle$ |  $(0.012\pm 0.004)$ $~{}\mathrm{GeV}^{4}$Shifman:1978bx
$\langle g_{s}^{3}G^{3}\rangle$ |  $(0.57\pm 0.29)$ $~{}\mathrm{GeV}^{6}$Narison:2015nxh

Then, the working regions of the obtained sum rules will be analysed. For
this, we should determine the working interval of the parameters $s_{0}$ and
$M^{2}$. The dominance of pole contribution and convergence of OPE must be
satisfied within the working interval of $s_{0}$ and $M^{2}$. Note that the
sum rules results must be independent of slight variations of these
parameters.

The OPE converges slowly for both tetraquark and molecular states in the QCD
sum rules analysis in the literature. The lower and maximum bounds of the pole
contributions in our calculation require the following condition:

$\mathrm{PC}=\frac{\Pi(M_{\mathrm{max}}^{2},\
s_{0})}{\Pi(M_{\mathrm{max}}^{2},\ \infty)}>(35-73)\%.$ (20)

Applying the mentioned above criteria, the working regions of the parameters
$M^{2}$ and $s_{0}$ are determined as follows:

$\displaystyle
M^{2}(\mathrm{GeV^{2}})\in[5,6],~{}~{}~{}~{}s_{0}(\mathrm{GeV^{2}})\in[21,22].$

Our analysis ensures that the perturbative contribution is dominant over the
non-perturbative ones. In this region, the dimensions six and eight
condensates become significant to compensate for the OPE convergence, as seen
in Figures [2, 2] which is compatible with the selection of
$s_{0}\simeq(m_{Y}+0.15$ GeV$)^{2}$ which is the energy of the first excited
state of the related hadron.

Figure 1: The relative contributions of perturbative and non-perturbative part
of the OPE side for the resonance $Y(4500)$.

Figure 2: The relative contributions of non-perturbative part of the OPE side
in terms of dimensions versus $M^{2}$ for the $Y(4500)$.

The analyses of mass and decay constant of the $Y$ state concerning $M^{2}$
and $s_{0}$ are plotted in Figures [6-6] within these working regions. Both of
them are stable in terms of variations of $M^{2}$ and $s_{0}$.

Figure 3: Mass of resonance $Y(4500)$ as a function of $M^{2}$.

Figure 4: Mass of $Y(4500)$ state versus continuum threshold $s_{0}$.

Figure 5: Decay constant of $Y(4500)$ state as a function of Borel mass
parameter $M^{2}$.

Figure 6: Decay constant of $Y(4500)$ resonance versus continuum threshold
$s_{0}$.

Our results and comparison with the other theoretical studies in the
literature are presented in Table 2. As seen, our mass value nicely agrees
with the BESIII and effective field theory values.

Table 2: The mass and decay constant of resonance $Y(4500)$. | $m_{Y}$ | $f_{Y}$
---|---|---
| $(\mathrm{MeV})$ | $(\mathrm{GeV^{4}})$
Present Work | $4488.35\pm 11.54$ | $(4.04\pm 0.36)\times 10^{-3}$
Experiment | $4484.7\pm 13.3\pm 24.1$ | -
BESIII:2022joj  |  |
Effective Field Theory | $4494.2_{-8.94}^{+5}$ | -
Peng:2022nrj  |  |

##  4 Conclusion

Recent discoveries of $XYZ$ resonances have broadened our perspective of the
hadron spectrum and serve as excellent laboratories to probe exotic
structures. In this paper, we have focused on the charmonium-like $Y(4500)$
state data from the $e^{+}e^{-}\rightarrow K^{+}K^{-}J/\psi$ annihilations. We
calculate the mass and decay constant of this hidden charm-strange meson with
quark composition $[c\bar{s}][\bar{c}s]$, and quantum numbers $J^{PC}=1^{--}$
in the framework of QCDSR. We find that the $Y(4500)$ state can be well
described by the $D_{s}\bar{D}_{s1}$ molecule as suggested in Refs.
Dong:2021juy ; Peng:2022nrj . More data and additional measurements are needed
to provide further information about the nature of the newly observed
structure and charmonium-like states.

## References

  * (1) H. X. Chen, W. Chen, X. Liu, Y. R. Liu and S. L. Zhu, Rept. Prog. Phys. 86, no.2, 026201 (2023)
  * (2) N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C. P. Shen, C. E. Thomas, A. Vairo and C. Z. Yuan, Phys. Rept. 873, 1-154 (2020)
  * (3) X. Liu, Chin. Sci. Bull. 59, 3815-3830 (2014)
  * (4) H. X. Chen, W. Chen, X. Liu and S. L. Zhu, Phys. Rept. 639, 1-121 (2016)
  * (5) S. Agaev, K. Azizi and H. Sundu, Turk. J. Phys. 44, no.2, 95-173 (2020)
  * (6) R. M. Albuquerque, S. Narison, D. Rabetiarivony and G. Randriamanatrika, Nucl. Phys. A 1007, 122113 (2021)
  * (7) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 106, no.1, 014019 (2022)
  * (8) S. S. Agaev, K. Azizi and H. Sundu, [arXiv:2211.14129 [hep-ph]].
  * (9) H. Mutuk, Eur. Phys. J. C 82, no.12, 1142 (2022)
  * (10) S. S. Agaev, K. Azizi and H. Sundu, Eur. Phys. J. Plus 131, no.10, 351 (2016)
  * (11) J. Y. Süngü, A. Türkan and E. Veli Veliev, Acta Phys. Polon. B 50, 1501 (2019)
  * (12) L. Maiani, A. D. Polosa and V. Riquer, Sci. Bull. 66, 1616-1619 (2021)
  * (13) Q. N. Wang, W. Chen and H. X. Chen, Chin. Phys. C 45, no.9, 093102 (2021)
  * (14) J. Y. Süngü, A. Türkan, H. Sundu and E. V. Veliev, Eur. Phys. J. C 82, no.5, 453 (2022)
  * (15) H. Mutuk, Y. Saraç, H. Gümüs and A. Ozpineci, Eur. Phys. J. C 78, no.11, 904 (2018)
  * (16) H. Sundu, SDU J. Nat. Appl. Sci. 20, no.3, 448 (2016)
  * (17) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 95, no.11, 114003 (2017)
  * (18) A. Türkan and H. Dag, Nucl. Phys. A 985, 38-65 (2019)
  * (19) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 106, no.1, 014025 (2022)
  * (20) A. Ali, C. Hambrock and W. Wang, Phys. Rev. D 88, no.5, 054026 (2013)
  * (21) J. Y. Süngü, A. Türkan, H. Dağ and E. Veli Veliev, Adv. High Energy Phys. 2019, 8091865 (2019)
  * (22) K. Azizi and N. Er, Phys. Lett. B 811, 135979 (2020)
  * (23) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 96, no.3, 034026 (2017)
  * (24) Z. G. Wang, X. S. Yang and Q. Xin, Int. J. Mod. Phys. A 36, no.27, 2150202 (2021)
  * (25) R. Tiwari and A. K. Rai, [arXiv:2206.04478 [hep-ph]].
  * (26) X. Liu, H. Huang, J. Ping, D. Chen and X. Zhu, Eur. Phys. J. C 81, no.10, 950 (2021)
  * (27) B. Aubert et al. [BaBar], Phys. Rev. Lett. 95, 142001 (2005)
  * (28) T. E. Coan et al. [CLEO], Phys. Rev. Lett. 96, 162003 (2006)
  * (29) C. Z. Yuan et al. [Belle], Phys. Rev. Lett. 99, 182004 (2007)
  * (30) M. Ablikim et al. [BESIII], Chin. Phys. C 46, no.11, 111002 (2022)
  * (31) M. Ablikim et al. [BESIII], [arXiv:2301.07321 [hep-ex]].
  * (32) F. Z. Peng, M. J. Yan, M. Sánchez Sánchez and M. Pavon Valderrama, Phys. Rev. D 107, no.1, 016001 (2023)
  * (33) X. K. Dong, F. K. Guo and B. S. Zou, Progr. Phys. 41, 65-93 (2021)
  * (34) M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B 147, 385-447 (1979)
  * (35) Z. S. Chen, Z. R. Huang, H. Y. Jin, T. G. Steele and Z. F. Zhang, Chin. Phys. C 46, no.6, 063102 (2022)
  * (36) L. J. Reinders, H. Rubinstein and S. Yazaki, Phys. Rept. 127, 1 (1985)
  * (37) R. L. Workman et al. [Particle Data Group], PTEP 2022, 083C01 (2022)
  * (38) H. G. Dosch, M. Jamin and S. Narison, Phys. Lett. B 220, 251-257 (1989)
  * (39) S. Narison, Nucl. Part. Phys. Proc. 270-272, 143-153 (2016)