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_jleft(frac{h_j e^{u_j}}{int h_j e^{u_j}}- {1} right)=0quadtext{in ,}M,quad iin 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_jin mathbb R^+$($jin 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 $rhoin 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-rhoin Gamma_N$, but also reveal new robustness relations of coefficient functions at different blowup points.