@kiess2; @kiessling2], in the non-abelian Chern-Simons-Higgs theory[@dunne; @jostlinwang; @yang] and other gauge field models [@dziar; @phys; @kimleelee]. For recent developments of these subjects or related Liouville systems in more general settings, we refer the readers to [@barto1; @barto3; @ChenLin1; @chenlin2; @chenlinnew; @clwang; @CSW; @CSW1; @licmp; @linduke; @linarch; @lwy; @lwz; @linzhang1; @nolasco2; @rubinstein; @spruck; @wolansky1; @wolansky2; @zhangcmp; @zhangccm] and the references therein.

For any solution $u$ of , clearly adding any constant to $u$ gives another solution. So it is nature to assume $u\in \,{^\text{\r{}}\hspace{-.33cm}H}^1(M)$, where $${^\text{\r{}}\hspace{-.33cm}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\}.$$ Corresponding to (\[mainsys\]) we set $${^\text{\r{}}\hspace{-.33cm}H}^{1,n}=\,{^\text{\r{}}\hspace{-.33cm}H}^1(M)\times\cdots\times\,{^\text{\r{}}\hspace{-.33cm}H}^1(M)$$ to be the space for solutions. For any $\rho=(\rho_1,\cdots, \rho_n)$, $\rho_i>0 (i\in I=\{1,...,n\})$, let $\varPhi_\rho$ be a nonlinear functional defined in ${^\text{\r{}}\hspace{-.33cm}H}^{1,n}$ by $$\varPhi_\rho(u)=\frac{1}{2}\sum_{i,j\in I} a^{ij}\int_M \nabla_g u_i\cdot \nabla_g u_j dV_g-\sum_{j\in I}\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}$. It is easy to see that equation is the Euler-Lagrangian equation of $\varPhi_\rho$.

If the matrix $A$ satisfies the following two assumptions: $$\begin{aligned}
&& (H1):\quad \mbox{ $A$ is symmetric, nonnegative, irreducible and invertible}. \\
&&(H2): \quad a^{ii}\leq 0,\,\, \forall i\in I, \quad a^{ij}\geq 0, \,\, \forall i\ne j,
\quad \sum_{j\in I}a^{ij}\geq 0, \,\, \forall i\in I,\end{aligned}$$ 