a_{ij} (\frac{h_je^{u_j}}{\int_Mh_je^{u_j}dV_g}-1)=0,\quad i\in I:=\{1,..,n\}$$ where $dV_g$ is the volume form, $A=(a_{ij})$ is a non-negative constant matrix, $\Delta_g$ is the Laplace-Beltrami operator ($-\Delta_g\ge 0$). When $n=1$ and $a_{11}=1$, equation is the mean field equation of the Liouville type: $$\label{equfromsys}
\Delta_g u+\rho\left(\frac{h e^u}{\int_M h e^udV_g}-1\right)=0\quad\text{in \,}M.$$ Therefore, the Liouville system is a natural extension of the classical Liouville equation, which has been extensively studied for the past three decades. Both the Liouville equation and the Liouville system are related to various fields of geometry, Physics, Chemistry and Ecology. For example in conformal geometry, when $\rho=8\pi$ and $M$ is the sphere $\mathbb S^2$, equation is equivalent to the famous Nirenberg problem. For a bounded domain in $\mathbb R^2$ and $n=1$, a variant of can be derived from the mean field limit of Euler flows or spherical Onsager vortex theory, as studied by Caglioti, Lions, Marchioro and Pulvirenti[@lion1; @lion2], Kiessling[@kiess], Chanillo and Kiessling[@chanillo1] and Lin [@linarch]. In classical gauge field theory, equation is closely related to the Chern-Simons-Higgs equation for the abelian case, see [@caffarelli; @hong; @jackiw; @yang]. Various Liouville systems are also used to describe models in the theory of self-gravitating systems [@aly], Chemotaxis[@childress; @keller], in the physics of charged particle beams[@bennet; @debye; 