\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{cite} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2015 (2015), No. 76, pp. 1--28.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2015 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2015/76\hfil Mixed interior and boundary peak solutions] {Mixed interior and boundary peak solutions of the Neumann problem for the H\'enon equation \\ in $\mathbb{R}^2$} \author[Y. Zhang, H. Yang \hfil EJDE-2015/76\hfilneg] {Yibin Zhang, Haitao Yang} \address{Yibin Zhang \newline College of Sciences, Nanjing Agricultural University, Nanjing 210095, China} \email{yibin10201029@njau.edu.cn} \address{Haitao Yang \newline Department of Mathematics, Zhejiang University, Hangzhou 310027, China} \email{htyang@css.zju.edu.cn} \thanks{Submitted October 14, 2014. Published March 26, 2015.} \subjclass[2000]{35J25, 35J20, 58K05} \keywords{Mixed interior and boundary peak solutions; H\'enon-type weight; \hfill\break\indent large exponent; Lyapunov-Schmidt reduction process} \begin{abstract} Let $\Omega$ be a bounded domain in $\mathbb{R}^2$ with smooth boundary and $0\in\overline{\Omega}$, we study the Neumann problem for the H\'enon equation \begin{gather*} -\Delta u+u=|x|^{2\alpha}u^p,\quad u>0 \quad \text{in } \Omega,\\ \frac{\partial u}{\partial\nu}=0\quad \text{on } \partial\Omega, \end{gather*} where $\nu$ denotes the outer unit normal vector to $\partial\Omega$, $-1<\alpha\not\in\mathbb{N}\cup\{0\}$ and $p$ is a large exponent. In a constructive way, we show that, as $p$ approaches $+\infty$, such a problem has a family of positive solutions with arbitrarily many interior and boundary spikes involving the origin. The same techniques lead also to a more general result on H\'enon-type weights. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{remark}[theorem]{Remark} \allowdisplaybreaks \section{Introduction} In this article, we consider the Neumann problem \begin{equation}\label{a1} \begin{gathered} -\Delta u+u=S(x)u^p\quad \text{in } \Omega,\\ u>0 \quad\text{in }\Omega,\\ \frac{\partial u}{\partial\nu}=0\quad \text{on }\partial\Omega, \end{gathered} \end{equation} where $\Omega$ is a smooth bounded domain in $\mathbb{R}^2$, $S$ is a nonnegative function on $\overline{\Omega}$, $p$ is a large exponent and $\nu$ denotes the outer unit normal vector to $\partial\Omega$. It is well known that problem \eqref{a1} with $S\equiv1$ has a strong biological meaning because it appears in the study of the stationary Keller-Segel system with the logarithmic sensitivity function from chemotaxis (see \cite{KS,N1}): \begin{equation}\label{a3} \begin{gathered} D_1\Delta v-\chi\nabla\cdot(v\nabla \log \omega)=0\quad\text{in }\Omega,\\ D_2\Delta \omega-a\omega+bv=0\quad\text{in } \Omega,\\ \frac{\partial \omega}{\partial \nu} =\frac{\partial v}{\partial \nu}=0\quad \text{on } \partial\Omega,\\ \frac1{|\Omega|}\int_{\Omega}v(x)dx=\bar{v}>0\quad\text{(prescribed),} \end{gathered} \end{equation} where $\Omega$ is a smooth bounded domain in $\mathbb{R}^N(N\geq2)$ and the constants $D_1$, $D_2$, $a$, $b$ and $\chi$ are positive. Indeed, it is easy to check that solutions of \eqref{a3} satisfy the relation $$ \int_{\Omega}v|\nabla(\log v-p\log\omega)|^2=0\quad\text{where }p=\chi/D_1, $$ so that $v=\lambda\omega^p$ for some positive constant $\lambda$. Thus, setting $\varepsilon^2=D_2/a$, $\gamma=(b\lambda/a)^{\frac1{p-1}}$ and $u=\gamma\omega$, we see that $u$ satisfies the singularly perturbed elliptic problem \begin{equation}\label{a4} \begin{gathered} -\varepsilon^2\Delta u+u=u^p\quad \text{in }\Omega,\\ u>0\quad \text{in } \Omega,\\ \frac{\partial u}{\partial\nu}=0\quad \text{on }\partial\Omega. \end{gathered} \end{equation} In the pioneering papers \cite{LNT,NT1,NT2}, Lin, Ni and Takagi proved that for $\varepsilon>0$ sufficiently small and for $p>1$ subcritical (more precisely, $10\quad \text{in } \Omega,\\ \frac{\partial u}{\partial\nu}=0\quad \text{on }\partial\Omega, \end{gathered} \end{equation} where $\Omega$ is a smooth bounded domain in $\mathbb{R}^2$, $0\in\overline{\Omega}$, $-1<\alpha\not\in\mathbb{N}\cup\{0\}$, $p$ is a large exponent and $\nu$ denotes the outer unit normal vector to $\partial\Omega$. Indeed, the presence of the H\'enon-type weight can produce significant influence on the existence of a solution as well as its asymptotic behavior. For this purpose, this paper is devoted to constructing solutions to problem \eqref{*} with spike-layer profiles at points inside $\Omega$ and on the boundary of $\Omega$ involving the origin when $0\in\overline{\Omega}$ and $p$ tends to $+\infty$. In particular, we recover the result in \cite{MW1} when $S\equiv1$. For \eqref{*} we obtain the following result. \begin{theorem} \label{thm1.1} Assume that $\Omega$ is a smooth bounded domain in $\mathbb{R}^2$ and $0\in\overline{\Omega}$. Then for any $m,\tilde{m}\in\mathbb{N}\cup\{0\}$ and for $p$ sufficiently large problem \eqref{*} has a positive solution $u_p$ which concentrates at $m+\tilde{m}+1$ different points of $\overline{\Omega}$; i.e., as $p\to +\infty$, $$ p|x|^{2\alpha}u^{p+1}_p\rightharpoonup d_1e\delta_{0}+ 8\pi e\sum_{i=2}^{m+1}\delta_{\tilde{q}_i}+ 4\pi e\sum_{i=m+2}^{m+\tilde{m}+1}\delta_{\tilde{q}_i} $$ weakly in the sense of measure in $\overline{\Omega}$, for some points $\tilde{q}_{2},\ldots,\tilde{q}_{m+1}\in\Omega$ and some points $\tilde{q}_{m+2},\ldots,\tilde{q}_{m+\tilde{m}+1} \in\partial\Omega$, where $d_1=8\pi$ for $0\in\Omega$, and $d_1=4\pi$ for $0\in\partial\Omega$. Furthermore, for any $\delta>0$ sufficiently small, \begin{gather*} u_p\to0\quad \text{uniformly in } \overline{\Omega}\setminus \cup_{i=2}^{m+\tilde{m}+1}B_{\delta}(\tilde{q}_i)\cup B_{\delta}(0),\\ \sup_{x\in\overline{\Omega}\cap B_{\delta}(0)}u_p(x)\to\sqrt{e},\quad \sup_{x\in\overline{\Omega}\cap B_{\delta}(\tilde{q}_i)}u_p(x)\to\sqrt{e}, \end{gather*} as $p\to+\infty$. \end{theorem} The above theorem is proved in a constructive way which also works for the more general case involving the H\'enon-type weight with mixed interior and boundary source points as follows: \begin{equation}\label{a2} S(x)=c(x)\prod_{i=1}^{n+\tilde{n}}|x-q_i|^{2\alpha_i}, \end{equation} where $n+\tilde{n}\geq1$, $q_1,\ldots,q_n\in\Omega$, $q_{n+1},\ldots,q_{n+\tilde{n}}\in\partial\Omega$, $-1<\alpha_1,\ldots,\alpha_{n+\tilde{n}} \not\in\mathbb{N}\cup\{0\}$ and $c: \overline{\Omega}\to \mathbb{R}$ is a continuous function satisfying $\inf_{\overline{\Omega}}c>0$, so that problem \eqref{a1} becomes \begin{equation}\label{a} \begin{gathered} -\Delta u+u=c(x)\prod_{i=1}^{n+\tilde{n}}|x-q_i|^{2\alpha_i}u^p\quad \text{in } \Omega,\\ u>0\quad \text{in }\Omega,\\ \frac{\partial u}{\partial\nu}=0 \quad \text{on }\partial\Omega. \end{gathered} \end{equation} Let us first define the corresponding Green's function for the Neumann problem \begin{equation}\label{f} \begin{gathered} -\Delta_x G(x,y)+G(x,y)=\delta_y(x)\quad \text{in }\Omega,\\ \frac{\partial G(x,y)}{\partial \nu_x}=0\quad\text{on }\partial\Omega. \end{gathered} \end{equation} The regular part of $G(x,y)$ is defined depending on whether $y$ lies in the domain or on its boundary as \begin{equation}\label{f1} H(x,y)=\begin{cases} G(x,y)+\frac1{2\pi}\log|x-y|&\text{for} y\in\Omega,\\ G(x,y)+\frac1{\pi}\log|x-y|&\text{for } y\in\partial\Omega. \end{cases} \end{equation} In this way, $H(\cdot,y)$ is of class $C^{1,\beta}$ in $\overline{\Omega}$. Next, for any nonnegative integers $m$ and $\tilde{m}$, we introduce \begin{gather*} J_1=\{1,\ldots,n\},\quad J_2=\{n+1,\ldots,n+\tilde{n}\},\\ J_3=\{n+\tilde{n}+1,\ldots,n+\tilde{n}+m\},\quad J_4=\{n+\tilde{n}+m+1,\dots,n+\tilde{n}+m+\tilde{m}\}. \end{gather*} Furthermore, we set \begin{gather*} \alpha_i=0\quad \text{for }i\in J_3\cup J_4,\\ c_i(x)=\frac{S(x)}{|x-q_i|^{2\alpha_i}}\quad \text{for }i\in\cup_{l=1}^4J_l,\\ d_i=\begin{cases} 8\pi(1+\alpha_i)&\text{for } i\in J_1\cup J_3,\\ 4\pi(1+\alpha_i)&\text{for } i\in J_2\cup J_4. \end{cases} \end{gather*} Also we define \begin{gather*} \Gamma_1=\{q_1,\ldots,q_n\},\quad \Gamma_2=\{q_{n+1},\ldots,q_{n+\tilde{n}}\}, \\ \Lambda_{m}^{\tilde{m}}=(\Omega\setminus\Gamma_1)^m\times (\partial\Omega\setminus\Gamma_2)^{\tilde{m}}\setminus\triangle_{m}^{\tilde{m}}, \end{gather*} where $\triangle_{m}^{\tilde{m}}$ denotes the diagonal set. Now, we fix $n+\tilde{n}$ different source points $q_i$, $i\in J_1\cup J_2$. For $\delta>0$ sufficiently small but fixed we define a configuration space \begin{align*} \Lambda_{m}^{\tilde{m}}(\delta) =\big\{&(q_{n+\tilde{n}+1},\ldots,q_{n+\tilde{n}+m+\tilde{m}})\in \Lambda_{m}^{\tilde{m}}: \operatorname{dist}(q_i,\partial\Omega) \geq 2 \delta\; \forall i\in J_3;\\ &\operatorname{dist}(q_i,q_j)\geq2\delta\; forall i,j\in\cup_{l=1}^4J_l,\;i\neq j \big\}. \end{align*} As a consequence, for points $q=(q_{n+\tilde{n}+1},\ldots,q_{n+\tilde{n}+m+\tilde{m}}) \in\Lambda_{m}^{\tilde{m}}(\delta)$, if we set \begin{equation}\label{b} \begin{aligned} \varphi_{m}^{\tilde{m}}(q) &= \sum_{i\in J_3\cup J_4}d_i\Big\{2\log c_i(q_i)+d_iH(q_i,q_i)+ \sum_{j\in J_1\cup J_2}2d_jG(q_i,q_j)\\ &\quad + \sum_{j\in J_3\cup J_4,\,\,j\neq i}d_jG(q_i,q_j) \Big\}, \end{aligned} \end{equation} we have the following theorem for \eqref{a}, which is the main result of this article. \begin{theorem} \label{thm1.2} Assume that $\Omega$ is a smooth bounded domain in $\mathbb{R}^2$, $n+\tilde{n}\geq1$ and $\inf_{\overline{\Omega}}c>0$. Then for any $m,\tilde{m}\in\mathbb{N}\cup\{0\}$ and for $p$ sufficiently large there exist different points $q^p_l\in\Omega\setminus\Gamma_1$, $l\in J_3$, and $q^p_l\in\partial\Omega\setminus\Gamma_2$, $l\in J_4$, so that \eqref{a} has a positive solution $u_p$ which possesses exactly $n+\tilde{n}+m+\tilde{m}$ local maximum points involving $q_1,\ldots,q_n$, $q^p_{n+\tilde{n}+1},\ldots,q^p_{n+\tilde{n}+m}\in\Omega$, and $q_{n+1},\ldots,q_{n+\tilde{n}}$, $q^p_{n+\tilde{n}+m+1},\ldots,q^p_{n+\tilde{n}+m+\tilde{m}} \in\partial\Omega$. Moreover, $u_p$ has the following concentration property: $$ pS(x)u^{p+1}_p\rightharpoonup e\sum_{l\in J_1\cup J_2}d_l\delta_{q_l}+e\sum_{l\in J_3 \cup J_4}d_l\delta_{\tilde{q}_l} \quad\text{as } p\to+\infty, $$ where $\tilde{q}=(\tilde{q}_{n+\tilde{n}+1},\ldots, \tilde{q}_{n+\tilde{n}+m+\tilde{m}})$ is a global minimum point of $\varphi^{\tilde{m}}_{m}$ in $\Lambda^{\tilde{m}}_{m}(\delta)$ such that for $l\in J_3\cup J_4$, $\operatorname{dist}(q^p_l,\,\tilde{q}_l)\to0$ as $p\to+\infty$. Furthermore, $$ u_p\to0\quad \text{uniformly in } \overline{\Omega}\setminus \cup_{l\in J_1\cup J_2}B_{\delta}(q_l)\cup \cup_{l\in J_3\cup J_4}B_{\delta}(q^p_l), $$ and for the points $q_l$, $l\in J_1\cup J_2$, and $q^p_l$, $l\in J_3\cup J_4$, $$ \sup_{x\in\overline{\Omega}\cap B_{\delta}(q_l)}u_p(x)\to\sqrt{e},\quad \sup_{x\in\overline{\Omega} \cap B_{\delta}(q^p_l)}u_p(x)\to\sqrt{e}, $$ as $p\to+\infty$. \end{theorem} \begin{remark} \label{rmk1.3}\rm The assumption $\inf_{\overline{\Omega}}c>0$ guarantees the existence of global minimum for the function $\varphi^{\tilde{m}}_{m}$ in $\Lambda^{\tilde{m}}_{m}(\delta)$, which follows from properties of the Green's function. The proof is similar to \cite[Lemma 6.1]{MW1}. \end{remark} \begin{remark} \label{rmk1.4}\rm Assume that $\inf_{\overline{\Omega}}c>0$ and $\partial\Omega\setminus\Gamma_2$ has at least one circle, by Ljusternik-Schnirelman theory, we can find another distinct solution satisfying Theorem \ref{thm1.2}. The proof is similar to the one in \cite{DDM}. \end{remark} Note that Theorem \ref{thm1.2} was partly proved in \cite{MW1} only when $S\equiv 1$. Moreover, comparing to the result of \cite{MW1}, this theorem provides a similar but more complex concentration phenomenon involving the presence of mixed interior and boundary peak solutions to \eqref{a}. Indeed, Theorem \ref{thm1.2} implies the existence of solutions for \eqref{a} concentrating at points $q_l$, $l\in J_1\cup J_2$, and $\tilde{q}_l$, $l\in J_3\cup J_4$. Unlike concentration set in \cite{MW1} only contains no source points in the domain and on the boundary, our concentration set also contains some interior source points $q_l$, $l\in J_1$, and boundary source points $q_l$, $l\in J_2$. This in return implies that the presence of mixed interior and boundary source points makes sure that some interior or boundary peak points of solutions of \eqref{a} always locate at these source points. For this reson, if we consider a very simple case of the H\'enon-type weight defined in \eqref{a2}, where \begin{equation} S(x)=|x-0|^{2\alpha}\quad \text{with $0\in\overline{\Omega}$ and $-1<\alpha\not\in\mathbb{N}\cup\{0\}$}, \end{equation} then the corresponding problem \eqref{*} always admits a family of positive solutions with arbitrarily many interior and boundary spikes involving the origin when $p$ tends to $+\infty$, which implies the result in Theorem \ref{thm1.1} hold. Besides, we also point out the interesting result in \cite{D} that solutions for the Liouville equation with the H\'enon-type weight only concentrate at interior points different from the location of the sources. Finally, it is necessary to mention the analogy existing between our results and those known for the Dirichlet problem in $\mathbb{R}^2$: \begin{equation}\label{a5} \begin{gathered} -\Delta u=S(x)u^p\quad \text{in } \Omega,\\ u>0 \quad \text{in } \Omega,\\ u=0 \quad \text{on } \partial\Omega. \end{gathered} \end{equation} Let us point out that \eqref{a5} does not allow any solution with boundary spike-layer profile to exist (see \cite{GNN}), which shows that the Dirichlet boundary condition is far more rigid than the Neumann boundary condition. For $S\equiv1$, asymptotic behavior of least energy solutions of \eqref{a5} is well understood after the works \cite{AG,EG,RW,RW1}: $pu^{p+1}$ approaches a Dirac mass at the harmonic center of $\Omega$ when $p$ tends to infinity. Construction of solutions with this behavior has been achieved in \cite{EMP1}, in which it is shown that for $S\equiv1$, problem \eqref{a5} has solutions with $m$ interior spikes if $\Omega$ is not simply connected. As for the case of the H\'enon-type weight \begin{equation}\label{a6} S(x)=c(x)\prod_{i=1}^n|x-q_i|^{2\alpha_i}, \end{equation} where $n\geq1$, $q_1,\ldots,q_n\in\Omega$, $-1<\alpha_1,\ldots,\alpha_n\not\in\mathbb{N}\cup\{0\}$ and $c: \Omega\to \mathbb{R}$ is a continuous function such that $c(q_i)>0$ for all $i=1,\ldots,n$, related constructions for problem \eqref{a5} have also been performed in \cite{CY2,EPW}, in which it is shown that under a $C^0$-stable critical point assumption there exists a family of positive solutions with exactly $n+m$ interior spikes involving source points $q_1,\ldots,q_n$ as $p$ tends to $+\infty$. The general strategy for proving our main results relies on a Lyapunov-Schmidt reduction procedure, which has appeared in many of the other results mentioned above, as in \cite{DKM,EGP,EMP1,EPW,MW1}. The sketch of this procedure is given as follows: in Section 2 we describe exactly the ansatz for the solution that we are searching for. Then we rewrite problem \eqref{a} in terms of a linearized operator for which a solvability theory, subject to suitable orthogonality conditions, is performed through solving a linearized problem and an intermediate nonlinear problem in Section 3. In Section 4 we reduce problem \eqref{a} to finding critical points of a finite-dimensional function and give its asymptotic expansion. Finally, the proof of Theorem \ref{thm1.2} is contained in Section 5. \section{Ansatz} In this section we describe the approximate solution for \eqref{a} and then we estimate the error of such approximation in appropriate norms. We first fix $n+\tilde{n}$ distinct source points $q_i$, $i\in J_1\cup J_2$, and for $\delta>0$ sufficiently small but fixed we choose points $q=(q_{n+\tilde{n}+1},\ldots,q_{n+\tilde{n}+m+\tilde{m}}) \in\Lambda_{m}^{\tilde{m}}(\delta)$. Moreover, we set \begin{equation} \varepsilon_p=e^{-p/4}, \end{equation} and consider positive numbers $\mu_i$ such that \begin{equation} \delta<\mu_i<\delta^{-1} \quad \forall i\in\cup_{l=1}^4J_l. \end{equation} We define the function \begin{equation}\label{2.0} u_{i}(x)=\log\frac{8(1+\alpha_i)^2\varepsilon_p^2\mu_i^2} {[\varepsilon^2_p\mu_i^2+|x-q_i|^{2(1+\alpha_i)}]^2}, \end{equation} and a correction term as the solution of \begin{equation}\label{2.1} \begin{gathered} -\Delta H_i+H_i=-u_i\quad\text{in }\Omega,\\ \frac{\partial H_i}{\partial \nu}=-\frac{\partial u_i}{\partial \nu} \quad \text{on }\partial\Omega. \end{gathered} \end{equation} \begin{lemma} \label{lem2.1} For any $0<\tau<1/2$ and for any $i\in\cup_{l=1}^4J_l$, then we have \begin{equation}\label{2.6e} H_i(x)=d_iH(x,q_i)+\frac{1}{2}p-\log8(1+\alpha_i)^2\mu_i^2+O(e^{-\tau\frac{p}{4}}), \end{equation} uniformly in $\overline{\Omega}$, where $H$ is the regular part of Green's function defined by \eqref{f1}. \end{lemma} \begin{proof} Note that, on the boundary, we have $$ \frac{\partial H_i}{\partial \nu}=-\frac{\partial u_i}{\partial \nu} =4(1+\alpha_i)|x-q_i|^{2\alpha_i}\frac{(x-q_i)\cdot\nu(x)} {\varepsilon^2_p\mu_i^2+|x-q_i|^{2(1+\alpha_i)}}. $$ Thus, $$ \lim_{p\to\infty}\frac{\partial H_i}{\partial\nu}=4(1+\alpha_i) \frac{(x-q_i)\cdot\nu(x)} {|x-q_i|^2}\quad \forall x\in\partial\Omega\setminus\{q_i\}. $$ On the other hand, by \eqref{f}-\eqref{f1}, the regular part of Green's function $H(x,q_i)$ satisfies \begin{equation}\label{2.11a} \begin{gathered} -\Delta H(x,q_i)+H(x,q_i)=\frac{4(1+\alpha_i)}{d_i}\log |x-q_i|\quad \text{in } \Omega,\\ \frac{\partial H(x,q_i)}{\partial \nu} =\frac{4(1+\alpha_i)}{d_i}\frac{(x-q_i)\cdot\nu(x)} {|x-q_i|^2}\quad \text{on } \partial\Omega. \end{gathered} \end{equation} So, if we set $s_i(x)=H_i(x)-d_i H(x,q_i)-\frac12p+\log8(1+\alpha_i)^2\mu_i^2$, we obtain \begin{gather*} -\Delta s_i+s_i=\log\frac{1}{|x-q_i|^{4(1+\alpha_i)}}- \log\frac1{[\varepsilon_p^2\mu_i^2+|x-q_i|^{2(1+\alpha_i)}]^2}\quad \text{in }\Omega,\\ \frac{\partial s_i}{\partial \nu}=\frac{\partial H_i}{\partial \nu} -4(1+\alpha_i)\frac{(x-q_i)\cdot\nu(x)} {|x-q_i|^2} \quad \text{on }\partial\Omega. \end{gather*} A direct computation shows that for any $\beta\in(1,+\infty)\cap(\frac2{1+\alpha_i},+\infty)$, there exists a positive constant $C$ such that \begin{gather*} \big\|\frac{\partial H_i}{\partial \nu}-4(1+\alpha_i)\frac{(x-q_i)\cdot\nu(x)} {|x-q_i|^2}\big\|_{L^\beta(\partial\Omega)}\leq Ce^{-\frac1{4\beta(1+\alpha_i)}p}, \\ \big\|\log\frac{1}{|x-q_i|^{4(1+\alpha_i)}}- \log\frac1{[\varepsilon_p^2\mu_i^2+|x-q_i|^{2(1+\alpha_i)}]^2} \big\|_{L^\beta(\Omega)}\leq Cpe^{-\frac1{2\beta(1+\alpha_i)}p}. \end{gather*} By $L^\beta$ theory $$ \|s_{i}\|_{W^{1+s,\beta}(\Omega)} \leq C(\|\Delta s_i\|_{L^{\beta}(\Omega)}+\|\frac{\partial s_i}{\partial \nu}\|_{L^\beta(\partial\Omega)}) \leq Ce^{-\frac1{4\beta(1+\alpha_i)}p} $$ for $00 \quad \text{in } \Omega_{\varepsilon},\\ \frac{\partial v}{\partial\nu}=0 \quad \text{on }\partial\Omega_{\varepsilon}. \end{gathered} \end{equation} Let us define the first approximation solution of \eqref{2.13} as \begin{equation}\label{2.13a} V_q(y)=\varepsilon_p^{\frac2{p-1}}U_q(\varepsilon_p y), \end{equation} where $U_q$ is defined by \eqref{2.6}. For $|\varepsilon_py-q_i|<\delta$ with $\delta$ sufficiently small but fixed, by using Lemmas \ref{lem2.1}--\ref{lem2.2} and \eqref{2.6a}, we have \begin{align*} %\label{2.7b} &V_q(y)\\ &= \frac1{p^{\frac{p}{p-1}}\mu_i^{\frac2{p-1}} c_i(q_i)^{\frac1{p-1}}}[u_i+H_i+\frac1p(\omega_{0i}+H_{0i}) +\frac1{p^2}(\omega_{1i}+H_{1i})] \\ &\quad +\sum_{j\neq i} \frac1{p^{\frac{p}{p-1}}\mu_j^{\frac2{p-1}} c_j(q_j)^{\frac1{p-1}}}[u_j+H_j+\frac1p(\omega_{0j}+H_{0j}) +\frac1{p^2}(\omega_{1j}+H_{1j})] \\ &=\frac1{p^{\frac{p}{p-1}}\mu_i^{\frac2{p-1}} c_i(q_i)^{\frac1{p-1}}}\Big\{ \big[U^i(\frac{x-q_i}{v_i\rho_i})+\frac12p-\log\mu_i^2\big]+ \big[d_iH(x,q_i)+\frac{1}{2}p \\ &\quad -\log8(1+\alpha_i)^2\mu_i^2\big] +\frac1p\big[\omega_{0i}-\frac{d_iC_{0i}}{4(1+\alpha_i)}H(x,q_i) +\frac{C_{0i}}{1+\alpha_i}\log\mu_i-\frac{pC_{0i}}{4(1+\alpha_i)}\big] \\ &\quad +\frac1{p^2}\big[\omega_{1i}-\frac{d_iC_{1i}}{4(1+\alpha_i)}H(x,q_i) +\frac{C_{1i}}{1+\alpha_i}\log\mu_i-\frac{pC_{1i}}{4(1+\alpha_i)}\big] +O(\varepsilon_p^\tau)\Big\}\\ &\quad +\sum_{j\neq i} \frac{d_jG(x,q_j)}{p^{\frac{p}{p-1}}\mu_j^{\frac2{p-1}}c_j(q_j)^{\frac1{p-1}}} \big[1-\frac{C_{0j}}{4p(1+\alpha_j)}-\frac{C_{1j}}{4p^2(1+\alpha_j)} +O(\varepsilon_p^\tau)\big] \\ &=\frac1{p^{\frac{p}{p-1}}\mu_i^{\frac2{p-1}}c_i(q_i)^{\frac1{p-1}}} \Big\{p+U^i(\frac{x-q_i}{v_i\rho_i})+\frac{1}p\omega_{0i} +\frac{1}{p^2}\omega_{1i}+O(|x-q_i|)+O(\varepsilon^\tau_p) \\ &\quad -\log8(1+\alpha_i)^2\mu_i^4+d_i \big[1-\frac{C_{0i}}{4p(1+\alpha_i)} -\frac{C_{1i}}{4p^2(1+\alpha_i)}\big]H(q_i,q_i) \\ &\quad + \frac{1}{p(1+\alpha_i)}\big(C_{0i}+\frac{1}pC_{1i}\big) \big(\log\mu_i-\frac14p\big) \\ &\quad +\sum_{j\neq i}\big[\frac{\mu_i^2c_i(q_i)}{\mu_j^2c_j(q_j)} \big]^{\frac1{p-1}}d_j \big[1-\frac{C_{0j}}{4p(1+\alpha_j)}-\frac{C_{1j}}{4p^2(1+\alpha_j)}\big] G(q_i,q_j)\Big\}. \end{align*} We now choose the parameters $\mu_i$: we assume they are defined by the relation \begin{align*}\label{2.6b} \log8(1+\alpha_i)^2\mu_i^4 &=d_i \big[1-\frac{C_{0i}}{4p(1+\alpha_i)}-\frac{C_{1i}}{4p^2(1+\alpha_i)}\big] H(q_i,q_i) \\ &\quad + \frac{1}{p(1+\alpha_i)}\big(C_{0i}+\frac{1}pC_{1i}\big) \big(\log\mu_i-\frac14p\big) \\ &\quad +\sum_{j\neq i}\big[\frac{\mu_i^2c_i(q_i)}{\mu_j^2c_j(q_j)} \big]^{\frac1{p-1}}d_j \big[1-\frac{C_{0j}}{4p(1+\alpha_j)}-\frac{C_{1j}}{4p^2(1+\alpha_j)}\big]G(q_i,q_j). \end{align*} Taking into account the explicit expression \eqref{2.5b} of the constant $C_{0i}$, we observe that $\mu_i$ bifurcates, as $p$ tends to $+\infty$, from the value \begin{equation}\label{2.6bbc} \bar{\mu}_i=e^{-\frac34}e^{\frac14[d_iH(q_i,q_i) +\sum_{j\neq i}d_jG(q_j,q_i)]}, \end{equation} solution of the equation \begin{equation}\label{2.6c} \log8(1+\alpha_i)^2\mu_i^4=d_iH(q_i,q_i) -\frac{C_{0i}}{4(1+\alpha_i)}+\sum_{j\neq i}d_jG(q_i,q_j). \end{equation} Thus, $\mu_i$ is a perturbation of order $\frac1p$ of the value $\bar{\mu}_i$, namely \begin{equation}\label{2.6bb} \mu_i=e^{-\frac34}e^{\frac14[d_iH(q_i,q_i) +\sum_{j\neq i}d_jG(q_j,q_i)]}\left(1+O(\frac1p)\right). \end{equation} From this choice of the parameters $\mu_i$, we deduce that for $|\varepsilon_py-q_i|<\delta$, \begin{equation}\label{2.7c} V_q(y) =\frac1{p^{\frac{p}{p-1}}\mu_i^{\frac2{p-1}}c_i(q_i)^{\frac1{p-1}}}\big[ p+U^i(z_i)+\frac{1}p\omega^{0i}(z_i)+\frac{1}{p^2}\omega^{1i}(z_i) +\theta(y)\big], \end{equation} with $$ z_i=\frac{1}{v_i\rho_i}(\varepsilon_py-q_i),\quad \theta(y)=O(\rho_i|z_i|)+O(\varepsilon^\tau_p). $$ In the rest of this paper, we will seek solutions of problem \eqref{2.13} in the form $v=V_q+\phi$, where $\phi$ will represent a lower order correction. In terms of $\phi$, problem \eqref{2.13} becomes \begin{equation}\label{2.23} \begin{gathered} L(\phi):=-\Delta \phi+\varepsilon_p^2\phi-W_q\phi=R_q+N(\phi)\quad \text{in }\Omega_{\varepsilon},\\ \frac{\partial \phi}{\partial\nu}=0 \quad \text{on }\partial\Omega_{\varepsilon}, \end{gathered} \end{equation} where \begin{equation}\label{2.13b} W_q=pS(\varepsilon_p y)V_q^{p-1}, \end{equation} $R_q$ is the ``error term'' \begin{equation}\label{2.13bb} R_q=\Delta V_q-\varepsilon_p^2V_q+S(\varepsilon_p y)V_q^p, \end{equation} and $N(\phi)$ stands for the ``nonlinear term'' \begin{equation}\label{2.24} N(\phi)=S(\varepsilon_p y) [(V_q+\phi)^p-V^p_q-pV_q^{p-1}\phi]. \end{equation} The main step in solving problem \eqref{2.23} for small $\phi$ is, under a suitable choice of the points $q_i$, that of a solvability theory for the linear operator $L$. In developing this theory, we introduce a weighted $L^\infty$-space $$ \mathcal {C}_{*}=\{h\in L^{\infty}(\Omega_\varepsilon): \|h\|_{*}<+\infty\}, $$ with the norm \begin{equation}\label{3.2a} \|h\|_{*}=\sup_{y\in\overline{\Omega}_{\varepsilon}} \frac{|h(y)|}{\varepsilon_p^2+\sum_{\alpha_i<0,\,i\in J_1\cup J_2}\frac{(\frac{\varepsilon_p}{v_i\rho_i})^2 |z_i|^{2\alpha_i}}{(1+|z_i|)^{4+2\alpha+2\alpha_i}}+\sum_{\alpha_i\geq 0,\,i\in\cup_{l=1}^{4}J_{l}} \frac{(\frac{\varepsilon_p}{v_i\rho_i})^2}{(1+|z_i|)^{4+2\alpha}}}\,, \end{equation} where we fix $-1<\alpha<\alpha_0$ with $\alpha_0=\min\{\alpha_i: i\in\cup_{l=1}^4J_l\}$. With respect to the norm \eqref{3.2a}, the error term $R_q$ defined in \eqref{2.13bb} can be estimated in the following way. \begin{lemma} \label{lem2.3} For $\delta>0$ sufficiently small but fixed, there exist $C$, $D_0>0$ and $p_0>1$ such that \begin{gather}\label{2.22} \|R_q\|_*\leq C/p^4, \\ \label{2.19} W_q(y)\leq D_0\sum_{i=1}^{n+\tilde{n}+m+\tilde{m}} (\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}e^{U^i(z_i)}, \end{gather} for any $q\in\Lambda_{m}^{\tilde{m}}(\delta)$ and $p\geq p_0$. Furthermore, \begin{equation}\label{2.18c} W_q(y)=(\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}e^{U^i(z_i)} \big[1+\frac{\omega^{0i}-U^i-\frac12(U^i)^2}p +O(\frac{\log^4(2+|z_i|)}{p^2})\big], \end{equation} for any $|z_i|\leq\delta(v_i\rho_i)^{-1/2}$, where $z_i=\frac{\varepsilon_py-q_i}{v_i\rho_i}$. \end{lemma} \begin{proof} Observe that by \eqref{2.1}, \begin{gather*} -\Delta (u_i+H_i)+\varepsilon_p^2(u_i+H_i) =(\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}e^{U^i(z_i)} \quad \text{in } \Omega_\varepsilon,\\ \frac{\partial}{\partial \nu}(u_i+H_i)=0 \quad \text{on } \partial\Omega_\varepsilon, \end{gather*} and by \eqref{2.3} and \eqref{2.2}, for $j=0,1$, \begin{gather*} -\Delta(\omega_{ji}+ H_{ji})+\varepsilon_p^2(\omega_{ji}+H_{ji}) =(\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i} (e^{U^i}\omega^{ji}-f^{ji})(z_i)\quad \text{in }\Omega_\varepsilon,\\ \frac{\partial}{\partial \nu}(\omega_{ji}+H_{ji})=0\quad \text{on }\partial\Omega_\varepsilon, \end{gather*} which combined with the definition of $V_q$ in \eqref{2.13a} can deduce directly that \begin{equation}\label{2.13c} \begin{aligned} &-\Delta V_q+\varepsilon_p^2V_q\\ &=\sum_{i=1}^{n+\tilde{n}+m+\tilde{m}} \frac{(\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}} {p^{\frac p{p-1}}\mu_i^{\frac2{p-1}} c_i(q_i)^{\frac1{p-1}}} \big[e^{U^i}+\frac{e^{U^i}\omega^{0i}-f^{0i}}p +\frac{e^{U^i}\omega^{1i}-f^{1i}}{p^2}\big](z_i). \end{aligned} \end{equation} From the definition \eqref{2.0a} of $U^i$ and the asymptotic behavior \eqref{2.5a} of $\omega^{ji}$, $j=0,1$, we obtain that in a region far away from the points $q_i$, namely for $|x-q_i|>\delta$; i.e. $|z_i|>\delta(v_i\rho_i)^{-1}$ for all $i$, $$ U^i(z_i)=-p+O(1),\quad \omega^{ji}(z_i)=\frac{C_{ji}}{4(1+\alpha_i)}p+O(1), $$ and so \begin{equation}\label{2.13d} \big[e^{U^i}+\frac{e^{U^i}\omega^{0i}-f^{0i}}p +\frac{e^{U^i}\omega^{1i}-f^{1i}}{p^2}\big](z_i)=O(p^2e^{-p}). \end{equation} This, together with \eqref{2.6a} and \eqref{2.13a} imply \begin{equation}\label{2.21} R_q=O(pe^{-p})+O(p^{-p-1}). \end{equation} Let us now fix the index $i$ and the region $|z_i|\leq\delta(v_i\rho_i)^{-1/2}$. From \eqref{2.7c} we obtain \begin{equation}\label{2.7a} \begin{aligned} S(\varepsilon_p y)V_q^p &=\frac{(\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}} {p^{\frac{p}{p-1}}\mu_i^{\frac{2}{p-1}} c_i(q_i)^{\frac{1}{p-1}}}(1+O(\rho_i|z_i|)) \Big[1+\frac{U^i(z_i)}{p}\\ &\quad +\frac{\omega^{0i}(z_i)}{p^2} +\frac{\omega^{1i}(z_i)}{p^3} +\frac{\theta(y)}{p}\Big]^p. \end{aligned} \end{equation} By the Taylor expansions of the exponential and logarithmic function, we have for $|z_i|\leq\delta(v_i\rho_i)^{-1/2}$, \begin{equation}\label{2.7f} \begin{aligned} &\Big(1+\frac{a}p+\frac{b}{p^2}+\frac{c}{p^3}\Big)^p\\ &=e^a\big[1+\frac1p\big(b-\frac{a^2}2\big)+\frac1{p^2}\big(c-ab+\frac{a^3}3 +\frac{b^2}2+\frac{a^4}8-\frac{a^2b}2\big) \\ &\quad +O(\frac{\log^6(|z_i|+1)}{p^3})\big], \end{aligned} \end{equation} provided $-5(1+\alpha_i)\log(|z_i|+2)\leq a(z_i)\leq C$ and $|b(z_i)|+|c(z_i)|\leq C\log(|z_i|+2)$. Thus, for $|z_i|\leq\delta(v_i\rho_i)^{-1/2}$, \begin{equation}\label{2.7h} \begin{aligned} S(\varepsilon_p y)V_q^p &=\frac{(\frac{\varepsilon_p}{v_i\rho_i})^2 |z_i|^{2\alpha_i}}{p^{\frac p{p-1}}\mu_i^{\frac2{p-1}} c_i(q_i)^{\frac1{p-1}}}\big[e^{U^i}+\frac{e^{U^i}\omega^{0i}-f^{0i}}p +\frac{e^{U^i}\omega^{1i}-f^{1i}}{p^2}\big](z_i)\\ &\quad + (\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}e^{U^i(z_i)} O\big(\frac{\log^6(|z_i|+2)}{p^4} +\varepsilon^\tau_p+p^2\rho_i|z_i|\big). \end{aligned} \end{equation} Joining together \eqref{2.13c}, \eqref{2.13d} and \eqref{2.7h}, in this region we obtain \begin{equation}\label{2.21b} R_q=(\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}e^{U^i(z_i)} O(\frac{\log^6(|z_i|+2)}{p^4} +\varepsilon^\tau_p+p^2\rho_i|z_i|)+O(pe^{-p}). \end{equation} On the other hand, if $\delta(v_i\rho_i)^{-1/2}\leq|z_i|\leq\delta(v_i\rho_i)^{-1}$, we have \begin{equation}\label{2.17a} -\Delta V_q+\varepsilon_p^2V_q =O(p(\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}e^{U^i(z_i)})+O(pe^{-p}), \end{equation} and by \eqref{2.7a}, \begin{equation}\label{2.7e} S(\varepsilon_p y)V_q^p=O(\frac1p(\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i} e^{U^i(z_i)}), \end{equation} since $(1+\frac{s}p)^s\leq e^s$. Thus, in this region \begin{equation}\label{2.21a} R_q=O(p(\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}e^{U^i(z_i)}) +O(pe^{-p}). \end{equation} So, estimate \eqref{2.22} can be easily derived from \eqref{2.21}, \eqref{2.21b} and \eqref{2.21a}. Finally, to prove the estimate over $W_q(y)=pS(\varepsilon_p y)V_q^{p-1}$, we first notice a slight modification of formula \eqref{2.7f} $$ \Big(1+\frac{a}p+\frac{b}{p^2}+\frac{c}{p^3}\Big)^{p-1} =e^a\big[1+\frac1p\big(b-a-\frac{a^2}2\big)+ +O(\frac{\log^4(|z_i|+2)}{p^2})\big]. $$ Thus, for $|z_i|\leq\delta(v_i\rho_i)^{-1/2}$, by \eqref{2.7c}, \begin{align*} W_q(y) &=pS(\varepsilon_p y)V_q^{p-1}\\ &=\frac{c_i(\varepsilon_p y)|\varepsilon_p y-q_i|^{2\alpha_i}}{\mu_i^2c_i(q_i)} \big[1+\frac{U^i}p+\frac{\omega^{0i}}{p^2} +\frac{\omega^{1i}}{p^3}+O(\frac{\rho_i|z_i|+\varepsilon^\tau_p}{p})\big]^{p-1}\\ &=(\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}e^{U^i(z_i)} \big[1+\frac{\omega^{0i}-U^i-\frac12(U^i)^2}p +O(\frac{\log^4(2+|z_i|)}{p^2})\big]. \end{align*} In addition, if $|z_i|\geq\delta(v_i\rho_i)^{-1}$ for all $i$, we obtain that $W_q(y)=O(p^{1-p})$, and if $\delta(v_i\rho_i)^{-1/2}<|z_i|<\delta(v_i\rho_i)^{-1}$, $W_q(y)=O((\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}e^{U^i(z_i)})$, which completes the proof. \end{proof} \begin{remark} \label{rmk2.4} \rm As for $W_q$, let us point out that if $|z_i|\leq\delta(v_i\rho_i)^{-1}$ for some $i$, $$ pS(\varepsilon_p y)\Big(V_q+O(\frac1{p^3})\Big)^{p-2} =O\Big((\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}e^{U^i(z_i)}\Big). $$ Since this estimate is true if $|z_i|>\delta(v_i\rho_i)^{-1}$ for all $i\in\cup_{l=1}^4J_l$, we have \begin{equation}\label{2.21e} pS(\varepsilon_p y)\Big(V_q+O(\frac1{p^3})\Big)^{p-2} \leq C\sum_{i=1}^{n+\tilde{n}+m+\tilde{m}} (\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}e^{U^i(z_i)}. \end{equation} \end{remark} \section{The linearized problem and the nonlinear problem} In this section we shall study first the following linear problem: given $h$ of class $\mathcal {C}_{*}$ and points $q \in\Lambda_{m}^{\tilde{m}}(\delta)$, we find a function $\phi$, scalars $c_{ij},\,i\in J_3$, $j=1,2$, and $c_{i1},\,i\in J_4$, such that \begin{equation}\label{3.3} \begin{gathered} L(\phi)=h+\sum_{i\in J_3\cup J_4}c_{i1}\chi_i\,Z_{i1} +\sum_{i\in J_3}c_{i2}\chi_i\,Z_{i2}\quad \text{in }\Omega_{\varepsilon},\\ \frac{\partial\phi}{\partial\nu}=0\quad \text{on } \partial\Omega_{\varepsilon},\\ \int_{\Omega_{\varepsilon}}\chi_i Z_{ij}\phi=0\quad \forall i\in J_3,\; j=1,2;\; i\in J_4,\; j=1, \end{gathered} \end{equation} where $\chi_i$ and $Z_{ij}$ are defined as follows. \begin{gather*} L_i(\phi)=\Delta\phi+\frac{8(1+\alpha_i)^2|z|^{2\alpha_i}}{[1+|z|^{2(1+\alpha_i)} ]^2}\phi \quad \forall i\in\cup_{l=1}^4J_l,\quad z_{i1}=\frac{4\operatorname{Re}(z)}{|z|^2+1} \quad \forall i\in J_3\cup J_4,\\ z_{i0}=\frac1{\mu_i}\frac{|z|^{2(1+\alpha_i)}-1}{|z|^{2(1+\alpha_i)}+1}\quad \forall i\in\cup_{l=1}^4J_l\quad z_{i2}=\frac{4\operatorname{Im}(z)}{|z|^2+1} \quad \forall i\in J_3\cup J_4. \end{gather*} It is well known that any bounded solution to $L_i(\phi)=0$ in $\mathbb{R}^2$ is a linear combination of $z_{i0}$, $z_{i1}$ and $z_{i2}$ for $i\in J_3\cup J_4$ (see \cite{BP,CL5}), or proportional to $z_{i0}$ for $i\in J_1\cup J_2$ (see \cite{CY1,E,EPW}). Then for $i\in J_2\cup J_4$, any bounded solution to $L_i(\phi)=0$ in $\mathbb{R}_+^2$ and $\frac{\partial}{\partial z_2}\phi(z_1,0)=0$ on $\partial\mathbb{R}_+^2$ is a linear combination of $z_{i0}$ and $z_{i1}$ if $i\in J_4$, or proportional to $z_{i0}$ if $i\in J_2$. For each point $q_i\in\partial\Omega$ with $i\in J_2\cup J_4$, we have to strengthen the boundary. Let us assume that $q_i=0$ and the unit outward normal at $q_i$ is $(0,-1)$. Let $G(x_1)$ be the defining function for the boundary $\partial\Omega$ in a neighborhood $B_\delta(q_i)$ of $q_i$, that is, $\Omega\cap B_\delta(q_i)=\{(x_1,x_2)\in B_\delta(q_i): x_2>G(x_1)\}$. Then let $F_i: B_\delta(q_i)\cap\overline{\Omega}\to \mathbb{R}^2$ be defined by $F_i=(F_{i1}, F_{i2})$, where $$ F_{i1}=x_1+\frac{x_2-G(x_1)}{1+|G'(x_1)|^2}G'(x_1) \quad\text{and} \quad F_{i2}=x_2-G(x_1). $$ We consider $\tilde{z}_i=F_i^\varepsilon(y)=\frac{1}{v_i\rho_i}F_i(\varepsilon_p y)$ and its inverse $y=\tilde{F}^\varepsilon_i(\tilde{z}_i) =\frac1{\varepsilon_p}F_i^{-1}(v_i\rho_i\tilde{z}_i)$. Set $$ \tilde{z}_i(y):=\begin{cases} z_i(y) & \forall i\in J_1\cup J_3,\\ F_i^\varepsilon(y) &\forall i\in J_2\cup J_4. \end{cases} $$ Furthermore, set \begin{gather*} Z_{i0}(y):= z_{i0}(\tilde{z}_i)\quad \forall i\in\cup_{l=1}^4J_l,\quad Z_{i1}(y):= z_{i1}(\tilde{z}_i)\quad \forall i\in J_3\cup J_4,\\ Z_{i2}(y):=z_{i2}(\tilde{z}_i)\quad \forall i\in J_3,\quad \chi_i(y):=\chi(|\tilde{z}_i|)\quad \forall i\in\cup_{l=1}^4J_l, \end{gather*} where $\chi(r)$ is a smooth, non-increasing cut-off function such that for a large but fixed number $R_0>0$, $\chi(r)=1$ if $r\leq R_0$, and $\chi(r)=0$ if $r\geq R_0+1$. It is important to note that $F_i$, $i\in J_2\cup J_4$, preserves the Neumann boundary condition and \begin{equation}\label{3.2b} \Delta_yZ_{i0}+W_{i0}Z_{i0}=\begin{cases} 0 & \forall \,i\in J_1\cup J_3,\\ O(\frac{\varepsilon_p^2}{v_i\rho_i}\frac{|z_i|^{2\alpha_i}}{(1+|z_i|)^{3+4\alpha_i}}) & \forall i\in J_2\cup J_4, \end{cases} \end{equation} where $$ W_{i0}:=\big(\frac{\varepsilon_p}{v_i\rho_i}\big)^2 \frac{8(1+\alpha_i)^2|z_i|^{2\alpha_i}} {[1+|z_i|^{2(1+\alpha_i)}]^2}. $$ The main result of this section is the following: \begin{proposition} \label{prop3.1} There exist $C>0$ and $p_0>1$ such that for any $p\geq p_0$, $h\in\mathcal {C}_{*}$ and $q\in\Lambda_{m}^{\tilde{m}}(\delta)$, there exists a unique solution $\phi\in L^{\infty}(\Omega_{\varepsilon})$, scalars $c_{ij},\,i\in J_3$, $j=1,2$, and $c_{i1},\,i\in J_4$, to \eqref{3.3}. Moreover, such solution satisfies \begin{equation}\label{3.5} \|\phi\|_{\infty}\leq Cp\|h\|_{*}. \end{equation} \end{proposition} The proof will be divided into a series of lemmas which we state and prove next. \begin{lemma} \label{lem3.2} For $p$ large enough, there exist $R_1>0$, and $$ \psi: \overline{\Omega}_{\varepsilon}\setminus \cup_{i=1}^{n+\tilde{n}+m+\tilde{m}}B_{i,R_1}\,\,\longmapsto\,\mathbb{R} $$ positive and uniformly bounded so that \begin{equation}\label{3.7a} \begin{gathered} L(\psi)\geq\sum_{i\in\cup_{l=1}^{4}J_l} (\frac{\varepsilon_p}{v_i\rho_i})^2\frac{1}{|z_i|^{4+2\alpha}}+\varepsilon_p^2 \quad \text{in }\Omega_{\varepsilon} \setminus\cup_{i=1}^{n+\tilde{n}+m+\tilde{m}}B_{i,R_1},\\ \frac{\partial}{\partial\nu}\psi\geq\sum_{i\in J_2\cup J_4} \frac{\varepsilon_p}{v_i\rho_i}\frac{1}{|z_i|^{3+2\alpha}}+\varepsilon_p \quad \text{on } \partial\Omega_{\varepsilon}\setminus\cup_{i\in J_2\cup J_4}B_{i,R_1}, \end{gathered} \end{equation} where $-1<\alpha<\alpha_0$ and $B_{i,R_1}=\{y\in\overline{\Omega}_{\varepsilon}:|z_i|0$, $r_i=|z_i|$ and $R_b=\frac{1}{b}3^{\frac{1}{2(1+\alpha)}}$. Then for $R_b\leq r_i\leq\delta(v_i\rho_i)^{-1}$, $$ -\Delta g_{1i}(y)+\varepsilon_p^2g_{1i}(y) >\big(\frac{\varepsilon_p}{v_i\rho_i}\big)^2 b^{-2(1+\alpha)}\frac{4(1+\alpha)^2}{r_i^{2\alpha+4}}, $$ and by \eqref{2.19}, $$ W_q(y)\leq C_1\big(\frac{\varepsilon_p}{v_i\rho_i}\big)^2\frac1{r_i^{2\alpha_{i}+4}}. $$ So, if $b>0$ is sufficiently small so that $4(1+\alpha)^2b^{-2(1+\alpha)}>C_1+1$, we have \begin{equation}\label{3.8d} L(g_{1i})(y)>\big(\frac{\varepsilon_p}{v_i\rho_i}\big)^2\frac{1}{r_i^{2\alpha+4}}, \end{equation} and for $i\in J_2\cup J_4$, \begin{equation}\label{3.8e} |\frac{\partial }{\partial \nu}g_{1i}(y)|\leq \frac{\varepsilon_p}{v_i\rho_{i}}\frac{C_{2}\delta} {r_i^{3+2\alpha}}. \end{equation} We also take $$ g_{2i}(y)=1-\frac{1}{r_i^{3+2\alpha}}\operatorname{Im}(z_i)\quad \text{for }i\in J_2\cup J_4,\; R_b\leq r_i\leq\delta(v_i\rho_i)^{-1}. $$ Then \begin{gather}\label{3.8f} L(g_{2i})(y)\geq (\frac{\varepsilon_p}{v_i\rho_i})^2 \big\{(1+2\alpha)(3+2\alpha)\frac{1}{r_i^{5+2\alpha}} \operatorname{Im}(z_i)-C_1\frac{1}{r_i^{4+2\alpha_{i}}}\big\}, \\ \label{3.8g} \frac{\partial}{\partial \nu}g_{2i}(y) \geq\frac{\varepsilon_p}{v_i\rho_i}\frac{C_{3}}{r_i^{3+2\alpha}}. \end{gather} Consider now $$ g_{3i}(y)=k_{1}g_{1i}(y)+g_{2i}(y)\quad \text{for } i\in J_2\cup J_4,\,\,R_b\leq r_i\leq\delta(v_i\rho_i)^{-1}, $$ where $k_1$ is chosen larger and $\delta$ is chosen smaller if necessary. It is easy to check that $|g_{3i}(y)|\leq C_{4}$ and \[ |\nabla g_{3i}(y)|=O(\frac{\varepsilon_p}{v_i\rho_i}\frac{1}{r_i^{3+2\alpha}}). \] By \eqref{3.8d}-\eqref{3.8g}, we find \begin{gather}\label{3.8h} L(g_{3i})(y)>(\frac{\varepsilon_p}{v_i\rho_i})^2\frac{1}{r_i^{4+2\alpha}}, \\ \label{3.8i} \frac{\partial}{\partial \nu}g_{3i}(y) >\frac{\varepsilon_p}{v_i\rho_i}\frac{C_{5}}{r_i^{3+2\alpha}}. \end{gather} Let $\eta_i(y)=\eta(v_i\rho_ir_i)$, where $\eta=\eta(t)$ is a smooth cut-off function in $\mathbb{R}^2$ such that $\eta=1$ if $t\leq\frac{1}{2}\delta$, $\eta=0$ if $t\geq\delta$. Obviously, for $\delta(2v_i\rho_i)^{-1}\leq r_i\leq\delta(v_i\rho_i)^{-1}$, $|\nabla\eta_i(y)|\leq C_6\varepsilon_p$ and $|\Delta\eta_i(y)|\leq C_6\varepsilon^2_p$. Besides, let $g_{0}(y)=\tilde{g}(\varepsilon_p y)$, where $\tilde{g}$ is a bounded positive solution of \begin{gather*} -\Delta\tilde{g}+\tilde{g}=1\quad \text{in } \Omega,\\ \frac{\partial\tilde{g}}{\partial\nu}=1 \quad \text{on }\partial\Omega, \end{gather*} so that $-\Delta g_0+\varepsilon_p^2g_0=\varepsilon_p^2$ in $\Omega_{\varepsilon}$, $\frac{\partial}{\partial\nu}g_{0}=\varepsilon_p$ on $\partial\Omega_{\varepsilon}$. Moreover, $g_{0}(y)$ is uniformly bounded on $\overline{\Omega}_{\varepsilon}$, i.e. $|g_{0}(y)|\leq C_7$. Thus, for numbers $k_2,k_3$ and $R_1$ such that $k_2>\max\{2,C_{5}^{-1}\}$, $k_3>2+k_2C_6+k_2C_4C_6$ and $R_{1}=\max\{R_b,(k_{3}C_{1}C_{7})^{\frac{1}{2(\alpha_{i}-\alpha)}}|\,i\in \cup_{l=1}^{4}J_l\}$. The function $$ \psi(y)=\sum_{i\in J_1\cup J_3}k_2\eta_{i}(y)g_{1i}(y) +\sum_{i\in J_2\cup J_4}k_{2}\eta_{i}(y)g_{3i}(y)+k_{3}g_{0}(y) $$ meets the requirements. The rest of the proof is similar to that of \cite[Lemma 3.4]{CY1}, so we omit it. \end{proof} \begin{lemma} \label{lem3.3} There exist $C>0$ and $p_0>1$ such that for all $p>p_0$, points $q\in\Lambda^{\tilde{m}}_{m}(\delta)$, $h\in\mathcal {C}_{*}$, and $\phi$ the solution to \begin{equation}\label{3.3r} \begin{gathered} L(\phi)=-\Delta \phi+\varepsilon_p^2\phi-W_q\phi=h\quad \text{in }\Omega_{\varepsilon},\\ \frac{\partial\phi}{\partial\nu}=0\quad \text{on } \partial\Omega_{\varepsilon}, \end{gathered} \end{equation} under the orthogonality conditions \begin{equation}\label{3.7} \begin{gathered} \int_{\Omega_{\varepsilon}}\chi_i\,Z_{i0}\phi=0 \quad \forall i\in J_1\cup J_2,\\ \int_{\Omega_{\varepsilon}}\chi_i Z_{ij}\phi=0 \quad \forall i\in J_3,\;j=0,1,2;\; i\in J_4,\; j=0,1, \end{gathered} \end{equation} one has \begin{equation}\label{3.8} \|\phi\|_{\infty}\leq C\|h\|_{*}. \end{equation} \end{lemma} \begin{proof} First we consider the ``inner norm'' $$ \|\phi\|_{R}=\sup_{y\in\overline{\Omega}_{\varepsilon} \cap(\cup_{i=1}^{n+\tilde{n}+m+\tilde{m}}B_{i,R})}|\phi(y)|, $$ we claim that there is a constant $C>0$ such that \begin{equation}\label{3.9} \|\phi\|_{\infty}\leq C(\|\phi\|_{R_1}+\|h\|_{*}), \end{equation} with $R_{1}$ given by Lemma \ref{lem3.2}. Set $$ \tilde{\phi}=C_1\psi(y)(\|\phi\|_{R_1}+\|h\|_{*}) \quad\forall y\in\overline{\Omega}_{\varepsilon}\setminus \cup_{i=1}^{n+\tilde{n}+m+\tilde{m}}B_{i,R_1}, $$ where $\psi$ is the positive, uniformly bounded barrier constructed by Lemma \ref{lem3.2} and $C_1>0$ is chosen larger if necessary. Then for $y\in\Omega_{\varepsilon}\setminus\cup_{i=1}^{n+\tilde{n}+m+\tilde{m}}B_{i,R_1}$, \begin{align*} L(\tilde{\phi}\pm\phi)(y) &\geq C_{1}\,\|h\|_{*}\Big\{\sum_{i\in\cup_{l=1}^{4} J_l}(\frac{\varepsilon_p}{v_i\rho_i})^2|\frac{\varepsilon_p y-q_i}{v_i\rho_i}|^{-(4+2\alpha)}+\varepsilon_p^2\Big\}\pm h(y)\\ &\geq|h(y)|\pm h(y)\geq0, \end{align*} for $y\in\partial\Omega_{\varepsilon}\setminus\cup_{i=1}^{n+\tilde{n}+m+\tilde{m}}B_{i,R_1}$, $$ \frac{\partial}{\partial\nu}(\tilde{\phi}\pm\phi)(y)\geq 0, $$ and for $y\in\overline{\Omega}_{\varepsilon}\cap(\cup_{i=1}^{n+\tilde{n} +m+\tilde{m}}\partial B_{i,R_1})$, $$ (\tilde{\phi}\pm\phi)(y)>\|\phi\|_{R_1}\pm\phi(y)\geq |\phi(y)|\pm\phi(y)\geq 0. $$ From the maximum principle (see \cite{PW}), it follows that $-\tilde{\phi}(y)\leq\phi(y)\leq\tilde{\phi}(y)$ on $\overline{\Omega}_{\varepsilon}\setminus \cup_{i=1}^{n+\tilde{n}+m+\tilde{m}}B_{i,R_1}$, which implies that \eqref{3.9} holds. Let us prove the a priori estimate \eqref{3.8} by contradiction. We assume the existence of sequences $p_j\to+\infty$, functions $h_j$ with $\|h_j\|_{*}\to 0$, points $q^j=(q^j_{n+\tilde{n}+1},\ldots,q^j_{n+\tilde{n}+m+\tilde{m}}) \in\Lambda^{\tilde{m}}_{m}(\delta)$, solutions $\phi_j$ with $\|\phi_j\|_\infty=1$ such that \eqref{3.3r}-\eqref{3.7} holds. From estimate \eqref{3.9}, $\sup_{\overline{\Omega}_{\varepsilon}\cap B_{i,R_1}}|\phi_j|\geq\kappa>0$ for some $i$ and $\kappa$. To simplify the notation, let us set $p=p_j$, $\varepsilon_p=\varepsilon_{p_j}$ and $q_i=q_i^j$. Set $\hat{\phi}_j(z)=\phi_j(\frac{v_i\rho_i}{\varepsilon_p}z+\frac{q_i}{\varepsilon_p})$ and $\hat{h}_j(z)=h_j(\frac{v_i\rho_i}{\varepsilon_p}z+\frac{q_i}{\varepsilon_p})$. While $q_i\in\Omega$, $i\in J_1\cup J_3$, by \eqref{2.18c}, $\hat{\phi}_j$ satisfies $$ -\Delta\hat{\phi}_j+(v_i\rho_i)^2\hat{\phi}_j -\frac{8(1+\alpha_i)^2|z|^{2\alpha_i}}{[1+|z|^{2(1+\alpha_i)}]^2} [1+O(\frac1p)\,]\hat{\phi}_j=(\frac{v_i\rho_i}{\varepsilon_p})^2\,\widehat{h}_j, $$ for $z\in B_{R_1}(0)$. Obviously, for any $\beta\in [1,-\frac1{\alpha}]$, $(\frac{v_i\rho_i}{\varepsilon_p})^2\hat{h}_j\to 0$ in $L^\beta(B_{R_1}(0))$. Since $\frac{8(1+\alpha_{i})^2|z|^{2\alpha_{i}}}{[1+|z|^{2(1+\alpha_i)}]^2}$ is bounded in $L^\beta(B_{R_1}(0))$ and $\|\hat{\phi}_j\|_{\infty}=1$, elliptic regularity theory readily implies that $\hat{\phi}_j$ converges uniformly over compact subsets near the origin to a bounded nontrivial solution $\hat{\phi}$ of $L_i(\hat{\phi})=0$ in $\mathbb{R}^2$. Then $\hat{\phi}$ is proportional to $z_{i0}$ for $i\in J_1$, or a linear combination of $z_{i0}$, $z_{i1}$ and $z_{i2}$ for $i\in J_3$. However, our assumed orthogonality conditions \eqref{3.7} on $i$ and $\phi_j$, pass to limit and yield $\int\chi z_{i0}\hat{\phi}dz=0$ for $i\in J_1$, or $\int\chi z_{il}\hat{\phi}dz=0$ for $i\in J_3,\,l=0,1,2$, which implies $\hat{\phi}\equiv0$. This is absurd because $\hat{\phi}$ is nontrivial. While $q_i\in\partial\Omega,\,i\in J_2\cup J_4$, $\hat{\phi}_j$ satisfies \begin{gather*} -\Delta\hat{\phi}_j+(v_i\rho_i)^2\hat{\phi}_j -\frac{8(1+\alpha_i)^2|z|^{2\alpha_i}}{[1+|z|^{2(1+\alpha_i)}]^2} [1+O(\frac1p)]\hat{\phi}_j=\big(\frac{v_i\rho_i}{\varepsilon_p}\big)^2\\ \hat{h}_j\quad \text{in } B_{R_1}(0)\cap\Omega_{v_i\rho_i},\\ \frac{\partial}{\partial\nu}\hat{\phi}_j=0\quad \text{on } B_{R_1}(0)\cap\partial\Omega_{v_i\rho_i}, \end{gather*} where $\Omega_{v_i\rho_i}=\frac{1}{v_i\rho_i}(\Omega-\{q_i\})$. Obviously, for any $\beta\in[1,-\frac1{\alpha}]$, $(\frac{v_i\rho_i}{\varepsilon_p})^2\hat{h}_j\to 0$ in $L^\beta(B_{R_1}(0)\cap\Omega_{v_i\rho_i})$. Since $\frac{8(1+\alpha_{i})^2|z|^{2\alpha_i}}{[1+|z|^{2(1+\alpha_i)}]^2}$ is bounded in $L^\beta(B_{R_1}(0)\cap\Omega_{v_i\rho_i})$ and $\|\hat{\phi}_j\|_{\infty}=1$, elliptic regularity theory readily implies that $\hat{\phi}_j$ converges uniformly on $B_{R_1}(0)\cap\overline{\Omega}_{v_i\rho_i}$ to a bounded nontrivial solution $\hat{\phi}$ of $L_i(\hat{\phi})=0$ in $\mathbb{R}^2_+$ and $\frac{\partial }{\partial z_2}\hat{\phi}(z_1,0)=0$ on $\partial\mathbb{R}^2_+$. Then $\hat{\phi}$ is proportional to $z_{i0}$ for $i\in J_2$, or a linear combination of $z_{i0}$ and $z_{i1}$ for $i\in J_4$. However, our assumed orthogonality conditions \eqref{3.7} on $i$ and $\phi_j$, pass to limit and easily yield $\int\chi z_{i0}\hat{\phi}=0$ for $i\in J_2$, or $\int\chi z_{il}\hat{\phi}dz=0$ for $i\in J_4,\,l=0,1$, which implies $\hat{\phi}\equiv0$. This is also absurd because $\hat{\phi}$ is nontrivial, which completes the proof. \end{proof} \begin{lemma} \label{lem3.4} There exist $C>0$ and $p_0>1$ such that for any $p>p_0$, points $q\in\Lambda^{\tilde{m}}_{m}(\delta)$, $h\in\mathcal {C}_{*}$, and solution $\phi$ of \eqref{3.3r} with \begin{equation}\label{3.10} \int_{\Omega_{\varepsilon}}\chi_iZ_{ij}\phi=0 \quad \forall i\in J_3,\; j=1,2;\; i\in J_4,\;j=1, \end{equation} one has \begin{equation}\label{3.11} \|\phi\|_{\infty}\leq Cp\|h\|_{*}. \end{equation} \end{lemma} \begin{proof} Let $R>R_0+1$ be large and fixed. For $i\in \cup_{l=1}^{4}J_l$, denote \begin{gather*} a_{i0}=\frac1{\mu_i[H(q_i,q_i)-\frac{1}{d_i}4(1+\alpha_i)\log(v_i\rho_i R)]}, \\ \widehat{Z}_{i0}(y)=Z_{i0}(y)-\frac1{\mu_i} +a_{i0}G(q_i,\varepsilon_p y). \end{gather*} Let $\eta_1$ and $\eta_2$ be radial smooth cut-off functions in $\mathbb{R}^2$ such that \begin{gather*} 0\leq\eta_1\leq1;\quad |\nabla\eta_1|\leq C\,\,\text{in}\,\mathbb{R}^2; \quad \eta_1\equiv1\text{ in }B_R(0);\quad \eta_1\equiv0 \text{ in } \mathbb{R}^2\setminus B_{R+1}(0);\\ 0\leq\eta_2\leq1;\quad |\nabla\eta_2|\leq C\text{ in }\mathbb{R}^2;\quad \eta_2\equiv1\text{ in } B_{\frac{\kappa}{4}\delta}(0);\quad \eta_2\equiv0\text{ in }\mathbb{R}^2\setminus B_{\frac{\kappa}3\delta}(0), \end{gather*} where $0<\kappa\leq1$ will be chosen later on. Set \begin{gather*} \eta_{1i}(y)=\eta_1(|\tilde{z}_i(y)|),\quad \eta_{2i}(y)=\eta_2(v_i\rho_i|\tilde{z}_i(y)|), \\ \widetilde{Z}_{i0}(y)=\eta_{1i}Z_{i0}+(1-\eta_{1i})\eta_{2i}\widehat{Z}_{i0}. \end{gather*} Now define \begin{equation}\label{3.11t} \widetilde{\phi}=\phi+\sum_{i=1}^{n+\tilde{n}+m+\tilde{m}}e_i\widetilde{Z}_{i0}, \end{equation} where $e_i$ is chosen such that for any $i\in\cup_{l=1}^4J_l$, $$ e_i\int_{\Omega_{\varepsilon}}\chi_i|Z_{i0}|^2 +\int_{\Omega_{\varepsilon}}\chi_iZ_{i0}\phi=0, $$ Note that $\frac{\partial}{\partial\nu}\widetilde{Z}_{i0}=0$ for any $i\in J_2\cup J_4$, which arises from $\frac{\partial}{\partial z_2}z_{i0}(z_1,0)=0$ and $F_i$ preserves the Neumann boundary condition. Thus \begin{equation}\label{3.12} \begin{gathered} L(\widetilde{\phi})=h+\sum_{i=1}^{n+\tilde{n}+m +\tilde{m}}e_i\,L(\widetilde{Z}_{i0}) \quad\text{in }\Omega_{\varepsilon},\\ \frac{\partial\widetilde{\phi}}{\partial\nu}=0 \quad \text{on } \partial\Omega_{\varepsilon}, \end{gathered} \end{equation} and $\widetilde{\phi}$ satisfies the orthogonality conditions \eqref{3.7}. By \eqref{3.8}, it implies that \begin{equation}\label{3.13} \|\widetilde{\phi}\|_\infty\leq\, C\big\{\|h\|_*+\sum_{i=1}^{n+\tilde{n}+m +\tilde{m}}|e_i|\cdot\|L(\widetilde{Z}_{i0})\|_{*}\big\}. \end{equation} Multiplying \eqref{3.12} by $\widetilde{Z}_{i0}$, $i\in\cup_{l=1}^4J_l$, and integrating by parts, it follows that \begin{equation}\label{3.134} \langle L(\widetilde{Z}_{i0}),\widetilde{\phi}\rangle =\langle\widetilde{Z}_{i0},h\rangle+e_i\langle L(\widetilde{Z}_{i0}),\widetilde{Z}_{i0}\rangle, \end{equation} where $\langle f,g\rangle=\int_{\Omega_\varepsilon}fg$. Then for $i\in\cup_{l=1}^{4}J_l$, by \eqref{3.13}-\eqref{3.134}, $$ |e_i\langle L(\widetilde{Z}_{i0}),\widetilde{Z}_{i0}\rangle| \leq C\|h\|_{*}(1+\|L(\widetilde{Z}_{i0})\|_{*}) +C\|L(\widetilde{Z}_{i0})\|_{*}\sum_{i=1}^{n+\tilde{n}+m +\tilde{m}}|e_j|\cdot\|L(\widetilde{Z}_{j0})\|_{*}. $$ Let us claim that for any $i\in\cup_{l=1}^{4}J_l$, \begin{gather} \|L(\widetilde{Z}_{i0})\|_{*}=O(\frac1p), \label{3.15}\\ \langle L(\widetilde{Z}_{i0}),\widetilde{Z}_{i0}\rangle =-\frac{C}{p}\big\{1+O(\frac{\log R}{R^{2(1+\alpha_i)}})\big\}.\label{3.16} \end{gather} Once these estimates are proven, it easily follows that \begin{equation}\label{3.17} |e_i|\leq Cp\,\|h\|_{*}. \end{equation} This, together with \eqref{3.11t}, \eqref{3.13} and \eqref{3.15}, implies estimate \eqref{3.11} for $\phi$. \noindent \textbf{Proof of \eqref{3.15}.} For sake of simplicity, we only prove estimate $\eqref{3.15}$ for any $i\in J_2\cup J_4$. Since $\nabla F_i(q_i)=Id$, it follows that \begin{gather} z_i=\frac{\varepsilon_p y-q_i}{v_i\rho_i} =\frac{\varepsilon_p\tilde{F}_i^\varepsilon(\tilde{z}_i)-q_i}{v_i\rho_i} =\tilde{z}_i\{1+O(v_i\rho_i\tilde{z}_i)\}, \label{3.53a}\\ \tilde{z}_i=F_i^\varepsilon(y) =\frac{F_i(v_i\rho_i z_i+q_i)}{v_i\rho_i}=z_i\{1+O(v_i\rho_iz_i)\}, \label{3.53a1}\\ \Delta_y=(\frac{\varepsilon_p}{v_i\rho_i})^2\{\Delta_{\tilde{z}_i} +O(\rho_i|\tilde{z}_i|) \nabla_{\tilde{z}_i}^2+O(\rho_i)\nabla_{\tilde{z}_i}\}. \label{3.53c} \end{gather} So, $\kappa\leq1$ can be chosen such that if $|\tilde{z}_i|\leq\frac13\kappa\delta(v_i\rho_i)^{-1}$, then $|z_i|\leq\delta(v_i\rho_i)^{-1}$. Note that \begin{equation}\label{3.14} \begin{aligned} L(\widetilde{Z}_{i0}) &=\eta_{1i}LZ_{i0}+(1-\eta_{1i})\eta_{2i}L\widehat{Z}_{i0} +2\nabla\eta_{1i}\nabla Z_{i0}+Z_{i0}\Delta\eta_{1i}\\ & \quad +2\nabla\widehat{Z}_{i0}\nabla[(1-\eta_{1i})\eta_{2i}] +\widehat{Z}_{i0}\Delta[(1-\eta_{1i})\eta_{2i}]. \end{aligned} \end{equation} For $r_i:=|\tilde{z}_i|\leq R+1$, by \eqref{2.18c}, \eqref{3.2b} and \eqref{3.53a}, $$ L(Z_{i0})=(W_q-W_{i0})Z_{i0}+O(\frac{\varepsilon_p^2}{v_i\rho_i}\frac{ |z_i|^{2\alpha_i}}{(1+|z_i|)^{3+4\alpha_i}})-\varepsilon_p^2Z_{i0}, $$ which implies \begin{equation}\label{3.15b} \|\eta_{1i}L(Z_{i0})\|_{*}=O(\frac1p). \end{equation} Note that for $ r_i\geq R$, by \eqref{3.53a}-\eqref{3.53a1}, \begin{gather} |Z_{i0}-\frac1{\mu_i}|\leq\,C(1+|\tilde{z}_i|)^{-2(1+\alpha_i)} =O(|z_i|^{-2(1+\alpha_i)}), \label{3.15c}\\ |a_{i0}G(q_i,\varepsilon_p y)-\frac1{\mu_i}| \leq \frac{C(1+\log\frac{|z_i|}R)} {H(q_i,q_i)-\frac1{d_i}4(1+\alpha_i)\log(v_i\rho_i R)}=O (\frac{1+\log\frac{|z_i|}R}p).\label{3.15d} \end{gather} For $R\leq r_i\leq\frac13\kappa\delta(v_i\rho_i)^{-1}$, by \eqref{3.2b}, $$ L(\widehat{Z}_{i0})=W_q[a_{i0}G(q_i,\varepsilon_p y)-\frac1{\mu_i}] +[W_q-W_{i0}]Z_{i0}-\varepsilon_p^2(Z_{i0}-\frac1{\mu_i}) +O(\frac{\varepsilon^2_p}{v_i\rho_i|z_i|^{3+2\alpha_i}}), $$ which, together with \eqref{2.19}, \eqref{3.15c} and \eqref{3.15d}, easily indicates that \begin{equation}\label{3.15e} \|(1-\eta_{1i})\eta_{2i}L(\widehat{Z}_{i0})\|_{*}=O(\frac1p). \end{equation} For $R\leq r_i\leq R+1$, \begin{equation}\label{3.15j} Z_{i0}-\widehat{Z}_{i0} =\frac1{\mu_i}-a_{i0}G(q_i,\varepsilon_p y)=\frac{H(q_i,q_i) -H(q_i,\varepsilon_p y)+ \frac{4(1+\alpha_i)}{d_i}\log\frac{|z_i|}R}{\mu_i[H(q_i,q_i) -\frac{4(1+\alpha_i)}{d_i}\log(v_i\rho_i R)]}, \end{equation} which, together with \eqref{3.53c}, implies that \begin{gather}\label{3.15f} (Z_{i0}-\widehat{Z}_{i0})\Delta\eta_{1i}= O(\frac1p)(\frac{\varepsilon_p}{v_i\rho_i})^2\{\eta''_1(r_i) +\frac1{r_i}\eta'_1(r_i)+O(\rho_i)\}, \\ \label{3.15g} \nabla\eta_{1i}\nabla(Z_{i0}-\widehat{Z}_{i0}) =-(\frac{\varepsilon_p}{v_i\rho_i})^2 \frac{\eta_1'}{\mu_i\log(v_i\rho_iR)}\frac1{r_i}\,\{1+O(\frac1p)\}. \end{gather} Thus, \begin{equation}\label{3.15a} \|2\nabla\eta_{1i}\nabla(Z_{i0}-\widehat{Z}_{i0})+(Z_{i0} -\widehat{Z}_{i0})\Delta\eta_{1i}\|_{*}=O(\frac1p). \end{equation} For $\frac14\kappa\delta(v_i\rho_i)^{-1}\leq r_i \leq \frac13\kappa\delta(v_i\rho_i)^{-1}$, it follows that $\widehat{Z}_{i0}=O(\frac1p)$ and $\nabla\widehat{Z}_{i0}=O(\frac{\varepsilon_p}p)$. Then \begin{equation}\label{3.15i} \|\widehat{Z}_{i0}\Delta\eta_{2i}+2\nabla\widehat{ Z}_{i0}\nabla\eta_{2i}\|_{*} =O(\frac1p). \end{equation} Joining \eqref{3.14}, \eqref{3.15b}, \eqref{3.15e}, \eqref{3.15a} and \eqref{3.15i}, we obtain that \eqref{3.15} holds for $i\in J_2\cup J_4$. \smallskip \noindent\textbf{Proof of \eqref{3.16}.} Note that for any $i\in\cup_{l=1}^{4}J_l$, $\langle L(\widetilde{Z}_{i0}),\widetilde{Z}_{i0}\rangle=I+K$, where \begin{align*} I&=\int_{\Omega_\varepsilon}\widetilde{Z}_{i0} \{\eta_{1i}L Z_{i0}+(1-\eta_{1i})\eta_{2i}L\widehat{Z}_{i0}\}\\ &=\int_{\Omega_\varepsilon}\eta_{2i}^2\big\{Z_{i0}+(1-\eta_{1i}) [a_{i0}G(q_i,\varepsilon_p y)-\frac1{\mu_i}]\big\} \Big\{(W_q-W_{i0})Z_{i0}+(\Delta Z_{i0}\\ &\quad +W_{i0}Z_{i0}) +(1-\eta_{1i})W_q[a_{i0}G(q_i,\varepsilon_p y)-\frac1{\mu_i}] -\varepsilon_p^2(Z_{i0}-\frac1{\mu_i}+\frac{1}{\mu_i}\eta_{1i})\Big\}, \end{align*} and \[ K=\int_{\Omega_{\varepsilon}}\widetilde{Z}_{i0}\{Z_{i0}\Delta\eta_{1i} +\widehat{Z}_{i0}\Delta[(1-\eta_{1i})\eta_{2i}]+2\nabla Z_{i0}\nabla\eta_{1i} +2\nabla\widehat{Z}_{i0}\nabla[(1-\eta_{1i})\eta_{2i}]\}. \] For $R\leq r_i=|\tilde{z}_i|\leq\frac13\kappa\delta(v_i\rho_i)^{-1}$, by \eqref{2.19}, \eqref{3.53a}, \eqref{3.53a1} and \eqref{3.15d}, \begin{align*} &\int_{\Omega_{\varepsilon}}\eta_{2i}^2(1-\eta_{1i})\{Z_{i0}+(1-\eta_{1i})[a_{i0}G(q_i,\varepsilon_p y)-\frac1{\mu_i}]\} W_q[a_{i0}G(q_i,\varepsilon_p y)-\frac1{\mu_i}]\\ &=O(\frac{\log R}{pR^{2(1+\alpha_i)}}). \end{align*} By \eqref{3.2b} and \eqref{3.15d}, $$ \int_{\Omega_{\varepsilon}}\eta_{2i}^2\{Z_{i0}+(1-\eta_{1i}) [a_{i0}G(q_i,\varepsilon_p y)-\frac1{\mu_i}]\}(\Delta Z_{i0} +W_{i0}Z_{i0})=O(\rho_i+\varepsilon_p^2). $$ By \eqref{3.15c}-\eqref{3.15d}, $$ \varepsilon_p^2\int_{\Omega_{\varepsilon}}\eta_{2i}^2\{Z_{i0} +(1-\eta_{1i}) [a_{i0}G(q_i,\varepsilon_p y)-\frac1{\mu_i}]\}(Z_{i0}-\frac1{\mu_i}) =O(\rho_i^2+\varepsilon_p^2). $$ Furthermore, \begin{equation}\label{3.16d} I=\int_{\Omega_{\varepsilon}}\eta_{2i}^2(W_q-W_{i0})Z^2_{i0} dy +O(\frac{\log R}{pR^{2(1+\alpha_i)}})+O(\rho_i+\varepsilon_p^2). \end{equation} By \eqref{2.5c}, a straightforward but tedious computation shows that \begin{gather}\label{3.16f} \int_{\mathbb{R}^2}(\frac{\varepsilon_p}{v_i\rho_i})^2|z_i|^{2\alpha_i}e^{U^i(z_i)} \frac{\omega^{0i}-U^i-\frac12(U^i)^2}pz^2_{i0}(z_i)dy =-\frac{8\pi(\alpha_i+1)}{p\mu_i^2}, \\ \label{3.16g} \int_{|z_i|\geq\delta(v_i\rho_i)^{-1/2}}(\frac{\varepsilon_p}{v_i\rho_i})^2 |z_i|^{2\alpha_i}e^{U^i(z_i)} \frac{\omega^{0i}-U^i-\frac12(U^i)^2}pz^2_{i0}(z_i)dy =O(\frac1p\varepsilon_p^{\frac12}). \end{gather} By \eqref{3.53a}-\eqref{3.53a1}, it easily follows that for $|z_i|\leq\delta(v_i\rho_i)^{-1/2}$, \begin{equation}\label{3.16g1} Z_{i0}(y)=z_{i0}(\tilde{z}_i)=z_{i0}(z_i)[1+O(\rho_i|z_i|)]. \end{equation} This, together with \eqref{2.18c} and \eqref{3.16f}-\eqref{3.16g1}, easily indicates that \begin{equation}\label{3.16h} \int_{|z_i|\leq\delta(v_i\rho_i)^{-1/2}}\eta_{2i}^2(W_q-W_{i0})Z^2_{i0}= \begin{cases} -\frac{8\pi(\alpha_i+1)}{p\mu_i^2} +O(\frac1{p^2}) & \forall \,i\in J_1\cup J_3,\\ -\frac{4\pi(\alpha_i+1)}{p\mu_i^2} +O(\frac1{p^2}) & \forall \,i\in J_2\cup J_4. \end{cases} \end{equation} Also by \eqref{2.19}, \begin{equation}\label{3.16e} \int_{|z_i|\geq\delta(v_i\rho_i)^{-1/2}}\eta_{2i}^2(W_q-W_{i0})Z^2_{i0} =\int_{|z_i|\geq\delta(v_i\rho_i)^{-1/2}}O(\frac1{|z_i|^{4+2\alpha_i}}) =O(\varepsilon_p). \end{equation} So, by \eqref{3.16d}, \eqref{3.16h} and \eqref{3.16e}, it follows that \begin{equation}\label{3.16i} I=\begin{cases} -\frac{8\pi(\alpha_i+1)}{p\mu_i^2} +O(\frac{\log R}{pR^{2(1+\alpha_i)}}) &\forall \,i\in J_1\cup J_3,\\ -\frac{4\pi(\alpha_i+1)}{p\mu_i^2} +O(\frac{\log R}{pR^{2(1+\alpha_i)}}) &\forall \,i\in J_2\cup J_4. \end{cases} \end{equation} Let us give the estimate of $K$. Integrating by parts the first two terms of $K$, \begin{align*} K&=\int_{\Omega_{\varepsilon}} \widetilde{Z}_{i0}\nabla Z_{i0}\nabla\eta_{1i} -Z_{i0}\nabla\widetilde{Z}_{i0}\nabla\eta_{1i} +\widetilde{Z}_{i0}\nabla\widehat{Z}_{i0}\nabla[(1-\eta_{1i})\eta_{2i}]\\ &\quad -\widehat{Z}_{i0}\nabla\widetilde{Z}_{i0}\nabla[(1-\eta_{1i})\eta_{2i}]\\ &=K_1+K_2, \end{align*} where \begin{align*} K_1&=\int_{R\leq|\tilde{z}_i|\leq R+1} \widetilde{Z}_{i0}\nabla Z_{i0}\nabla\eta_{1i}-Z_{i0}\nabla\widetilde{Z}_{i0} \nabla\eta_{1i} -\widetilde{Z}_{i0}\nabla\widehat{Z}_{i0}\nabla\eta_{1i} +\widehat{Z}_{i0}\nabla\widetilde{Z}_{i0}\nabla\eta_{1i}\\ &=\int_{R\leq|\tilde{z}_i|\leq R+1}\,\,Z_{i0}\nabla(Z_{i0} -\widehat{Z}_{i0})\nabla\eta_{1i} +(\widehat{Z}_{i0}-Z_{i0})\nabla(Z_{i0}-\widehat{Z}_{i0})\nabla\eta_{1i}\\ &\quad -(Z_{i0}-\widehat{Z}_{i0})\nabla\widehat{Z}_{i0}\nabla\eta_{1i}-(Z_{i0} -\widehat{Z}_{i0})^2|\nabla\eta_{1i}|^2, \end{align*} and \[ K_2= \int_{\frac{\kappa\delta}{4v_i\rho_i}\leq|\tilde{z}_i| \leq \frac{\kappa\delta}{3v_i\rho_i}} (\widetilde{Z}_{i0}\nabla\widehat{Z}_{i0} -\widehat{Z}_{i0}\nabla\widetilde{Z}_{i0})\nabla\eta_{2i} =-\int_{\Omega_{\varepsilon}} |\widehat{Z}_{i0}|^2|\nabla\eta_{2i}|^2=O(\frac1{p^2}). \] Note that for $R\leq|\tilde{z}_i|\leq R+1$, by \eqref{3.53a}-\eqref{3.53a1}, $|\nabla\widehat{Z}_{i0}|=\frac{\varepsilon_p}{v_i\rho_i} \{O(\frac1{R^{2(1+\alpha_i)+1}})+O(\frac1{pR})\}$. By \eqref{3.15j} and \eqref{3.15g}, it follows that \begin{gather*} \int_{R\leq|\tilde{z}_i|\leq R+1}(Z_{i0} -\widehat{Z}_{i0})^2|\nabla\eta_{1i}|^2dy=O(\frac1{p^2}), \\ \int_{R\leq|\tilde{z}_i|\leq R+1}(\widehat{Z}_{i0}-Z_{i0}) \nabla(Z_{i0}-\widehat{Z}_{i0})\nabla\eta_{1i}dy=O(\frac1{p^2}), \\ \int_{R\leq|\tilde{z}_i|\leq R+1}(Z_{i0}-\widehat{Z}_{i0}) \nabla\widehat{Z}_{i0}\nabla\eta_{1i}dy =O(\frac1p\cdot\frac1{R^{2(1+\alpha_i)+1}}). \end{gather*} Thus $$ K=K_1+K_2=\int_{R\leq|\tilde{z}_i|\leq R+1}Z_{i0} \nabla(Z_{i0}-\widehat{Z}_{i0})\nabla\eta_{1i}dy +O(\frac1p\cdot\frac1{R^{2(1+\alpha_i)+1}}), $$ which, together with \eqref{3.53a}, \eqref{3.53a1} and \eqref{3.15g}, implies that \begin{equation}\label{3.16j1} K=\begin{cases} -\frac{2\pi}{\mu_i^2|\log(v_i\rho_iR)|}\{1+O(\frac{1}{R^{2(1+\alpha_i)}})\} &\forall i\in J_1\cup J_3,\\ -\frac{\pi}{\mu_i^2|\log(v_i\rho_iR)|}\{1+O(\frac{1}{R^{2(1+\alpha_i)}})\} & \forall i\in J_2\cup J_4. \end{cases} \end{equation} As a consequence, estimate \eqref{3.16} can be derived from \eqref{3.16i}-\eqref{3.16j1}. \end{proof} \begin{proof}[Proof of Proposition \ref{prop3.1}] Let us prove \eqref{3.5} by contradiction. Assume that as $p_l\to + \infty$, \begin{equation}\label{3.5j} \|\phi_l\|_{\infty}=1,\quad p_l\|h_l\|_*\to0,\quad p_l\Big(\sum_{i\in J_3\cup J_4}|c^l_{i1}| +\sum_{i\in J_3}|c^l_{i2}|\Big)\geq\kappa>0. \end{equation} For convenience, we omit the dependence on $l$. By \eqref{3.11}, \begin{equation}\label{3.5f} \|\phi\|_{\infty}\leq Cp\big\{\|h\|_{*}+ \sum_{i\in J_3\cup J_4}|c_{i1}| +\sum_{i\in J_3}|c_{i2}|\big\}. \end{equation} Testing equation \eqref{3.3} against $\eta_{2i}Z_{ij}$ with $i\in J_3$, $j=1,2$, or $i\in J_4$, $j=1$, and integrating by parts, it implies that \begin{equation}\label{3.5b} |c_{ij}|\leq C\big\{|\langle \phi,\,L(\eta_{2i}Z_{ij})\rangle|+|\langle h,\,\eta_{2i}Z_{ij}\rangle|\big\}. \end{equation} By \eqref{3.53c}, it follows that for $|\tilde{z}_i|\leq\frac13\kappa\delta(v_i\rho_i)^{-1}$, \begin{align*} L(\eta_{2i}Z_{ij}) &=\eta_2(v_i\rho_i\tilde{z}_i)z_{ij}(\tilde{z}_i) \big\{W_q-(\frac{\varepsilon_p}{v_i\rho_i})^2\frac{8(1+\alpha_i)^2 |\tilde{z}_i|^{2\alpha_i}}{[1+|\tilde{z}_i|^{2(1+\alpha_i)}]^2}\big\}\\ &\quad +\big(\frac{\varepsilon_p}{v_i\rho_i}\big)^2 \big\{O(\rho_i|\tilde{z}_i|)\nabla_{\tilde{z}_i}^2z_{ij} +O(\rho_i)\nabla_{\tilde{z}_i}z_{ij}+O(\rho^2_i)z_{ij}\big\}(\tilde{z}_i). \end{align*} Note that $|z_{ij}|=O((1+|\tilde{z}|)^{-1-\alpha_i})$, $|\nabla_{\tilde{z}_i} z_{ij}|=O((1+|\tilde{z}|)^{-2-\alpha_i})$ and $|\nabla^2_{\tilde{z}_i}z_{ij}|=O((1+|\tilde{z}|)^{-3-\alpha_i})$. This combined \eqref{2.19}-\eqref{2.18c} and \eqref{3.53a}-\eqref{3.53a1} implies that for $|\tilde{z}_i|\leq\frac12\delta(v_i\rho_i)^{-1/2}$, \begin{align*} &L(\eta_{2i}Z_{ij})\\ &=\big(\frac{\varepsilon_p}{v_i\rho_i}\big)^2| \tilde{z}_i|^{2\alpha_i}z_{ij}(\tilde{z}_i)e^{U^i(\tilde{z}_i)} \big\{\frac{[\omega^{0i}-U^i-\frac12(U^i)^2](|\tilde{z}_i|)}p +O(\frac{\log^4(2+|\tilde{z}_i|)}{p^2})\big\}\\ &\quad +\,(\frac{\varepsilon_p}{v_i\rho_i})^2 \big\{O(\frac{|\tilde{z}_i|^{2\alpha_i}}{(1+|\tilde{z}_i|)^{5(1+\alpha_i)}}) +O(\frac{\rho_i}{(1+|\tilde{z}_i|)^{2+\alpha_i}}) +O(\frac{\rho^2_i}{(1+|\tilde{z}_i|)^{1+\alpha_i}})\big\}, \end{align*} and for $\frac12\delta(v_i\rho_i)^{-1/2}\leq|\tilde{z}_i| \leq\frac13\kappa\delta(v_i\rho_i)^{-1}$, $$ L(\eta_{2i}Z_{ij}) =(\frac{\varepsilon_p}{v_i\rho_i})^2O(\frac{1}{(1+|\tilde{z}_i|)^{5+3\alpha_i}} +\frac{\rho_i}{(1+|\tilde{z}_i|)^{2+\alpha_i}} +\frac{\rho^2_i}{(1+|\tilde{z}_i|)^{1+\alpha_i}}). $$ This implies that for any $i\in J_3$, $j=1,2$, or $i\in J_4$, $j=1$, \begin{equation}\label{3.5d} \langle \phi, L(\eta_{2i}Z_{ij})\rangle =\frac1p\,E_{ij}(\hat{\phi}_i)+O(\frac{\|\phi\|_{\infty}}{p^2}), \end{equation} where \[ \hat{\phi}_i(\tilde{z}_i)=\begin{cases} \phi(\frac{v_i\rho_iz_i+q_i}{\varepsilon_p}) & \forall i\in J_3,\\ \phi(\tilde{F}_i^\varepsilon(\tilde{z}_i)) & \forall i\in J_4, \end{cases} \] and for $i\in J_3$, \begin{gather*} E_{i1}(\hat{\phi}_i)=\int_{B(0,\frac12\delta(v_i\rho_i)^{-1/2})} \frac{32\operatorname{Re}(z_i)[\omega^{0i}-U^i-\frac{1}2(U^i)^2]} {[1+|z_i|^2]^3}\hat{\phi}_i dz_i,\\ E_{i2}(\hat{\phi}_i)=\int_{B(0,\frac12\delta(v_i\rho_i)^{-1/2})} \frac{32\operatorname{Im}(z_i)[\omega^{0i}-U^i-\frac{1}2(U^i)^2]} {[1+|z_i|^2]^3}\hat{\phi}_i dz_i, \end{gather*} and for $i\in J_4$, by \eqref{3.53a}-\eqref{3.53a1}, \begin{align*} E_{i1}(\hat{\phi}_i) &=\int_{\mathbb{R}_+^2\cap B(0,\frac12\delta(v_i\rho_i)^{-1/2})} \frac{32\operatorname{Re}(\tilde{z}_i)[\omega^{0i}-U^i-\frac{1}2(U^i)^2]} {[1+|\tilde{z}_i|^2]^3}\hat{\phi}_i\det[D\tilde{F}_i(\tilde{z}_i)]d\tilde{z}_i\\ &=\int_{\mathbb{R}_+^2\cap B(0,\frac12\delta(v_i\rho_i)^{-1/2})} \frac{32\operatorname{Re}(\tilde{z}_i)[\omega^{0i}-U^i-\frac{1}2(U^i)^2]} {[1+|\tilde{z}_i|^2]^3}\hat{\phi}_id\tilde{z}_i+O(\varepsilon_p\|\phi\|_{\infty}). \end{align*} Also \begin{equation}\label{3.5c} \langle h,\,\eta_{2i}Z_{ij}\rangle=O(\|h\|_*)=o(\frac1p). \end{equation} Substituting \eqref{3.5d}-\eqref{3.5c} in \eqref{3.5b}, it follows that \begin{equation}\label{3.5e} |c_{ij}|=O(\|h\|_*)+O(\frac{|\phi\|_\infty}p). \end{equation} Furthermore, $|c_{ij}|=O(\frac1p)$. Thus, similar to the blowup analysis in the proof of Lemma \ref{lem3.3}, it follows that there exists the constant $C_i$, $i\in J_3\cup J_4$, such that $$ \hat{\phi}_i(\tilde{z}_i)\to C_i\frac{|\tilde{z}_i|^2-1}{|\tilde{z}_i|^2+1} \quad\text{uniformly in } C^0_{loc}(\mathbb{R}^2). $$ From Lebesgue's theorem and the radial properties of $\omega^{0i}$ and $U^i$, it follows that for $i\in J_3$, \begin{gather*} E_{i1}(\hat{\phi}_i)\to C_i\int_{\mathbb{R}^2} \frac{32(|z_i|^2-1)\operatorname{Re}(z_i)} {[1+|z_i|^2]^4}[\omega^{0i}-U^i-\frac{(U^i)^2}2]dz_i=0,\\ E_{i2}(\hat{\phi}_i)\to C_i\int_{\mathbb{R}^2} \frac{32(|z_i|^2-1)\operatorname{Im}(z_i)} {[1+|z_i|^2]^4}[\omega^{0i}-U^i-\frac{(U^i)^2}2]dz_i=0, \end{gather*} and for $i\in J_4$, $$ E_{i1}(\hat{\phi}_i)\to C_i\int_{\mathbb{R}_+^2} \frac{32(|\tilde{z}_i|^2-1)\operatorname{Re}(\tilde{z}_i)} {[1+|\tilde{z}_i|^2]^4}[\omega^{0i}-U^i-\frac{(U^i)^2}2]d\tilde{z}_i=0, $$ which together with $\eqref{3.5d}$ implies $\langle \phi,\,L(\eta_{2i}Z_{ij})\rangle=o(\frac1p)$. Then by \eqref{3.5b}-\eqref{3.5c}, it implies that $p(\sum_{i\in J_3\cup J_4}|c_{i1}| +\sum_{i\in J_3}|c_{i2}|)=o(1)$, which contradicts with the assumption \eqref{3.5j}, and so estimate \eqref{3.5} is established. Moreover, by \eqref{3.5e}, it also implies that \begin{equation}\label{3.5g} |c_{ij}|\leq C\|h\|_{*}, \end{equation} which implies that there exists a unique trivial solution to problem \eqref{3.3} with $h\equiv 0$. Thus, from Fredholm's alternative, there exists a unique solution $\phi$, scalars $c_{ij}$, $i\in J_3$, $j=1,2$, and $c_{i1}$, $i\in J_4$, of problem \eqref{3.3} for any $h\in \mathcal {C}_{*}$, which completes the proof. \end{proof} \begin{remark} \label{rmk3.5}\rm Given $h\in L^{\infty}(\Omega_\varepsilon)$ with $\|h\|_*<\infty$, let $\phi$ be the solution to \eqref{3.3} given by Proposition \ref{prop3.1}. Multiplying the first equation in \eqref{3.3} by $\phi$ and integrating by parts, we obtain $$ \|\phi\|_{H^1}^2:=\int_{\Omega_\varepsilon}|\nabla\phi|^2 +\varepsilon_p^2\phi^2=\int_{\Omega_\varepsilon}W_q\phi^2 +\int_{\Omega_\varepsilon}h\phi. $$ Moreover, using Lemma \ref{lem2.3}, we can prove that $$ \big|\int_{\Omega_\varepsilon}W_q\phi^2\big|\leq C\|\phi\|_{\infty}^2, $$ and therefore $$ \|\phi\|_{H^1}\leq C(\|h\|_*+\|\phi\|_{\infty}). $$ \end{remark} Now, we solve the intermediate nonlinear problem: for any points $q\in\Lambda^{\tilde{m}}_{m}(\delta)$, we find a function $\phi$, scalars $c_{ij}$, $i\in J_3$, $j=1,2$, and $c_{i1}$, $i\in J_4$, such that \begin{equation}\label{4.1} \begin{gathered} L(\phi)=R_q+N(\phi)+\sum_{i\in J_3\cup J_4}c_{i1}\chi_iZ_{i1} +\sum_{i\in J_3}c_{i2}\chi_iZ_{i2} \quad \text{in } \Omega_{\varepsilon},\\ \frac{\partial}{\partial\nu}\phi=0\quad \text{on } \partial\Omega_{\varepsilon},\\ \int_{\Omega_{\varepsilon}}\chi_iZ_{ij}\phi=0 \quad \forall i\in J_3,\;j=1,2;\; i\in J_4,\;j=1. \end{gathered} \end{equation} \begin{proposition} \label{prop3.6} There exist $C>0$ and $p_0>1$ such that for any $p>p_0$ and $q\in\Lambda^{\tilde{m}}_{m}(\delta)$, problem \eqref{4.1} admits a unique solution $\phi\in L^{\infty}(\Omega_{\varepsilon})$, scalars $c_{ij}$, $i\in J_3$, $j=1,2$, and $c_{i1}$, $i\in J_4$, such that \begin{equation}\label{4.2} \|\phi\|_{\infty}\leq C/p^3,\quad \|\phi\|_{H^1}\leq C/p^3. \quad |c_{ij}|\leq C/p^4. \end{equation} \end{proposition} The proof of the above proposition can be done along the lines of those of \cite[Lemma 4.1]{MW1}; we omit it here. \begin{remark} \label{rmk3.7}\rm Using the fixed point characterization of the solution $\phi=\phi(q)$ to \eqref{4.1}, the Implicit Function Theorem and Remark \ref{rmk3.5}, we can easily verify that $\phi(q)$ is differentiable with respect to $q\in\Lambda^{\tilde{m}}_{m}(\delta)$, in $L^{\infty}(\Omega_\varepsilon)$ and $H^1(\Omega_\varepsilon)$. \end{remark} \begin{remark} \label{rmk3.8}\rm The function $V_q+\phi$, with $\phi$ given by Proposition \ref{prop3.6}, is positive in $\overline{\Omega}_\varepsilon$. In fact, we observe first that $p|\phi|\to0$ uniformly over compacts of $\overline{\Omega}_\varepsilon$. In addition, from \eqref{2.7c} we argue that, in the region close to some point $q_i$, $V_q+\phi$ is positive. Outside this region, we may conclude the same from \eqref{2.6a} and \eqref{2.13a}. \end{remark} \section{Variational reduction} After problem \eqref{4.1} has been solved, we find a solution of \eqref{2.23} with $m+\tilde{m}\geq 1$, and hence for \eqref{a} if $q\in\Lambda^{\tilde{m}}_{m}(\delta)$ satisfies \begin{equation}\label{5.1} c_{ij}(q)=0\quad \forall i\in J_3,\; j=1,2;\; i\in J_4,\; j=1. \end{equation} To solve it we consider the energy functional of \eqref{a}, \begin{equation}\label{6.1} J_p(u)=\frac1{2}\int_{\Omega}(|\nabla u|^2+u^2)dx-\frac1{p+1}\int_{\Omega}S(x)u^{p+1}dx, \end{equation} and its finite-dimensional restriction \begin{equation}\label{6.1aa} F_p(q)=J_p(U_q+\tilde{\phi})\quad \forall q\in\Lambda^{\tilde{m}}_{m}(\delta), \end{equation} where $U_q$ is defined in \eqref{2.6} and \begin{equation}\label{6.1a} \tilde{\phi}(q)(x)=\varepsilon_p^{-\frac{2}{p-1}}\phi(q)(\varepsilon_p^{-1}x). \end{equation} The following proposition tells us that critical points of $F_p$ correspond to solutions of \eqref{5.1}. \begin{proposition} \label{prop4.1} The functional $F_p$ is of class $C^1$. Moreover, for $p$ large enough, if $D_{q}F_p(q)=0$, then $q$ satisfies \eqref{5.1}. \end{proposition} \begin{proof} A direct consequence of the results in Remark \ref{rmk3.7} is the fact that $F_p(q)$ is a $C^1$-function of $q$ since the map $q\to\phi$ is a $C^1$-map in $H^1(\Omega_\varepsilon)$. Define \begin{equation}\label{6.2} I_p(v)=\frac1{2}\int_{\Omega_{\varepsilon}}(|\nabla v|^2+\varepsilon_p^2v^2)dy-\frac1{p+1}\int_{\Omega_{\varepsilon}} S(\varepsilon_p y)v^{p+1}dy. \end{equation} Then, making a change of variable, we have \begin{equation}\label{6.2a} \varepsilon_p^{\frac4{p-1}}F_p(q)=I_p(V_q+\phi). \end{equation} Since $D_qF_p(q)=0$, we have that \begin{equation} \label{5.5} \begin{aligned} 0&=D I_p(V_q+\phi)[D_q V_q+D_q\phi] \\ &=\sum_{i,j}c_{ij}\int_{\Omega_{\varepsilon}}\chi_i Z_{ij} (D_q V_q+D_q\phi) \\ &=\sum_{i,j}c_{ij}\int_{\Omega_{\varepsilon}}\chi_i Z_{ij} D_q V_q-\phi D_q[\chi_i Z_{ij}], \end{aligned} \end{equation} because $\int_{\Omega_{\varepsilon}}\chi_iZ_{ij}\phi=0$. By the expression of $V_q$ in \eqref{2.13a}, a direct computation shows that for any $l\in J_3$, $s=1,2$, and $l\in J_4$, $s=1$, $$ \partial_{(q_l)_s}V_q=\frac{\frac{1}{\varepsilon_p\mu_l}}{p^{\frac{p}{p-1}}\mu_l^{\frac2{p-1}} c_l(q_l)^{\frac1{p-1}}}\big\{\frac{4(z_l)_s} {1+|z_l|^2}-\frac{\partial_{(z_l)_s}\omega^{0l}}{p} -\frac{\partial_{(z_l)_s}\omega^{1l}}{p^2}\big\}+O(\frac{\varepsilon^\tau_p}p). $$ Consequently, \eqref{5.5} can be written as, for each $l\in J_3$, $s=1,2$, and $l\in J_4$, $s=1$, \begin{equation} \Big[\frac{\mu_lA_l}{p^{\frac p{p-1}}\mu_l^{\frac2{p-1}} c_l(q_l)^{\frac1{p-1}}}\int_0^{R_0+1}\chi(r)\frac{16r^{3}}{(\,1+r^2)^2} dr\Big]c_{ls}+\sum_{i,j}c_{ij}O(\frac1{p^2})=0, \end{equation} where $A_l=\pi$ for $l\in J_3$, and $A_l=\frac12\pi$ for $l\in J_4$. This is a strictly diagonal dominant system for $p$ sufficiently large. We thus get that $c_{ls}=0$ for any $l\in J_3$, $s=1,2$, or $l\in J_4$, $s=1$. \end{proof} Next, we need to give the expansion of $F_p$ in terms of $\varphi^{\tilde{m}}_{m}$ defined in \eqref{b} as $p$ goes to $+\infty$. \begin{proposition} \label{prop4.2} There exist constants $k_1$, $k_2$ and $k_3>0$ depending only on the points $q_i$, $i\in J_1\cup J_2$, such that \begin{equation}\label{5.3} F_p(q)=\frac{k_1}{p}-\frac{2k_1}{p^2}\log p +\frac{k_2}{p^2}-\frac{k_3}{p^2}\varphi^{\tilde{m}}_{m}(q)+O(\frac{\log^2p}{p^3}), \end{equation} uniformly for all points $q\in\Lambda_{m}^{\tilde{m}}(\delta)$. \end{proposition} \begin{proof} Multiplying the first equation in \eqref{4.1} by $V_q+\phi$ and integrating by parts, we obtain $$ \int_{\Omega_{\varepsilon}}|\nabla (V_q+\phi)|^2+\varepsilon_p^2(V_q+\phi)^2 =\int_{\Omega_{\varepsilon}}S(\varepsilon_p y)(V_q+\phi)^{p+1} +\sum_{i,j}c_{ij}\int_{\Omega_\varepsilon}\chi_iZ_{ij}V_q. $$ Since $V_q$ is a bounded function, by \eqref{4.2} we obtain that $$ \int_{\Omega_{\varepsilon}}|\nabla (V_q+\phi)|^2+\varepsilon_p^2(V_q+\phi)^2 =\int_{\Omega_{\varepsilon}}S(\varepsilon_p y)(V_q+\phi)^{p+1}+O(\frac1{p^4}), $$ uniformly for $q\in\Lambda_{m}^{\tilde{m}}(\delta)$. Hence, by \eqref{6.2} and \eqref{6.2a} we have \begin{equation}\label{5.4} F_p(q)=(\frac1{2}-\frac{1}{p+1})\varepsilon_p^{-\frac4{p-1}} \int_{\Omega_\varepsilon}|\nabla (V_q+\phi)|^2+\varepsilon_p^2(V_q+\phi)^2+O(\frac1{p^4}). \end{equation} We expand the term $\int_{\Omega_\varepsilon}|\nabla V_q|^2+\varepsilon_p^2V_q^2$: in view of \eqref{2.7c} and \eqref{2.13c} we obtain \begin{align*} &\int_{\Omega_\varepsilon}|\nabla V_q|^2+\varepsilon_p^2V_q^2\\ &=\int_{\Omega_\varepsilon}V_q(-\Delta V_q+\varepsilon_p^2V_q)\\ &=\sum_{i\in\cup_{l=1}^4J_l}\frac{1}{p^{\frac{2p}{p-1}}\mu_i^{\frac4{p-1}} c_i(q_i)^{\frac2{p-1}}}\int_{|z_i|\leq\delta(v_i\rho_i)^{-1}}|z_i|^{2\alpha_i}\\ &\quad\times\Big\{e^{U^i}+\frac{e^{U^i}\omega^{0i}-f^{0i}}p +\frac{e^{U^i}\omega^{1i}-f^{1i}}{p^2}\Big\} \big\{p+U^i+\frac{1}p\omega^{0i}+\frac{1}{p^2}\omega^{1i}\big\}dz_i +O(\frac1{p^3})\\ &=\sum_{i\in\cup_{l=1}^4J_l}\frac{A_i}{p^{\frac{2p}{p-1}}\mu_i^{\frac4{p-1}} c_i(q_i)^{\frac2{p-1}}}\int_{\mathbb{R}^2}\Big\{p|z_i|^{2\alpha_i} e^{U^i}+|z_i|^{2\alpha_i}[U^i-\frac{(U^i)^2}{2}+\omega^{0i}]e^{U^i}\Big\}\\ &\quad +O(\frac1{p^3})\\ &=\sum_{i\in\cup_{l=1}^4J_l}\big\{\frac1p [1-\frac{2\log p}{p}-\frac{2}{p}\log c_i(q_i)]d_i -\frac4{p^2}d_i\log\mu_i+\frac1{p^2}A_ie_i\big\} +O(\frac{\log^2p}{p^3}), \end{align*} where $A_i=1$ for $i\in J_1\cup J_3$, $A_i=\frac12$ for $i\in J_2\cup J_4$, and the last equality is due to the following relations: \begin{gather*} p^{-\frac{2p}{p-1}}=\frac1{p^2}-\frac{2}{p^3}\log p+O(\frac{\log^2p}{p^4}),\\ \mu_i^{-\frac4{p-1}}=1-\frac4p\log\mu_i+O(\frac1{p^2}),\\ c_i(q_i)^{-\frac{2}{p-1}}=1-\frac{2}{p}\log c_i(q_i)+O(\frac1{p^2}),\\ e_i=\int_{\mathbb{R}^2}|z_i|^{2\alpha_i}[U^i-\frac12(U^i)^2+\omega^{0i}]e^{U^i} \quad \forall i\in\cup_{l=1}^4 J_l. \end{gather*} From the expansion of $\mu_i$ in \eqref{2.6bb}, we have \begin{equation} \label{5.2} \begin{aligned} &\int_{\Omega_\varepsilon}|\nabla V_q|^2+\varepsilon_p^2V_q^2 \\ &=\sum_{i\in\cup_{l=1}^{4}J_l}\Big\{(1-\frac{2\log p}p +\frac3p)\frac{d_i}p+\frac{A_ie_i}{p^2}-\frac{d_i}{p^2}\ \Big[2\log c_i(q_i)\\ &\quad +d_i H(q_i,q_i)+\sum_{j\neq i}d_j G(q_i,q_j)\Big]\Big\} +O(\frac{\log^2p}{p^3})\\ &=-\frac{1}{p^2}\varphi^{\tilde{m}}_{m}(q) +\sum_{i\in\cup_{l=1}^{4}J_l}\big\{(1-\frac{2\log p}p+\frac3p) \frac{d_i}p+\frac{A_ie_i}{p^2}\big\} \\ &\quad -\sum_{i\in J_1\cup J_2}\frac{d_i}{p^2}\Big\{2\log c_i(q_i) +d_iH(q_i,q_i)+\sum_{j\in J_1\cup J_2,\,j\neq i}d_jG(q_i,q_j)\Big\} +O(\frac{\log^2p}{p^3}), \end{aligned} \end{equation} uniformly for $q\in\Lambda_{m}^{\tilde{m}}(\delta)$. In particular, \begin{equation}\label{5.2a} \|V_q\|_{H^1}^2=\int_{\Omega_\varepsilon}|\nabla V_q|^2+\varepsilon_p^2V_q^2=O(\frac1p). \end{equation} Furthermore, by \eqref{4.2} and \eqref{5.2a}, we obtain \begin{equation}\label{5.2aa} 2\int_{\Omega_\varepsilon}[\nabla V_q\nabla\phi+\varepsilon_p^2V_q\phi]+\int_{\Omega_\varepsilon}[|\nabla \phi|^2+\varepsilon_p^2\phi^2]=O(\frac1{p^{7/2}}). \end{equation} Also \begin{equation}\label{5.6} \varepsilon_p^{-\frac4{p-1}}=e^{\frac{p}{p-1}}=e+\frac{e}p+O(\frac1{p^2}). \end{equation} Thus, inserting \eqref{5.2}, \eqref{5.2aa} and \eqref{5.6} in \eqref{5.4}, we obtain that \eqref{5.3} holds. \end{proof} \begin{proof}[Proof of Theorem \ref{thm1.2}] First of all, from Proposition \ref{prop4.1}, we can provide a solution to problem \eqref{a} if we adjust $q\in\Lambda_{m}^{\tilde{m}}(\delta)$ so that it is a critical point of $F_p(q)$ defined by \eqref{6.1aa}. This is equivalent to finding a critical point of $$ \tilde{F}_p(q):=\frac1{k_3}\left(k_1p-2k_1\log p+k_2-p^2F_p(q)\right), $$ for suitable constants $k_1$, $k_2$ and $k_3$. On the other hand, from Proposition \ref{prop4.2}, we have \begin{equation}\label{r4} \tilde{F}_p(q)=\varphi^{\tilde{m}}_{m}(q)+O(\frac{\log^2 p}{p}), \end{equation} uniformly for all points $q\in\Lambda_{m}^{\tilde{m}}(\delta)$ as $p\to\infty$, where $\varphi^{\tilde{m}}_{m}$ is given by \eqref{b}. Next, as in \cite[Lemma 6.1]{MW1}, we have \begin{equation}\label{6.1b} \min_{q\in\partial\Lambda_{m}^{\tilde{m}}(\delta)}\varphi_{m}^{\tilde{m}}(q) \to +\infty \quad\text{as }\delta\to0. \end{equation} Thus, for $\delta$ small enough, $\varphi_{m}^{\tilde{m}}$ has a global minimum $M$ in $\Lambda_{m}^{\tilde{m}}(\delta)$. This, together with \eqref{r4}, implies that $\tilde{F}_p(q)$ has also a global minimum point $q^p\in\Lambda_{m}^{\tilde{m}}(\delta)$ such that $\varphi_{m}^{\tilde{m}}(q^p)\to M$ as $p\to\infty$. Moreover, up to a subsequence, there exists a global minimum point $\tilde{q}$ of $\varphi_{m}^{\tilde{m}}$ in $\Lambda^{\tilde{m}}_{m}(\delta)$ such that $q^p\to\tilde{q}$ as $p\to\infty$. The function $u_p=U_{q^p}+\tilde{\phi}(q^p)$, where $\tilde{\phi}(q^p)$ is defined in \eqref{6.1a}, is therefore a solution to \eqref{a} with the qualitative properties stated in the theorem. \end{proof} \subsection*{Acknowledgments} This research was supported by the National Natural Science Foundation of China (No. 11171214). \begin{thebibliography}{00} \bibitem{AG} Adimurthi, M. Grossi; \emph{Asymptotic estimates for a two-dimensional problem with polynomial nonlinearity}, Proc. Amer. Math. Soc., \textbf{132} (2004), 1013--1019. \bibitem{BP} S. Baraket, F. Parcard; \emph{Construction of singular limits for a semilinear elliptic equation in dimension 2}, Calc. Var. Partial Differential Equations, \textbf{6} (1998), 1--38. \bibitem{CL5} C.-C. Chen, C.-S. Lin; \emph{Sharp estimates for solutions of multi-bubbles in compact Riemann surfaces}, Comm. Pure Appl. Math., \textbf{55} (2002), 728--771. \bibitem{CY1} Y.-B. Chang, H.-T. Yang, \emph{Singular limits for 2-dimensional elliptic Neumann problem with exponential nonlinearity and singular data}, Commun. Contemp. Math., \textbf{15} (2013), 1350024 (43 pages). \bibitem{CY2} Y.-B. Chang, H.-T. Yang; \emph{Concentrating solutions for a two-dimensional elliptic problem with large exponent in weighted nonlinearity}, Complex Var. Elliptic Equ., \textbf{59} (2014), 1096--1117. \bibitem{D} T. D'Aprile; \emph{Multiple blow-up solutions for the Liouville equation with singular data}, Comm. Partial Differential Equations, \textbf{38} (2013), 1409--1436. \bibitem{DDM} J. D\'{a}vila, M. del Pino, M. Musso; \emph{Concentrating solutions in a two-dimensional elliptic problem with exponential Neumann data}, J. Funct. Anal., \textbf{227} (2005), 430--490. \bibitem{DKM} M. del Pino, M. Kowalczyk, M. Musso; \emph{Singular limits in Liouville-type equations}, Calc. Var. Partial Differential Equations, \textbf{24} (2005), 47--81. \bibitem{DFW} M. del Pino, P. L. Felmer, J.-C. Wei; \emph{On the role of mean curvature in some singularly perturbed Neumann problems}, SIAM J. Math. Anal., \textbf{31} (1999), 63--79. \bibitem{E} P. Esposito; \emph{Blowup solutions for a Liouville equation with singular data}, SIAM J. Math. Anal., \textbf{36} (2005), 1310--1345. \bibitem{EGP} P. Esposito, M. Grossi, A. Pistoia, \emph{On the existence of blowing-up solutions for a mean field equation}, Ann. Inst. H. Poincar\'e Analyse Non Lin\'eaire, \textbf{22} (2005), 227--257. \bibitem{EMP1} P. Esposito, M. Musso, A. Pistoia, \emph{Concentrating solutions for a planar elliptic problem involving nonlinearities with large exponent}, J. Differential Equations, \textbf{227} (2006), 29--68. \bibitem{EPW} P. Esposito, A. Pistoia, J.-C. Wei; \emph{Concentrationg solutions for the H\'enon equation in $\mathbb{R}^2$}, J. Anal. Math., \textbf{100} (2006), 249--280. \bibitem{EG} K. El Mehdi, M. Grossi; \emph{Asymptotic estimates and qualitative properties of an elliptic problem in dimension two}, Adv. Nonlinear Stud., \textbf{4} (2004), 15--36. \bibitem{GNN} B. Gidas, W.-M. Ni, L. Nirenberg; \emph{Symmetry and related properties via the maximum principle}, Comm. Math. Phys., \textbf{68} (1979), 209--243. \bibitem{GW} C.-F. Gui, J.-C. Wei; \emph{Multiple interior spike solutions for some singular perturbed Neumann problems}, J. Differential Equations, \textbf{158} (1999), 1--27. \bibitem{GW1} C.-F. Gui, J.-C. Wei; \emph{On multiple mixed interior and boundary peak solutions for some singularly perturbed Neumann problems}, Canad. J. Math., \textbf{52} (2000), 522--538. \bibitem{GWW} C.-F. Gui, J.-C. Wei, M. Winter; \emph{Multiple boundary peak solutions for some singularly perturbed Neumann problems}, Ann. Inst. H. Poincar\'e Anal. Non Lin\'eaire, \textbf{17} (2000), 249--289. \bibitem{H} M. H\'enon; \emph{Numerical experiments on the spherical stellar systems}, Astronomy and Astrophysics, \textbf{24} (1973), 229--238. \bibitem{KS} E. F. Keller, L. A. Segel, \emph{Initiation of slime mold aggregation viewed as an instability}, J. Theoret. Biol., \textbf{26} (1970), 399--415. \bibitem{LNT} C.-S. Lin, W.-M. Ni, I. Takagi; \emph{Large amplitude stationary solutions to a chemotaxis system}, J. Differential Equations, \textbf{72} (1988), 1--27. \bibitem{L} Y.-Y. Li; \emph{On a singularly perturbed equation with Neumann boundary condition}, Comm. Partial Differential Equations, \textbf{23} (1998), 487--545. \bibitem{MW1} M. Musso, J.-C. Wei; \emph{Stationary solutions to a Keller-Segel chemotaxis system}, Asymptot. Anal., \textbf{49} (2006), 217--247. \bibitem{N1} V. Nanjundiah; \emph{Chemotaxis, signal relaying and aggregation morphology}, J. Theoret. Biol., \textbf{42} (1973), 63--105. \bibitem{NT1} W.-M. Ni, I. Takagi; \emph{On the shape of least-energy solutions to a semilinear Neumann problem}, Comm. Pure Appl. Math., \textbf{44} (1991), 819--851. \bibitem{NT2} W.-M. Ni, I. Takagi; \emph{Locating the peaks of least-energy solutions to a semilinear Neumann problem}, Duke Math. J., \textbf{70} (1993), 247--281. \bibitem{PW} M. H. Protter, H. F. Weinberger; \emph{Maximum Principles in Differential Equations}, Corrected Reprint of the 1967 Original, Springer-Verlag, New York, 1984. \bibitem{RW} X.-F. Ren, W.-J. Wei; \emph{On a two-dimensional elliptic problem with large exponent in nonlinearity}, Trans. Amer. Math. Soc., \textbf{343} (1994), 749--763. \bibitem{RW1} X.-F. Ren, W.-J. Wei; \emph{Single-point condensation and least-energy solutions}, Proc. Amer. Math. Soc., \textbf{124} (1996), 111--120. \end{thebibliography} \end{document}