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

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.