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

In this paper, we continue to consider the generalized Liouville system: $${\mathrm{\Delta }}_g{u}_i+\sum _{j=1}^n{a}_{ij}{\rho }_j(\frac{h_j{e}^{u_j}}{\int h_j{e}^{u_j}}-1)=0\text{in }M,i\in I=\{1,\cdots ,n\},$$ where $(M,g)$ is a Riemann surface $M$ with volume $1$, $h_1,\mathrm{.}\mathrm{.},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 ${\mathrm{\Gamma }}_N$ and proved a priori estimates for $\rho =({\rho }_1,\mathrm{.}\mathrm{.},{\rho }_n)$ in areas separated by ${\mathrm{\Gamma }}_N$. Later Lin-Zhang also calculated the leading term of ${\rho }^k-\rho$ where $\rho\in \Gamma_1$isthelimitof$\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 ${\mathrm{\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 {\mathrm{\Gamma }}_N$, but also reveal new robustness relations of coefficient functions at different blowup points.