the authors prove in [@linzhang2] that for $\rho$ satisfying $$\label{gammaN}
8\pi N\sum_{i\in I}\rho_i<\sum_{i,j\in I}a_{ij}\rho_i\rho_j<8\pi (N+1)\sum_{i\in I}\rho_i,$$ there is a priori estimate for all solutions $u$ to (\[mainsys\]), and the Leray-Schauder degree $d_{\rho}$ for equation (\[mainsys\]) is $$d_{\rho}=
\frac{1}{N!}\bigg ((-\chi_M+1)...(-\chi_M+N) \bigg )\quad \text{if } \mbox{ (\ref{gammaN}) holds }$$ where $\chi_M$ is the Euler characteristic of $M$. Moreover, if $\rho^k$ tends to the hyper-surface $\{\rho;\quad 8\pi N\sum_{i\in I}\rho_i=\sum_{i,j\in I}a_{ij}\rho_i\rho_j\}$, there exist exactly $N$ disjoint blowup points (see [@linzhang2]).

The proof of the a priori bound in [@linzhang2] relies on the sharp estimate for a sequence of bubbling solutions to (\[mainsys\]). Let $u^k$ be the blowup solutions corresponding to $\rho^k$ and $B(p_t,\delta_0)$ ($t=1,..,N$) be disjoint balls around distinct blowup points in $M$. Then under assumptions $(H1)$ and $(H2)$, the behavior of $u^k$ around any $p_t$ is fully bubbling, that is, the maximum values of any components of $u^k$ in any of the balls are of the same magnitude: $$\max_{B(p_t,\delta_0)} u_i^k = \max_{B(p_t,\delta_0)} u_j^k +O(1), \quad \forall i, j\in I.$$ Moreover, after a suitable scaling around each blowup point $p_t$, $u^k$ converges to an entire solution $\tilde U=(\tilde U_1,..,\tilde U_n)$ of the following Liouville system: $$\label{12apr23e2}
\left\{\begin{array}{ll}
\Delta \tilde U_i+\sum_{j=1}^n a_{ij}e^{\tilde U_j}=0,\quad \mbox{in }\,\, 