Let $u_1^k(0)=\max_{B_1}u_1^k$ and $\epsilon_k=e^{-\frac 12u_1^k(0)}$, and $$\label{vikdef}
v_i^k(y)=u_i^k(\epsilon_ky)-u_1^k(0),\quad y\in \Omega_k:=B(0,\epsilon_k^{-1}).$$ Then $v^k=(v_1^k,...,v_n^k)$ converges in $C^2_{loc}(\mathbb R^2)$ to $v=(v_1,..,v_n)$, which satisfies $$\label{817e1}
\left\{\begin{array}{ll}-\Delta v_i=\sum_{j}a_{ij}\tilde h_j(0)e^{v_j},\quad \mathbb R^2,\quad i=1,..,n\\
\\
\int_{\mathbb R^2}e^{v_i}<\infty,\quad i=1,..,n,\quad v_1(0)=0,
\end{array}
\right.$$ where $\tilde h_j(0)=\lim_{k\to \infty}h_j^k(0)$.

For the rest of the paper we set $$m_i:=\frac{1}{2\pi}\int_{\mathbb R^2}\sum_{j=1}^na_{ij} \tilde h_j(0)e^{v_j}>2,\quad i\in I$$ and $m=\min\{m_1,...,m_n\}$. In this section we derive a first order estimate for $v^k$ in $\Omega_k$. In [@linzhang1] the authors prove that there is a sequence of global radial solutions $U^k=(U_1^k,..,U_n^k)$ of (\[817e1\]) such that $$\label{11mar9e1}
 |u_i^k(\epsilon_ky)-U_i^k(y)|\le C,\quad \mbox{for } |y|\le r_0\epsilon_k^{-1}.$$ From (\[11mar9e1\]) we have the following spherical Harnack inequality: $$\label{11mar9e3}
 |u_i^k(\epsilon_ky)-u_i^k(\epsilon_k y')|\le C$$ for all $|y|=|y'|=r\le r_0\epsilon_k^{-1}$ and $C$ is a constant independent of $k,r$. (\[11mar9e3\]) will play an essential role in the first order estimate. To improve (\[11mar9e1\]) we introduce $\phi_i^k$ to be a harmonic function: $$\label{816e5}
\left\{\begin{array}{ll}
-\Delta \phi_i^k=0,\quad B_1,\\
\\
\phi_i^k=u_i^k-\frac 1{2\pi}\int_{\partial B_1}u_i^k,\quad \mbox{on}\quad \partial B_1.
\end{array}
\right.$$ Obviously $\phi_i^k(0)=0$ by the mean value theorem and $\phi_i^k$ is uniformly bounded on $B_{1/2}$ because of (\[88e1\]). Later in section 6, when the results in section 3,4,5 will be used to prove the main theorems, the function $\phi_i^k$ will be specified when we consider the system 