\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2016 (2016), No. 103, pp. 1--21.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2016 Texas State University.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2016/103\hfil Concentrating solutions] {Multi-peak solutions for a planar Robin nonlinear elliptic problem with large exponent} \author[Y. Zhang, L. Shi \hfil EJDE-2016/103\hfilneg] {Yibin Zhang, Lei Shi} \address{Yibin Zhang \newline College of Sciences, Nanjing Agricultural University, Nanjing 210095, China} \email{yibin10201029@njau.edu.cn} \address{Lei Shi \newline College of Sciences, Nanjing Agricultural University, Nanjing 210095, China} \email{shileijsxh@163.com} \thanks{Submitted November 26, 2015. Published April 21, 2016.} \subjclass[2010]{35J25, 35B25, 35B38} \keywords{Concentrating solutions; large exponent; Robin boundary condition; \hfill\break\indent finite-dimensional reduction} \begin{abstract} We consider the elliptic equation $\Delta u+u^p=0$ in a bounded smooth domain $\Omega$ in $\mathbb{R}^2$ subject to the Robin boundary condition $\frac{\partial u}{\partial\nu} +\lambda b(x)u=0$. Here $\nu$ denotes the unit outward normal vector on $\partial\Omega$, $b(x)$ is a smooth positive function defined on $\partial\Omega$, $0<\lambda<+\infty$, and $p$ is a large exponent. For any fixed $\lambda$ large we find topological conditions on $\Omega$ which ensure the existence of a positive solution with exactly $m$ peaks separated by a uniform positive distance from the boundary and each from other as $p\to+\infty$ and $\lambda\to+\infty$. In particular, for a nonsimply connected domain such solution exists for any $m\geq1$. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{remark}[theorem]{Remark} \newtheorem{definition}[theorem]{Definition} \allowdisplaybreaks \section{Introduction} In this article we consider the boundary-value problem \begin{equation}\label{a1} \begin{gathered} \Delta u+u^p=0 \quad \text{in }\Omega,\\ u>0 \quad \text{in } \Omega,\\ \frac{\partial u}{\partial\nu} +\lambda b(x)u=0 \quad \text{on } \partial\Omega, \end{gathered} \end{equation} where $\Omega$ is a bounded smooth domain in $\mathbb{R}^2$, $\nu$ denotes the unit outward normal vector on $\partial\Omega$, $b(x)$ is a smooth positive function defined on $\partial\Omega$, $0<\lambda<+\infty$, and $p$ is a large exponent. The boundary condition in problem \eqref{a1} is called Robin boundary condition. Such an Robin boundary condition is particularly interesting in various branches of biological models (see \cite{DMO,N}). When $\lambda=0$, from integration by parts it is trivial to observe that \eqref{a1} has no solution. On the other hand, if $0<\lambda\leq+\infty$, it is easy to prove via standard variational methods that \eqref{a1} always has a least energy solution. Moreover, in the case $\lambda=+\infty$, problem \eqref{a1} is reduced to the problem \begin{equation}\label{a2} \begin{gathered} \Delta u+u^p=0\quad \text{in } \Omega,\\ u>0\quad\text{in } \Omega,\\ u=0\quad \text{on }\partial\Omega. \end{gathered} \end{equation} Ren and Wei \cite{RW1,RW2} showed that the least energy solution $u_p$ of \eqref{a2} develops one interior peak, namely $u_p$ approaches zero except one interior point where it has an $L^\infty$-norm bounded and bounded away from zero, uniformly in $p$ as $p\to+\infty$. More precisely, the authors prove that, up to a subsequence, both $p|\nabla u_p|^2$ and $pu_p^{p+1}$ behave as a Dirac mass near a critical point of the Robin function $H_\infty(x,x)$, where $H_\infty$ is the regular part of Green's function $G_\infty$ of the Dirichlet Laplacian in $\Omega$, i.e. $H_\infty(x,y)=G_\infty(x,y)+\frac1{2\pi}\log|x-y|$. Successively, in \cite{AG,EG} the authors give a further description of the asymptotic behavior of $u_p$, as $p\to+\infty$, by identifying a limit profile problem of Liouville-type $\Delta u+e^u=0$ in $\mathbb{R}^2$, $\int_{\mathbb{R}^2}e^u<+\infty$, and showing that $\|u_p\|_{\infty}\to\sqrt{e}$ as $p\to+\infty$. Furthermore, Esposito, Musso and Pistoia \cite{EMP} prove that if $\Omega$ is not simply connected, \eqref{a2} can have many other positive solutions which, as $p$ tends to infinity, concentrate at $m$ different points of $\Omega$, i.e. \begin{equation}\label{a21} pu_p^{p+1}\rightharpoonup8\pi e\sum_{i=1}^m\delta_{\xi_i} \quad \text{weakly in the sense of measure in } \overline{\Omega}, \end{equation} where points $\xi=(\xi_1,\ldots,\xi_m)\in\Omega^m$ corresponds to a critical point of the function \begin{equation}\label{a24} \varphi_m^\infty\big(\xi_1,\ldots,\xi_m\big) =\sum_{j=1}^mH_\infty(\xi_j,\xi_j)+\sum_{j\neq k}G_\infty(\xi_j,\xi_k). \end{equation} In contrast, Grossi and Takahashi \cite{GT1} prove that when $\Omega$ is convex, problem \eqref{a2} has no multi-peak solutions satisfying \eqref{a21}. Thus the assumption on the domain in \cite{EMP} is sharp for the construction of multiple concentrating solutions of \eqref{a2}. The purpose of our research is to give the construction of multi-peak solutions to the so called Robin problem \eqref{a1} with sufficiently large $p$ and $\lambda$, and to point out that in general the set of multi-peak solutions of this problem exhibits a richer structure than the problem with Dirichlet boundary condition, which we will finish in this paper and in \cite{ZS}. In this paper we prove that if $\Omega$ is not simple connected, then given any $m\geq1$, for $p$ and $\lambda$ large enough problem \eqref{a1} has a positive solution $u_{p,\lambda}$ concentrating at exactly $m$ points that stay uniformly separated from the boundary and from each other as $p\to+\infty$ and $\lambda\to+\infty$. In particular, we recover existence results already known in \cite{EMP} when $\lambda=+\infty$ and $p$ is large enough. To state our results, we need to introduce some notation. Let $G_\lambda(x,y)$ be the Green's function satisfying \begin{equation}\label{a5} \begin{gathered} -\Delta_x G_\lambda(x,y)=\delta_y(x)\quad x\in\Omega,\\ \frac{\partial G_\lambda}{\partial\nu}(x,y)+ \lambda b(x)G_\lambda(x,y)=0\quad x\in\partial\Omega, \end{gathered} \end{equation} then its regular part can be decomposed as \begin{equation}\label{a6} H_\lambda(x,y)=G_\lambda(x,y)-\frac1{2\pi}\log\frac1{|x-y|}. \end{equation} Furthermore, let \begin{equation}\label{5.8} \varphi_m^\lambda\big(\xi_1,\ldots,\xi_m\big) =\sum_{j=1}^mH_\lambda(\xi_j,\xi_j)+\sum_{j\neq k}G_\lambda(\xi_j,\xi_k). \end{equation} Our main result reads as follows. \begin{theorem} \label{thm1.1} Assume that $\Omega$ is not simply connected. Then given any $m\geq1$, there exist $p_m>0$ and $\lambda_m>0$ such that for any $p>p_m$ and $\lambda>\lambda_m$, problem \eqref{a1} has a solution $u_{p,\lambda}$ with $m$ concentration points $\xi_{1,p,\lambda},\ldots,\xi_{m,p,\lambda}$ separated at a uniform positive distance from the boundary and each other as $p\to+\infty$ and $\lambda\to+\infty$. More precisely, \[ u_{p,\lambda}(x)=\sum_{j=1}^{m}\frac{1}{\gamma\mu_j^{2/(p-1)}} \Big[\log \frac1{(\delta_j^2+|x-\xi_{j,p,\lambda}|^2)^2} +8\pi H_\lambda(x,\xi_{j,p,\lambda})\Big]+O(\frac1p), \] where the parameters $\gamma$, $\delta_j$ and $\mu_j$ satisfy \[ \gamma=p^{\frac{p}{p-1}}\rho^{\frac2{p-1}},\quad \delta_j=\mu_j\rho, \quad \rho=e^{-\frac14p}, \quad \frac1{C}<\mu_j0$, and $\xi_{p,\lambda}=(\xi_{1,p,\lambda},\ldots, \xi_{m,p,\lambda})\in\Omega^m$ satisfies $$ \lim_{p\to+\infty,\,\,\lambda\to+\infty} \nabla\varphi_m^\lambda\big(\xi_{1,p,\lambda},\ldots,\xi_{m,p,\lambda}\big)=0, $$ and $$ \operatorname{dist}(\xi_{j,p,\lambda},\partial\Omega)\geq2\varepsilon, \quad |\xi_{j,p,\lambda}-\xi_{k,p,\lambda}|\geq2\varepsilon \quad \forall j,k=1,\ldots,m;\; j\neq k, $$ for any $\varepsilon>0$ small. In particular, as $p\to+\infty$ and $\lambda\to+\infty$, \begin{gather*} pu_{p,\lambda}^{p+1}-8\pi e\sum_{j=1}^m\delta_{\xi_{j,p,\lambda}}\rightharpoonup0 \quad\text{weakly in the sense of measure in } \overline{\Omega}, \\ u_{p,\lambda}\to 0 \quad\text{uniformly in } \overline{\Omega}\setminus\cup_{j=1}^mB_{\varepsilon}(\xi_{j,p,\lambda}), \\ \sup_{B_{\varepsilon}(\xi_{j,p,\lambda})}u_{p,\lambda}\to\sqrt{e}. \end{gather*} \end{theorem} The rest of this article is devoted to the proof of Theorem \ref{thm1.1}. Our proof relies on a Lyapunov-Schmidt process as in \cite{DKM1,EGP,EMP,MW1}, but we now have to confront some difficulties that are brought by the presence of Robin boundary condition, which can be successfully overcome by making use of some versions of the maximum principle with Robin boundary condition. This is the delicate ingredient during we construct multi-peak solutions of problem \eqref{a1} through performing the finite-dimensional reduction and using the notion of a nontrivial critical level. This article is organized as follows. In Section $2$ we exactly describe the ansatz for the solution of problem \eqref{a1} and estimate the error. Then we rewrite problem \eqref{a1} in terms of a linearized operator for which a solvability theory, subject to suitable orthogonality conditions, is performed through solving a linearized problem in Section $3$. In Section $4$ we solve an auxiliary nonlinear problem. In Section $5$ we reduce \eqref{a1} to a finite system, as we will see in Section $5$. In the last section, we use the notion of a nontrivial critical level to give the proof of Theorem \ref{thm1.1}. \section{A first approximation of the solution} In this section we provide an ansatz for solutions of problem \eqref{a1}. A key ingredient to describe an approximate solution of problem \eqref{a1} is given by the standard bubble: \begin{equation}\label{2.1} U_{\delta,\xi}(x)=\log\frac{8\delta^2}{(\delta^2+|x-\xi|^2)^2},\quad \delta>0,\quad \xi\in\mathbb{R}^2. \end{equation} It is well known (see \cite{CL}) that those are all the solutions of the Liouville-type equation \begin{equation}\label{2.2} \begin{gathered} \Delta u+e^u=0\quad \text{in } \mathbb{R}^2,\\ \int_{\mathbb{R}^2}e^u<+\infty. \end{gathered} \end{equation} Let us introduce the configuration space in which the concentration points belong to \begin{equation}\label{2.3} \begin{aligned} \mathcal{O}_\varepsilon:=\Big\{&\xi=(\xi_1,\ldots,\xi_m)\in\Omega^m: \operatorname{dist}(\xi_j,\partial\Omega)\geq2\varepsilon, \; |\xi_j-\xi_k|\geq2\varepsilon,\\ & j,k=1,\ldots,m;\; j\neq k\Big\}, \end{aligned} \end{equation} where $\varepsilon>0$ is a sufficiently small but fixed number. Furthermore, we set, for each $j=1,\ldots,m$, \begin{equation}\label{2.4} \gamma=p^{\frac{p}{p-1}}\rho^{\frac2{p-1}},\quad \delta_j=\mu_j\rho,\quad \rho=e^{-\frac14p},\quad \frac1C<\mu_j0$, where the choice of $\mu_j$ will be determined later. Define now \begin{equation}\label{2.5} U_j(x)=\frac{1}{\gamma\mu_j^{2/(p-1)}}\Big[U_{\delta_j,\xi_j}(x) +\frac1p\omega_1\Big(\frac{x-\xi_j}{\delta_j}\Big) +\frac1{p^2}\omega_2\Big(\frac{x-\xi_j}{\delta_j}\Big) \Big]. \end{equation} Here, $\omega_1$ and $\omega_2$ are radial solutions of \begin{equation}\label{2.6} \Delta\omega_i+\frac{8}{(1+|y|^2)^2}\omega_i=\frac{8}{(1+|y|^2)^2}f_i(y) \quad\text{in }\mathbb{R}^2, \end{equation} for $i=1,2$, respectively, with \begin{equation}\label{2.7} f_1=\frac12U_{1,0}^2,\quad f_2=\omega_1U_{1,0}-\frac13U_{1,0}^3 -\frac12\omega_1^2-\frac18U_{1,0}^4+\frac{1}{2}\omega_1U_{1,0}^2, \end{equation} having asymptotic properties \begin{equation}\label{2.8} \begin{gathered} \omega_i(y)=C_i\log |y|+O\big(\frac1{|y|}\big) \quad \text{as } |y|\to+\infty, \\ \nabla\omega_i(y)=C_i\frac{y}{1+|y|^2}+O\Big(\frac{1}{1+|y|^2}\Big) \quad \text{for all }y\in\mathbb{R}^2, \end{gathered} \end{equation} for $i=1,2$, where \begin{equation}\label{2.9} C_i=8\int_{0}^{\infty}t\frac{t^2-1}{(t^2+1)^3}f_i(t)dt, \end{equation} in particular, \begin{equation}\label{2.10} \begin{aligned} \omega_1(y)&=\frac12U_{1,0}^2(y)+6\log(|y|^2+1) +\frac{2\log8-10}{|y|^2+1}+\frac{|y|^2-1}{|y|^2+1}\\ &\quad \times\Big\{-\frac12\log^28 +2\log^2(|y|^2+1) +4\int_{|y|^2}^{\infty}\frac{ds}{s+1} \log\frac{s+1}{s}\\ &\quad -8\log|y|\log(|y|^2+1)\Big\}, \end{aligned} \end{equation} and \begin{equation}\label{2.11} C_1=12-4\log8 \end{equation} (see \cite{CI,EMP}). Our final ansatz for a solution of \eqref{a1} is \begin{equation}\label{2.12} U_\xi(x)=\sum_{j=1}^m\big[U_j(x)+H_j(x)\big], \end{equation} where $H_j$ is a correction term defined as the solution of \begin{equation}\label{2.13} \begin{gathered} -\Delta H_j=0\quad \text{in } \Omega,\\ \frac{\partial H_j}{\partial \nu}+\lambda b(x)H_j = -\frac{\partial U_j}{\partial \nu} -\lambda b(x)U_{j} \text{on } \partial\Omega. \end{gathered} \end{equation} \begin{lemma} \label{lem2.1} For any set of points $\xi=(\xi_1,\ldots,\xi_m)\in\mathcal{O}_\varepsilon$, \begin{equation} \label{2.14} \begin{aligned} H_j(x)&=\frac{1}{\gamma\mu_j^{2/(p-1)}}\Big[ \Big(1-\frac{C_1}{4p}-\frac{C_2}{4p^2} \Big)8\pi H_\lambda(x,\xi_j) -\log(8\mu_j^2\rho^2)\\ &\quad +\Big(\frac{C_1}{p}+\frac{C_2}{p^2} \Big)\log(\mu_j\rho)+ O\big(\frac{\rho}p\big)\Big] \end{aligned} \end{equation} in $C(\overline{\Omega})$ and in $C^2_{\rm loc}(\Omega)$ as $p$ and $\lambda$ go to $+\infty$, where $H_\lambda$ is the regular part of Green's function defined in \eqref{a6}. \end{lemma} \begin{proof} First, on the boundary, by \eqref{2.1} and \eqref{2.8} we have \begin{align*} &\frac{\partial H_j}{\partial \nu}+\lambda b(x)H_j\\ &=-\frac{1}{\gamma\mu_j^{2/(p-1)}}\Big\{ \Big(-4+\frac{C_1}{p}+\frac{C_2}{p^2} \Big)\Big[ \frac{(x-\xi_j)\cdot\nu(x)}{|x-\xi_j|^2} -\lambda b(x) \log\frac1{|x-\xi_j|}\Big] \\ &\quad +\lambda b(x)\Big[\log(8\mu_j^2\rho^2) -\Big(\frac{C_1}{p}+\frac{C_2}{p^2} \Big)\log(\mu_j\rho)\Big] + O\big(\frac{\lambda\rho}{p}\big) \Big\}. \end{align*} The regular part of Green's function with Robin boundary condition $H_\lambda(x,\xi_j)$ satisfies \begin{equation}\label{2.15} \begin{gathered} -\Delta H_\lambda(x,\xi_j)=0 \quad \text{in }\Omega,\\ \frac{\partial H_\lambda(x,\xi_j)}{\partial \nu}+\lambda b(x)H_\lambda(x,\xi_j) =\frac1{2\pi}\frac{(x-\xi_j)\cdot\nu(x)}{|x-\xi_j|^2} -\frac{1}{2\pi}\lambda b(x)\log\frac{1}{|x-\xi_j|}\quad \text{on }\partial\Omega. \end{gathered} \end{equation} So, if we set \begin{align*} \widetilde{H}_j(x)&=\gamma\mu_j^{2/(p-1)}H_j(x)-\Big[ \Big(1-\frac{C_1}{4p}-\frac{C_2}{4p^2} \Big)8\pi H_\lambda(x,\xi_j) \\ &\quad -\log(8\mu_j^2\rho^2)+\Big(\frac{C_1}{p}+\frac{C_2}{p^2} \Big)\log(\mu_j\rho) \Big], \end{align*} then $\widetilde{H}_j(x)$ satisfies \begin{gather*} -\Delta\widetilde{H}_j=0\quad \text{in } \Omega,\\ \frac{\partial \widetilde{H}_j}{\partial \nu}+\lambda b(x)\widetilde{H}_j= O\Big(\frac{\lambda\rho}{p}\Big) \quad\text{on } \partial\Omega. \end{gather*} From the maximum principle with Robin boundary condition (see \cite[Lemma 2.6]{DKM}), we deduce \[ \max_{\overline{\Omega}}\big|\widetilde{H}_j(x)\big|+ \max_{\overline{\Omega}}\big|\operatorname{dist}(x,\partial\Omega)\nabla \widetilde{H}_j(x)\big|\leq \frac{C}{\lambda}\big\|\frac{\partial \widetilde{H}_j}{\partial \nu} +\lambda \widetilde{H}_j\big\|_{L^\infty(\partial\Omega)} =O\Big(\frac{\rho}{p}\Big). \] By the interior estimate of derivative of harmonic function, we derive estimate \eqref{2.14} in $C(\overline{\Omega})$ and in $C^2_{\rm loc}(\Omega)$. \end{proof} From Lemma \ref{lem2.1}, away from the points $\xi_j$, namely $|x-\xi_j|\geq\varepsilon $ for any $j=1,\ldots,m$, one has \begin{equation}\label{2.16} U_\xi(x) =\sum_{j=1}^{m} \frac{1}{\gamma\mu_j^{2/(p-1)}} \Big[\Big(1 -\frac{C_1}{4p}-\frac{C_2}{4p^2}\Big)8\pi G_\lambda(x,\xi_j) +O\Big(\frac{\rho}p \Big)\Big]. \end{equation} While for $|x-\xi_j|\leq\varepsilon $ with some $j$, if we write $x=\xi_j+\delta_j y$, then, by \eqref{2.4}, \eqref{2.5}, \eqref{2.14} and \eqref{2.16} we deduce \begin{align*} &U_\xi(x)\\ &=\frac{1}{\gamma\mu_j^{2/(p-1)}} \Big\{p+ U_{1,0}(y)+\frac1p\omega_1(y) +\frac1{p^2}\omega_2(y) +\Big(1-\frac{C_1}{4p}-\frac{C_2}{4p^2} \Big)8\pi H_\lambda(\xi_j,\xi_j)\\ &\quad -\log(8\mu_j^4) +\Big(\frac{C_1}{p}+\frac{C_2}{p^2} \Big)\log(\mu_j\rho) +O\big(\rho|y|\big) +O\big(\frac{\rho}p \big)\Big\} \\ &\quad +\sum_{k\neq j}\frac{1}{\gamma\mu_k^{2/(p-1)}} \Big[\Big(1 -\frac{C_1}{4p}-\frac{C_2}{4p^2}\Big)8\pi G_\lambda(\xi_j,\xi_k) +O\big(\rho |y|\big) +O\big(\frac{\rho}p\big)\Big]. \end{align*} We now choose the parameters $\mu_j$: we assume they are defined by the relation \begin{align*} \log(8\mu_j^4) &=\Big(1-\frac{C_1}{4p}-\frac{C_2}{4p^2}\Big) \Big[8\pi H_\lambda(\xi_j,\xi_j) +\sum_{k\neq j} \frac{\mu_j^{2/(p-1)}}{\mu_k^{2/(p-1)}} 8\pi G_\lambda(\xi_j,\xi_k)\Big] \\ &\quad +\Big(\frac{C_1}{p}+\frac{C_2}{p^2} \Big) \log\big(\mu_je^{-4p/4}\big). \end{align*} Thus, by the explicit expression \eqref{2.11} of the constant $C_1$, we observe that for $p$ large, the parameters $\mu_j$ satisfies \begin{equation}\label{2.17} \mu_j=e^{-\frac34}e^{2\pi H_\lambda(\xi_j,\xi_j) +2\pi\sum_{k\neq j}G_\lambda(\xi_j,\xi_k)} \big[1+O\big(\frac1p\big)\big]. \end{equation} From this choice of the parameters $\mu_j$, we deduce that for $|x-\xi_j|=\delta_j|y|\leq\varepsilon$, \begin{equation}\label{2.21} U_\xi(x)=\frac{1}{\gamma\mu_j^{2/(p-1)}}\Big[ p+U_{1,0}(y)+\frac1p\omega_1(y)+\frac1{p^2}\omega_2(y)+O( \rho|y|) +O\big(\frac{\rho}p\big)\Big]. \end{equation} \begin{remark} \label{rmk2.2}\rm Let us remark that $U_\xi$ is a positive, uniformly bounded function. Observe that for $|y|\leq\varepsilon/\delta_j$, $$ p+U_{1,0}(y)+\frac1p\omega_1(y)+\frac1{p^2}\omega_2(y) \geq 4\log\frac1{\varepsilon}+ \log(8\mu^4_j)+O\big(\frac{1}{p}\big). $$ Then it is easily checked that choosing $\varepsilon>0$ smaller if necessary, $U_\xi>0$ in $B(\xi_j, \varepsilon )$, and $\sup_{B(\xi_j, \varepsilon )}U_\xi\to\sqrt{e}$ as $p$ and $\lambda$ go to $+\infty$. Moreover, by the maximum principle, we see that $G_\lambda(x,\xi_j)>0$ in $\overline{\Omega}$ and thus by \eqref{2.16}, $U_\xi$ is a positive, uniformly bounded function in $\overline{\Omega}$. In conclusion, $00$ be fixed. There exist $C>0$, $D_0>0$, $p_0>0$ and $\lambda_0>0$ such that \begin{gather}\label{2.27} \|R_\xi\|_{*}\leq C/{p^4}, \\ \label{2.28} W_\xi(x)\leq D_0\sum_{j=1}^{m}e^{U_{\delta_j,\xi_j}(x)}, \end{gather} for any set of point $\xi=(\xi_1,\ldots,\xi_m)\in\mathcal{O}_{\varepsilon}$, any $p\geq p_0$ and $\lambda\geq\lambda_0$. Furthermore, \begin{equation}\label{2.29} W_{\xi}(x)=\frac{8}{\delta_j^2(1+|y|^2)^2} \Big[1+\frac1p\Big(\omega_1-U_{1,0}-\frac12U_{1,0}^2\Big)(y) +O\Big(\frac{\log^4(|y|+2)}{p^2}\Big)\Big], \end{equation} for any $|x-\xi_j|\leq \varepsilon\sqrt{\delta_j}$, where $y=\frac{1}{\delta_j}(x-\xi_j)$. \end{lemma} Since the proof of the above lemma is similar to those of \cite[Prop. 2.1 and Lemma 3.1]{EMP}, we omit it. \section{Linear and nonlinear problems} In this section, we shall study first bounded invertibility of the operator $L$ defined in \eqref{2.23}. Set \begin{equation}\label{3.1**} z_0(y)=\frac{|y|^2-1}{|y|^2+1},\quad z_i(y)=4\frac{y_i}{|y|^2+1},\quad i=1,\,2. \end{equation} It is well known \cite{BP} that any bounded solution to \begin{equation}\label{3.1*} \Delta\phi+\frac{8}{(1+|y|^2)^2}\phi=0 \quad\text{in }\mathbb{R}^2 \end{equation} is a linear combination of $z_i$, $i=0,1,2$. Let us consider the linear problem: given $h\in C(\overline{\Omega})$ and the set of points $\xi=(\xi_1,\ldots,\xi_m)\in\mathcal{O}_\varepsilon$, we find a function $\phi$ and scalars $c_{ij}$, $i=1,2$, $j=1,\ldots,m$, such that \begin{equation}\label{3.1} \begin{gathered} L(\phi)= \Delta\phi+W_{\xi}\phi =h+\sum_{i=1}^2\sum_{j=1}^mc_{ij}e^{U_{\delta_j,\xi_j}}Z_{ij} \quad \text{in } \Omega,\\ \frac{\partial \phi}{\partial\nu}+\lambda b(x)\phi=0 \quad \text{on } \partial\Omega,\\ \int_{\Omega}e^{U_{\delta_j,\xi_j}}Z_{ij}\phi=0 \quad \text{for }i=1,2;\; j=1,\ldots,m. \end{gathered} \end{equation} Here, for $i=0,1,2$ and $j=1,\ldots,m$, we denote \begin{equation}\label{3.2} Z_{ij}(x):=z_i\Big(\frac{x-\xi_j}{\delta_j}\Big) =\begin{cases} \frac{|x-\xi_j|^2-\delta_j^2}{|x-\xi_j|^2+\delta_j^2} & \text{if }i=0,\\[4pt] \frac{4\delta_j(x-\xi_j)_i}{|x-\xi_j|^2+\delta_j^2} &\text{if }i=1,2. \end{cases} \end{equation} \begin{proposition} \label{prop3.1} Let $\varepsilon>0$ be fixed. There exist $p_0>0$, $\lambda_0>0$ and $C>0$ such that for any $h\in C(\overline{\Omega})$, any the set of points $\xi=(\xi_1,\ldots,\xi_m)\in\mathcal{O}_\varepsilon$, any $p>p_0$ and $\lambda>\lambda_0$, there is a unique solution $\phi$, scalars $c_{ij}$, $i=1,2$, $j=1,\ldots,m$, to problem \eqref{3.1}, which satisfies \begin{equation}\label{3.3} \|\phi\|_{\infty}\leq Cp\|h\|_{*}. \end{equation} \end{proposition} \begin{proof} The proof of this result will be divided into six steps. \smallskip \noindent\textbf{Step 1:} The operator $L$ satisfies the maximum principle in $\widetilde{\Omega}:=\Omega\setminus\cup_{j=1}^mB(\xi_j,R\delta_j)$ for $R$ large, independent on $p$ and $\lambda$. Specifically, if $\psi$ satisfies \begin{gather*} L(\psi)= \Delta\psi+W_{\xi}\psi\leq0 \quad \text{in }\widetilde{\Omega},\\ \psi\geq0 \quad\text{on }\cup_{j=1}^m\partial B(\xi_j,R\delta_j) \quad\text{and}\quad \frac{\partial \psi}{\partial\nu}+\lambda b(x)\psi\geq0 \quad\text{on }\partial\Omega, \end{gather*} then $\psi\geq0$ in $\widetilde{\Omega}$. To prove this, it suffices to construct a positive function $Z$ on $\widetilde{\Omega}$ such that \begin{gather*} L(Z)= \Delta Z+W_{\xi}Z <0\quad \text{in }\widetilde{\Omega}, \\ Z>0 \quad \text{on }\cup_{j=1}^m\partial B(\xi_j,R\delta_j) \quad\text{and}\quad \frac{\partial Z}{\partial\nu}+\lambda b(x)Z>0 \quad \text{on }\partial\Omega. \end{gather*} Indeed, let $$ Z(x)=\sum_{j=1}^mz_0\Big( \frac{a(x-\xi_j)}{\delta_j} \Big),\quad a>0. $$ First, observe that, if $|x-\xi_j|\geq R\delta_j$ for $R>\frac1a$, then $Z(x)>0$. On the other hand, since $Z(x)\leq m$, $$ W_\xi(x)Z(x)\leq D_0Z(x)\sum_{j=1}^{m}e^{U_{\delta_j,\xi_j}(x)} \leq D_0Z(x)\sum_{j=1}^m\frac{8\delta_j^2}{|x-\xi_j|^4}\leq mD_0\sum_{j=1}^m\frac{8\delta_j^2}{|x-\xi_j|^4}, $$ where $D_0$ is the constant in Lemma \ref{lem2.3}. Further, by the definition of $z_0$, \begin{align*} -\Delta Z(x) &=\sum_{j=1}^ma^2\frac{8\delta_j^2(a^2|x-\xi_j|^2 -\delta_j^2)}{(a^2|x-\xi_j|^2+\delta_j^2)^3}\\ &\geq\frac13\sum_{j=1}^m\frac{8a^2\delta_j^2}{(a^2|x-\xi_j|^2 +\delta_j^2)^2} \geq\frac4{27} \sum_{j=1}^m\frac{8\delta_j^2}{a^2|x-\xi_j|^4}, \end{align*} provided $R>\sqrt{3}/a$. Thus, if $a$ is taken small and fixed, but independent of $p$ and $\lambda$, and $R$ is chosen sufficiently large depending on this $a$, then we have that $$ L(Z)=\Delta Z+W_{\xi}Z\leq\Big(-\frac4{27}\frac1{a^2}+mD_0\Big) \sum_{j=1}^m\frac{8\delta_j^2}{|x-\xi_j|^4}<0. $$ Moreover, \begin{gather*} \big|\frac{\partial}{\partial\nu} Z(x)\big| \leq \sum_{j=1}^m\frac{C\delta_j^2}{a^2|x-\xi_j|^3} =O\Big(\frac{\rho^2}{a^2\varepsilon^3}\Big)\quad \text{on }\partial\Omega, \\ Z(x)\geq\frac12 \quad\text{on } \partial\Omega\cup\big(\cup_{j=1}^m\partial B(\xi_j,R\delta_j)\big), \end{gather*} which, together with \eqref{2.4}, imply that on $\partial\Omega$, \begin{equation}\label{3.4} \frac{\partial Z}{\partial\nu}+\lambda b(x)Z \geq O\Big(\frac{1}{a^2\varepsilon^3}\rho^2\Big) +\frac12\lambda b(x)\geq O\big(e^{-p/2}\big)+\frac12\lambda\min_{x\in\partial\Omega} b(x)>0 \end{equation} provided that $p$ is chosen sufficiently large. The function $Z(x)$ is what we want. \smallskip \noindent\textbf{Step 2:} Let $R$ be as before. We define the ``inner norm'' of $\phi$ as $$ \|\phi\|_i=\sup_{x\in\cup_{j=1}^mB(\xi_j,R\delta_j)}|\phi(x)| $$ and claim that there is a constant $C>0$ such that if $L(\phi)=h$ in $\Omega$, $\frac{\partial \phi}{\partial\nu}+\lambda b(x)\phi=g$ on $\partial\Omega$, then $$ \|\phi\|_{L^{\infty}(\Omega)} \leq C\Big(\|\phi\|_i+ \|h\|_{*}+\frac1{\lambda}\|g\|_{L^{\infty}(\partial\Omega)} \Big), $$ for any $h\in C^{0,\alpha}(\overline{\Omega})$ and $g\in C^{0,\alpha}(\partial\Omega)$. We will establish this estimate with the use of suitable barriers. Let $M=2\operatorname{diam}\Omega$. Consider the problem \begin{gather*} -\Delta\psi_j=\frac{2\delta_j}{|x-\xi_j|^3} \quad \text{in } R\delta_j<|x-\xi_j|2$ is chosen larger if necessary. First of all, observe that for $x\in\cup_{j=1}^m\partial B(\xi_j,R\delta_j)$, by the definition of $Z$, $$ \widetilde{\phi}(x) \geq 2C_0\|\phi\|_iZ(x)\geq\|\phi\|_i\geq|\phi(x)|, $$ for $x\in\partial\Omega$, by \eqref{3.4}, \eqref{3.5} and the positivity of $Z(x)$ and $\psi_j(x)$, \begin{align*} &\frac{\partial \widetilde{\phi}}{\partial\nu}(x) +\lambda b(x)\widetilde{\phi}(x)\\ &\geq \Big[O\Big( e^{-\frac{1}{2}p} +\frac{1}{p}\Big) +C_0\lambda\min_{x\in\partial\Omega} b(x)\Big]\Big(\|\phi\|_i +\|h\|_{*}+\frac1{\lambda}\|g\|_{L^{\infty}(\partial\Omega)}\Big)\\ &\geq \frac12C_0\lambda\Big(\min_{x\in\partial\Omega} b(x)\Big) \Big(\|\phi\|_i+ \|h\|_{*}+\frac1{\lambda}\|g\|_{L^{\infty}(\partial\Omega)} \Big) \geq |g(x)|, \end{align*} and for $x\in\Omega\setminus\cup_{j=1}^mB(\xi_j,R\delta_j)$, by \eqref{2.28}, \eqref{3.5*} and the definition of $\|\cdot\|_*$ in \eqref{2.26}, \begin{align*} L(\widetilde{\phi}) &\leq C_0\Big(\|\phi\|_i+ \|h\|_{*}+\frac1{\lambda}\|g\|_{L^{\infty}(\partial\Omega)} \Big) \sum_{j=1}^m L(\psi_j)(x)\\ &=C_0\Big(\|\phi\|_i+ \|h\|_{*}+\frac1{\lambda}\|g\|_{L^{\infty}(\partial\Omega)} \Big)\sum_{j=1}^m\big(-\frac{2\delta_j}{|x-\xi_j|^3} +W_\xi(x)\psi_j(x) \big)\\ &\leq C_0\Big(\|\phi\|_i+ \|h\|_{*}+\frac1{\lambda}\|g\|_{L^{\infty}(\partial\Omega)} \Big) \sum_{j=1}^m\Big(-\frac{2\delta_j}{|x-\xi_j|^3} +\frac{2mD_0}{R}e^{U_{\delta_j,\xi_j}(x)} \Big)\\ &\leq -C_0 \|h\|_{*} \sum_{j=1}^m\frac{\delta_j}{(\delta_j^2+|x-\xi_j|^2)^{3/2}}\\ &\leq-|h(x)|\leq-|L(\phi)|(x), \end{align*} provided $R>16mD_0$, $p$ and $\lambda$ large enough. Hence, by the maximum principle in Step 1, we obtain $$ |\phi(x)|\leq\widetilde{\phi}(x) \quad\text{for }x\in\widetilde{\Omega}, $$ and therefore, since $Z(x)\leq m$ and $\psi_j(x)\leq\frac2{R}$, $$ \|\phi\|_{L^{\infty}(\Omega)} \leq C\Big(\|\phi\|_i+ \|h\|_{*}+\frac1{\lambda}\|g\|_{L^{\infty}(\partial\Omega)} \Big). $$ \smallskip \noindent\textbf{Step 3:} We prove uniform a priori estimates for solutions $\phi$ of the problem $L(\phi)=h$ in $\Omega$, $\frac{\partial \phi}{\partial\nu}+\lambda b(x)\phi=g$ on $\partial\Omega$, where $h\in C^{0,\alpha}(\overline{\Omega})$, $g\in C^{0,\alpha}(\partial\Omega)$ and in addition we prove the orthogonality conditions: \begin{equation}\label{3.14} \int_{\Omega}e^{U_{\delta_j,\xi_j}}Z_{ij}\phi=0 \quad \text{for }i=0,1,2;\; j=1,\ldots,m. \end{equation} Namely, we prove that there exists $C>0$ such that for $\xi=(\xi_1,\ldots,\xi_m)\in\mathcal{O}_\varepsilon$, $h\in C^{0,\alpha}(\overline{\Omega})$ and $g\in C^{0,\alpha}(\partial\Omega)$, $$ \|\phi\|_{L^{\infty}(\Omega)}\leq C\Big( \|h\|_{*}+\frac1{\lambda}\|g\|_{L^{\infty}(\partial\Omega)} \Big), $$ for $p$ and $\lambda$ sufficiently large. By contradiction, assume the existence of sequences $p_n\to+\infty$, $\lambda_n\to+\infty$, points $\xi^n=(\xi_1^n,\ldots,\xi_m^n)\in\mathcal{O}_\varepsilon$, functions $h_n$, $g_n$ and associated solutions $\phi_n$ such that $\|h_n\|_{*}\to0$, $\frac1{\lambda_n}\|g_n\|_{L^{\infty}(\partial\Omega)}\to0$ and $\|\phi_n\|_{L^{\infty}(\Omega)}=1$. Since $\|\phi_n\|_{L^{\infty}(\Omega)}=1$, Step 2 shows that $\liminf_{n\to+\infty}\|\phi_n\|_i>0$. Set $\widehat{\phi}^n_j(y)=\phi_n(\delta_j^ny+\xi_j^n)$ for $j=1,\ldots,m$. By \eqref{2.29} elliptic estimates imply that $\widehat{\phi}^n_j$ converges uniformly over compact sets to a bounded solution $\widehat{\phi}^{\infty}_j$ of equation \eqref{3.1*}. Furthermore, $\widehat{\phi}^{\infty}_j$ is a linear combination of the functions $z_i$, $i=0,1,2$, defined in \eqref{3.1**}. Since $\|\widehat{\phi}^n_j\|_{L^{\infty}(\Omega)}\leq1$, by Lebesgue's theorem, the orthogonality conditions on $\widehat{\phi}^n_j$ pass to the limit and give $$ \int_{\mathbb{R}^2}\frac{8}{(1+|y|^2)^2}z_i(y)\widehat{\phi}_j^{\infty}dy=0 \quad\text{for }i=0,1,2. $$ Hence, $\widehat{\phi}_j^{\infty}\equiv0$ for any $j=1,\ldots,m$ contradicting $\liminf_{n\to+\infty}\|\phi_n\|_i>0$. \smallskip \noindent\textbf{Step 4:} We prove that there exists a positive constant $C>0$ such that any solution $\phi$ of equation $L(\phi)=h$ in $\Omega$, $\frac{\partial \phi}{\partial\nu}+\lambda b(x)\phi=0$ on $\partial\Omega$ and in addition the orthogonality conditions: \begin{equation}\label{3.15} \int_{\Omega}e^{U_{\delta_j,\xi_j}}Z_{ij}\phi=0 \quad \text{for }i=1,2;\; j=1,\ldots,m, \end{equation} satisfies $$ \|\phi\|_{L^{\infty}(\Omega)}\leq Cp \|h\|_{*}, $$ for $h\in C^{0,\alpha}(\overline{\Omega})$. Proceeding by contradiction as in Step 3, we can suppose further that \begin{equation}\label{3.6} \|\phi_n\|_{L^{\infty}(\Omega)}=1, \quad p_n\|h_n\|_{*}\to 0 \quad\text{as }n\to+\infty. \end{equation} but we lose the condition $\int_{\mathbb{R}^2}\frac{8}{(1+|y|^2)^2}z_0(y)\widehat{\phi}_j^{\infty}=0$ in the limit. Hence, we have that \begin{equation}\label{3.7} \widehat{\phi}_j^{n}\to \widehat{\phi}_j^{\infty}=C_j\frac{|y|^2-1}{|y|^2+1} \quad \text{in }C_{\rm loc}^0(\mathbb{R}^2_+), \end{equation} for some constants $C_j$. To reach a contradiction, we have to show that $C_j=0$ for any $j=1,\ldots,m$. We will obtain it from the stronger condition \eqref{3.6} on $h_n$. To this end, we perform the following construction. According to \cite{CI,EMP}, there exist radial solutions $\omega$ and $\zeta$ respectively of equations $$ \Delta\omega+\frac{8}{(1+|y|^2)^2}\omega= \frac{8}{(1+|y|^2)^2}z_0(y), \quad \Delta\zeta+\frac{8}{(1+|y|^2)^2}\zeta= \frac{8}{(1+|y|^2)^2} \quad \text{in }\mathbb{R}^2, $$ such that \begin{gather*} \omega(y)= \frac43\log|y|+O\Big(\frac1{|y|}\Big),\quad \zeta(y)= O\big(\frac1{|y|}\big)\quad\text{as }|y|\to+\infty, \\ \nabla\omega(y) =\frac43\cdot\frac{y}{1+|y|^2}+ O\big(\frac1{1+|y|^2}\big),\quad \nabla\zeta(y)= O\big(\frac1{1+|y|^2}\big) \quad\text{for all } y\in\mathbb{R}^2, \end{gather*} since $8\int_0^{+\infty}r\frac{(r^2-1)^2}{(r^2+1)^4}dr=\frac43$ and $8\int_0^{+\infty}r\frac{r^2-1}{(r^2+1)^3}dr=0$. For simplicity in the rest of this article, we omit the dependence on $n$. For $j=1,\ldots,m$, define $$ u_j(x)=\omega\Big(\frac{x-\xi_j}{\delta_j}\Big) +\frac43(\log\delta_j)Z_{0j}(x) +\frac{8\pi}{3}H_\lambda(\xi_j,\xi_j) \zeta\Big(\frac{x-\xi_j}{\delta_j}\Big) $$ and denote its projection $Pu_j=u_j+\widetilde{H}_j$, where $\widetilde{H}_j$ is a correction term defined as the solution of \begin{gather*} -\Delta\widetilde{H}_j=0 \quad\textrm {in }\Omega,\\ \frac{\partial\widetilde{H}_j}{\partial \nu}+\lambda b(x)\widetilde{H}_j= -\frac{\partial u_j}{\partial \nu} -\lambda b(x)u_j \quad \text{on } \partial\Omega. \end{gather*} Observe that on $\partial\Omega$, $$ \Big(\frac{\partial}{\partial \nu}+\lambda b(x)\Big) \Big(\widetilde{H}_j+\frac{8\pi}3H_\lambda(x,\xi_j)\Big) =O\big(\lambda\rho\big) +(\log\delta_j)O\big(\lambda\rho^2\big) +H_\lambda(\xi_j,\xi_j)O\big(\lambda\rho\big). $$ From the maximum principle with Robin boundary condition we obtain \begin{equation}\label{3.8} \begin{gathered} Pu_j=u_j-\frac{8\pi}{3}H_\lambda(x,\xi_j) +O\big(\rho\big) \quad \text{in }C(\overline{\Omega}),\\ Pu_j=-\frac{8\pi}{3}G_\lambda(x,\xi_j) +O(\rho) \quad\text{in } C_{\rm loc}(\overline{\Omega}\setminus\{\xi_j\}). \end{gathered} \end{equation} The function $Pu_j$ solves \begin{equation}\label{3.9} \begin{gathered} \Delta Pu_j+W_\xi Pu_j=e^{U_{\delta_j,\xi_j}}Z_{0j} +\big(W_\xi-e^{U_{\delta_j,\xi_j}}\big)Pu_j+R_j \quad \text{in }\Omega,\\ \frac{\partial }{\partial \nu}Pu_j+\lambda b(x)Pu_j=0 \quad \text{on }\partial\Omega, \end{gathered} \end{equation} where \begin{equation}\label{3.9*} R_j(x)=\Big(Pu_j-u_j +\frac{8\pi}{3}H_\lambda(\xi_j,\xi_j) \Big)e^{U_{\delta_j,\xi_j}}. \end{equation} Multiplying \eqref{3.9} by $\phi$ and integrating by parts we obtain \begin{equation}\label{3.10} \int_{\Omega}e^{U_{\delta_j,\xi_j}}Z_{0j}\phi+ \int_{\Omega}\big(W_\xi-e^{U_{\delta_j,\xi_j}}\big)Pu_j\phi =\int_{\Omega}Pu_jh-\int_{\Omega}R_j\phi. \end{equation} We estimate each term of \eqref{3.10}. First of all, by Lebesgue's theorem and \eqref{3.7} we obtain \begin{equation}\label{3.11} \int_{\Omega}e^{U_{\delta_j,\xi_j}}Z_{0j}\phi \longrightarrow C_j\int_{\mathbb{R}^2}\frac{8(|y|^2-1)^2}{(|y|^2+1)^4}dy=\frac{8\pi}{3}C_j. \end{equation} From \eqref{2.15} and the maximum principle with Robin boundary condition, we deduce that $|\nabla H_\lambda(x,\xi_j)|=O(1)$ holds uniformly in $\Omega$. Thus, by \eqref{2.28}, \eqref{2.29} and \eqref{3.8}, we have \begin{align*} &\int_{\Omega}\big(W_\xi-e^{U_{\delta_j,\xi_j}}\big)Pu_j\phi\\ &= \int_{B(\xi_j,\varepsilon\sqrt{\delta_j})} \big(W_\xi-e^{U_{\delta_j,\xi_j}}\big)Pu_j\phi -\frac{8\pi}{3}\sum_{k\neq j}G_\lambda(\xi_k,\xi_j) \int_{B(\xi_k,\varepsilon \sqrt{\delta_k})}W_\xi\phi +O\big(\sqrt{\rho} \big)\\ &=\int_{B(0,\varepsilon/\sqrt{\delta_j})} \frac{8}{(1+|y|^2)^2}\frac1p\Big(\omega_1-U_{1,0}-\frac12U_{1,0}^2\Big) \frac43(\log\delta_j)z_0(y)\widehat{\phi}_j \\ &\quad -\frac{8\pi}{3}\sum_{k\neq j}G_\lambda(\xi_k,\xi_j) \int_{B(0,\varepsilon/\sqrt{\delta_k})}\frac{8}{(1+|y|^2)^2}\widehat{\phi}_k +O\big(\frac1p\big)\\ &=-\frac{C_j}{3}\int_{\mathbb{R}^2} \frac{8(|y|^2-1)^2}{(1+|y|^2)^4}\Big(\omega_1-U_{1,0}-\frac12U_{1,0}^2\Big)(y)+o(1)\,. \end{align*} Lebesgue's theorem and \eqref{3.7} imply \begin{align*} &\int_{B(0,\varepsilon/\sqrt{\delta_j})} \frac{8}{(1+|y|^2)^2}\Big(\omega_1-U_{1,0}-\frac12U_{1,0}^2\Big) z_0(y)\widehat{\phi}_j \\ &\to C_j\int_{\mathbb{R}^2}\frac{8(|y|^2-1)^2}{(1+|y|^2)^4} \Big(\omega_1-U_{1,0}-\frac12U_{1,0}^2\Big), \end{align*} and $$ \int_{B(0,\varepsilon/\sqrt{\delta_k})}\frac{8}{(1+|y|^2)^2}\widehat{\phi}_k \to C_k\int_{\mathbb{R}^2}\frac{8}{(1+|y|^2)^2}\frac{|y|^2-1}{|y|^2+1}=0. $$ In a straightforward but tedious comptation, by \eqref{2.10} we obtain $$ \int_{\mathbb{R}^2}\frac{8(|y|^2-1)^2}{(1+|y|^2)^4} \Big(\omega_1-U_{1,0}-\frac12U_{1,0}^2\Big)(y)dy=-8\pi. $$ Therefore \begin{equation}\label{3.12} \int_{\Omega}\big(W_\xi-e^{U_{\delta_j,\xi_j}}\big)Pu_j\phi =\frac{8\pi}{3}C_j+o(1). \end{equation} As for the right-hand side of \eqref{3.10}, that by \eqref{2.26} and \eqref{3.8}, we have \begin{equation}\label{3.13} \begin{aligned} \big|\int_{\Omega}Pu_jh\big| &\leq C\|h\|_{*}\sum_{k=1}^m \int_{\Omega_{\delta_k}} \frac{1}{(1+|y|^2)^{3/2}}|Pu_j(\delta_k y+\xi_k)| dy \\ &\leq C\|h\|_{*} \int_{\mathbb{R}^2} \frac{\log(|y|+2)}{(1+|y|^2)^{3/2}} dy + Cp\|h\|_{*} \int_{\mathbb{R}^2} \frac{dy}{(1+|y|^2)^{3/2}}\\ &\leq Cp\|h\|_{*}, \end{aligned} \end{equation} where $\Omega_{\delta_k}:=\frac{1}{\delta_k}(\Omega-\{\xi_k\})$. Finally, by \eqref{3.8} and \eqref{3.9*} we deduce \begin{equation}\label{3.16} \int_{\Omega}R_j\phi=O\Big(\int_{\Omega} e^{U_{\delta_j,\xi_j}} \big(|x-\xi_j|+\rho\big)dx\Big) =O\big(e^{-\frac14p}\big). \end{equation} Hence, inserting \eqref{3.11}-\eqref{3.16} in \eqref{3.10} and taking into account \eqref{3.6}, we conclude that $$ \frac{16\pi}{3}C_j=o(1) \quad \text{for any }j=1,\ldots,m. $$ Necessarily, $C_j=0$ by contradiction and the claim is proved. \smallskip \noindent\textbf{Step 5:} We establish the validity of the a priori estimate \begin{equation}\label{3.17} \|\phi\|_{\infty}\leq Cp\|h\|_{*} \end{equation} for solutions of problem \eqref{3.1} and $h\in C^{0,\alpha}(\overline{\Omega})$. Step 4 gives $$ \|\phi\|_{L^{\infty}(\Omega)} \leq Cp\Big( \|h\|_{*}+\sum_{i=1}^2\sum_{j=1}^m|c_{ij}|\cdot\|e^{U_{\delta_j,\xi_j}}Z_{ij}\|_* \Big) \leq Cp\Big(\|h\|_{*}+\sum_{i=1}^2\sum_{j=1}^m|c_{ij}|\Big). $$ As before, arguing by contradiction of \eqref{3.17}, we can proceed as in Step 3 and suppose further that \begin{equation}\label{3.18} \|\phi_n\|_{L^{\infty}(\Omega)}=1, \quad p_n\|h_n\|_{*}\to0, \quad p_n\sum_{i=1}^2\sum_{j=1}^m|c_{ij}^n|\geq\delta>0 \quad \text{as }n\to+\infty. \end{equation} We omit the dependence on $n$. It suffices to estimate the values of the constants $c_{ij}$. For this aim, we define $PZ_{ij}$ as the projection of $Z_{ij}$ under homogeneous Robin boundary condition, namely \begin{equation}\label{3.19} \begin{gathered} \Delta PZ_{ij}=\Delta Z_{ij}=-e^{U_{\delta_j,\xi_j}}Z_{ij} \text{in }\Omega,\\ \frac{\partial PZ_{ij}}{\partial\nu}+\lambda b(x)PZ_{ij}=0 \quad\text{on }\partial\Omega. \end{gathered} \end{equation} As in the proof of Lemma \ref{lem2.1}, for $i=1,2$ and $j=1,\ldots,m$ we have the expansions: \begin{equation}\label{3.20} PZ_{ij}=Z_{ij}+8\pi\delta_j\partial_{(\xi_j)_i}H_\lambda(\cdot,\xi_j) +O\big(\rho^3\big), \quad PZ_{0j}=Z_{0j}-1 +O\big(\rho^2\big), \end{equation} in $C(\overline{\Omega})$ and in $C^2_{\rm loc}(\Omega)$, and \begin{equation}\label{3.22} PZ_{ij}=8\pi\delta_j\partial_{(\xi_j)_i}G_\lambda(\cdot,\xi_j) +O(\rho^3), \quad PZ_{0j}=O(\rho^2), \end{equation} in $C(\overline{\Omega}\setminus\{\xi_j\})$ and in $C^2_{\rm loc}(\Omega\setminus\{\xi_j\})$. By \eqref{3.20}, \eqref{3.22} and that $|\partial_{(\xi_j)_i}H_\lambda(x,\xi_j)|=O(1)$ uniformly holds in $\Omega$, we can easily deduce the following ``orthogonality'' relations: for each $i,l=1,2$ and $j,k=1,\ldots,m$, \begin{equation}\label{3.23} \int_{\Omega}e^{U_{\delta_j,\xi_j}}Z_{ij}PZ_{lk} =\Big(64\int_{\mathbb{R}^2} \frac{|y|^2}{(1+|y|^2)^4}\Big)\delta_{jk}\delta_{il} +O\big(\rho\big), \end{equation} uniformly for any set of points $\xi=(\xi_1,\ldots,\xi_m)\in\mathcal{O}_\varepsilon$, where $\delta_{jk}$ and $\delta_{il}$ denote the Kronecker's symbols. Multiplying equation \eqref{3.1} by $PZ_{ij}$, $i=1,2$, $j=1,\ldots,m$, and integrating by parts we find \begin{equation}\label{3.25} \sum_{l=1}^2\sum_{k=1}^mc_{lk}\int_{\Omega}e^{U_{\delta_k,\xi_k}}Z_{lk} PZ_{ij}+\int_{\Omega}hPZ_{ij}= \int_{\Omega}W_{\xi}\phi PZ_{ij}-\int_{\Omega}e^{U_{\delta_j,\xi_j}}Z_{ij}\phi. \end{equation} By \eqref{2.28}, \eqref{2.29} and \eqref{3.23}, a direct computation shows \begin{equation} \begin{aligned} &Dc_{ij}+O\Big(e^{-\frac{p}{2}} \sum_{l,\,k}|c_{lk}|+\|h\|_*\Big)\\ &=\frac1p \int_{B(0,\varepsilon/\sqrt{\delta_j})}\frac{32y_i}{(1+|y|^2)^3} \Big(\omega_1-U_{1,0}-\frac12U_{1,0}^2\Big)\widehat{\phi}_j +O\Big(\frac{\|\phi\|_{\infty}}{p^2}\Big), \end{aligned} \end{equation} where $D=64\int_{\mathbb{R}^2} \frac{|y|^2}{(1+|y|^2)^4}$ and $\widehat{\phi}_j(y)=\phi(\delta_jy+\xi_j)$. Hence, we obtain \begin{equation}\label{3.28} \sum_{l=1}^2\sum_{k=1}^m|c_{lk}| =O\Big(\|h\|_{*}+\frac{1}{p}\|\phi\|_{\infty}\Big)=o(1). \end{equation} As in Step 4, we conclude that for each $j=1,\ldots,m$, $$ \widehat{\phi}_j\to C_j\frac{|y|^2-1}{|y|^2+1} \quad \text{in }C_{\rm loc}^0(\mathbb{R}^2), $$ with some constant $C_j\in\mathbb{R}$ and thus \begin{align*} &\int_{B(0,\varepsilon/\sqrt{\delta_j})}\frac{32y_i}{(1+|y|^2)^3} \Big(\omega_1-U_{1,0}-\frac12U_{1,0}^2\Big)\widehat{\phi}_j \\ &\to C_j\int_{\mathbb{R}^2}\frac{32y_i(|y|^2-1)}{(1+|y|^2)^4} \Big(\omega_1-U_{1,0}-\frac12U_{1,0}^2\Big)=0. \end{align*} Therefore, $$ \sum_{l=1}^2\sum_{k=1}^m|c_{lk}|=o(\frac1p)+O(\|h\|_{*}), $$ which is impossible because of \eqref{3.18}. \smallskip \noindent\textbf{Step 6:} We prove the solvability of \eqref{3.1}. To this purpose, we consider the spaces: \begin{gather*} K_{\xi}=\big\{\sum_{i=1}^2\sum_{j=1}^mc_{ij}PZ_{ij}: c_{ij}\in\mathbb{R}\text{ for } i=1,2;\; j=1,\ldots,m \big\},\\ K_{\xi}^{\bot}=\big\{\phi\in L^2(\Omega): \int_{\Omega}e^{U_{\delta_j,\xi_j}}Z_{ij}\phi=0 \text{ for } i=1,2;\; j=1,\ldots,m \big\}. \end{gather*} Define $\Pi_{\xi}: L^2(\Omega)\to K_{\xi}$ by $$ \Pi_{\xi}\phi=\sum_{i=1}^2\sum_{j=1}^mc_{ij}PZ_{ij}, $$ where the coefficients $c_{ij}$ are uniquely determined (as it follows by \eqref{3.23}) by the system $$ \int_{\Omega}e^{U_{\delta_k,\xi_k}}Z_{lk}\Big( \phi-\sum_{i=1}^2\sum_{j=1}^mc_{ij}PZ_{ij} \Big)=0 \quad \text{for any } l=1,2;\; k=1,\ldots,m. $$ Let $\Pi_{\xi}^\bot=Id-\Pi_{\xi}: L^2(\Omega)\to K^\bot_\xi$. Moreover, the Hilbert space $K_{\xi}^{\bot}\cap H^1(\Omega)$ is endowed with the inner product $$ \langle\phi,\,\psi\rangle_H = \int_{\Omega}\nabla\phi\nabla\psi +\lambda\int_{\partial\Omega} b(x)\phi\psi. $$ Problem \eqref{3.1}, expressed in a weak form, is equivalent to find $\phi\in K_{\xi}^{\bot}\cap H^1(\Omega)$ such that $$ \langle\phi,\,\psi\rangle_H=\int_{\Omega}\left(W_\xi\phi -h\right)\psi \quad \text{for all } \psi\in K_{\xi}^{\bot}\cap H^1(\Omega). $$ With the aid of Riesz's representation theorem, this equation gets rewritten in $K_{\xi}^{\bot}\cap H^1(\Omega)$ in the operator form $\phi=K(\phi)+\tilde{h}$, where \begin{gather*} \tilde{h}=-\Pi_{\xi}^\bot\big[\big(-\Delta\big)|_{\Omega}+ \big(\frac{\partial}{\partial\nu}+\lambda b(x)\big)|_{\partial\Omega}\big]^{-1}h, \\ K(\phi)=\Pi_{\xi}^\bot\big[\big(-\Delta\big)|_{\Omega}+ \big(\frac{\partial}{\partial\nu}+\lambda b(x)\big)|_{\partial\Omega} \big]^{-1}(W_\xi\phi) \end{gather*} is a linear compact operator in $K_{\xi}^{\bot}\cap H^1(\Omega)$. By the Fredholm's alternative with Robin boundary condition (see \cite{D,GT}), we obtain the unique solvability of this problem for any $\tilde{h}\in K^\bot_\xi$ provided that the homogeneous equation $\phi=K(\phi)$ has only the trivial solution in $K_{\xi}^{\bot}\cap H^1(\Omega)$, which in turn follows from the a priori estimate \eqref{3.17} in Step 5. Finally, by density we obtain the validity of \eqref{3.3} also for $h\in C(\overline{\Omega})$ (not only for $h\in C^{0,\alpha}(\overline{\Omega})$). \end{proof} \begin{remark} \label{rmk3.2} \rm Given $h\in C(\overline{\Omega})$, let $\phi$ be the solution of problem \eqref{3.1} given by Proposition \ref{prop3.1}. Multiplying \eqref{3.1} against $\phi$ and integrating by parts, we obtain $$ \|\phi\|_H^2:=\int_{\Omega}|\nabla\phi|^2 +\lambda\int_{\partial\Omega} b(x)\phi^2=\int_{\Omega}W_\xi\phi^2-\int_{\Omega}h\phi. $$ By Lemma \ref{lem2.3}, we obtain $$ \|\phi\|_{H}\leq C\left(\|h\|_{*}+\|\phi\|_{\infty}\right). $$ \end{remark} Let us solve the nonlinear auxiliary problem: for any set of points $\xi=(\xi_1,\ldots,\xi_m)$ in $\mathcal{O}_\varepsilon$, we find a function $\phi$, and scalars $c_{ij}$, $i=1,2$, $j=1,\ldots,m$, such that \begin{equation}\label{4.1} \begin{gathered} \Delta(U_\xi+\phi)+(U_\xi+\phi)^p= \sum_{i=1}^2\sum_{j=1}^mc_{ij}e^{U_{\delta_j,\xi_j}}Z_{ij} \quad \text{in } \Omega,\\ U_\xi+\phi>0 \quad \text{in }\Omega,\\ \frac{\partial \phi}{\partial\nu}+\lambda b(x)\phi=0 \quad \text{on } \partial\Omega,\\ \int_{\Omega}e^{U_{\delta_j,\xi_j}}Z_{ij}\phi=0 \quad \text{for } i=1,2;\; j=1,\ldots,m. \end{gathered} \end{equation} \begin{proposition} \label{prop3.3} Let $\varepsilon>0$ be fixed and small. There exist $C>0$, $p_0>0$ and $\lambda_0>0$ such that for any set of points $\xi=(\xi_1,\ldots,\xi_m)\in\mathcal{O}_\varepsilon$, any $p>p_0$ and $\lambda>\lambda_0$, problem \eqref{4.1} has a unique solution $\phi_\xi$, scalars $c_{ij}(\xi)$, $i=1,2$, $j=1,\ldots,m$, such that \begin{equation}\label{4.2} \|\phi_\xi\|_{\infty}\leq\frac{C}{p^3}, \quad \sum_{i=1}^2\sum_{j=1}^m|c_{ij}(\xi)|\leq\frac{C}{p^4},\quad \|\phi_\xi\|_{H}\leq\frac{C}{p^3}. \end{equation} Furthermore, the map $\xi\to\phi_\xi$ is a $C^1$-function in $C(\overline{\Omega})$ and $H^1(\Omega)$. \end{proposition} \begin{proof} Proposition \ref{prop3.1} allows us to apply the contraction mapping principle to find a solution for problem \eqref{4.1} satisfying \eqref{4.2}. Since it is a standard procedure, we shall not present the detailed proof, see \cite[Lemma 4.1]{EMP}. \end{proof} \section{Variational reduction} After problem \eqref{4.1} has been solved, we find a solution of \eqref{2.22} and hence to the original problem \eqref{a1} if $\xi=(\xi_1,\ldots,\xi_m)\in\mathcal{O}_\varepsilon$ satisfies \begin{equation}\label{5.1} c_{ij}(\xi)=0\quad \text{for all } i=1,2;\; j=1,\ldots,m. \end{equation} Equation \eqref{a1} is the Euler-Lagrange equation of the functional $J^\lambda_p: H^1(\Omega)\to\mathbb{R}$ defined by \begin{equation}\label{5.2b} J^\lambda_p(u)=\frac1{2}\int_{\Omega}|\nabla u|^2-\frac1{p+1}\int_{\Omega}u_+^{p+1} +\frac{\lambda}{2}\int_{\partial\Omega} b(x)u^2. \end{equation} We introduce the finite-dimensional restriction $F_p^\lambda:\,\mathcal{O}_\varepsilon\to\mathbb{R}$ given by \begin{equation}\label{5.3} F_p^\lambda(\xi)=J^\lambda_p(U_\xi+\phi_\xi), \end{equation} where $\phi_\xi$ is the unique solution to problem \eqref{4.1} given by Proposition \ref{prop3.3}. \begin{proposition} \label{prop4.1} The function $F_p^\lambda: \mathcal{O}_\varepsilon\to\mathbb{R}$ is of class $C^1$. Moreover, for all sufficiently large $p$ and $\lambda$, if $D_{\xi}F_p^\lambda(\xi)=0$, then $\xi$ satisfies \eqref{5.1}. \end{proposition} \begin{proof} The function $F_p^\lambda$ is of class $C^1$ since $\xi\to\phi_\xi$ is a $C^1$-map into $H^1(\Omega)$. Then $D_{\xi}F_p^\lambda(\xi)=0$ is equivalent to \begin{equation}\label{5.4} \begin{aligned} 0&=(DJ^\lambda_p)'(U_\xi+\phi_\xi) (D_\xi U_\xi+D_\xi\phi_\xi) \\ &=-\sum_{i=1}^2\sum_{j=1}^mc_{ij}(\xi)\int_{\Omega} e^{U_{\delta_j,\xi_j}}Z_{ij}D_{\xi}U_\xi +\sum_{i=1}^2\sum_{j=1}^mc_{ij}(\xi)\int_{\Omega}D_\xi \big(e^{U_{\delta_j,\xi_j}}Z_{ij}\big) \phi_\xi, \end{aligned} \end{equation} where the seconde equality is due to $\int_{\Omega}e^{U_{\delta_j,\xi_j}}Z_{ij}\phi_\xi=0$. From the definition of $U_\xi$ in \eqref{2.12}, we obtain \begin{align*} \partial_{(\xi_k)_l}U_\xi &=\sum_{j=1}^{m}\frac{1}{\gamma\mu_j^{2/(p-1)}\,} \Big\{\partial_{(\xi_k)_l}\Big[U_{\delta_j,\xi_j}(x) +\frac1p\omega_1\big(\frac{x-\xi_j}{\delta_j}\big) +\frac1{p^2}\omega_2\big(\frac{x-\xi_j}{\delta_j}\big) \\ &\quad +\gamma\mu_j^{2/(p-1)}H_j(x)\Big]+O(1) \Big\}. \end{align*} As in the proof of Lemma \ref{lem2.1}, by the maximum principle with Robin boundary condition we can prove that \begin{align*} \partial_{(\xi_k)_l}\big[\gamma\mu_j^{2/(p-1)}H_j(x)\big] &= \delta_{kj}\Big( 1-\frac{C_1}{4p}-\frac{C_2}{4p^2} \Big)8\pi\partial_{(\xi_k)_l} H_\lambda(x,\xi_j) \\ &\quad -\Big(2-\frac{C_1}{p}-\frac{C_2}{p^2} \Big)\partial_{(\xi_k)_l} \log\mu_j+O(\frac{\rho}{ p}), \end{align*} where $\delta_{kj}$ denote the Kronecker's symbol. Thus, by \eqref{2.1}, \eqref{2.4}, \eqref{2.8}, \eqref{3.2} and \eqref{3.20} we have that \begin{align*} \partial_{(\xi_k)_l}U_\xi &=\frac{1}{\delta_k\gamma\mu_k^{2/(p-1)}} \Big\{\Big( 1-\frac{C_1}{4p}-\frac{C_2}{4p^2}\Big)P Z_{lk} +O\Big(\rho^3 + \frac1p\frac{\delta_k^2}{|x-\xi_k|^2+\delta_k^2}\Big) \Big\}\\ &\quad +O( \frac{1}{\gamma}). \end{align*} On the other hand, it can be shown that $\|D_{\xi}(e^{U_{\delta_j,\xi_j}}Z_{ij})\|_{L^{\infty}(\Omega)}=O(1/\delta_j)$ by computing directly. Consequently, \eqref{5.4} can be written as, for each $l=1,2$ and $k=1,\ldots,m$, \begin{align*} &-\sum_{i,\,j}\frac{c_{ij}(\xi)\big[1+O(\frac1p)\big]}{\delta_k\gamma\mu_k^{2/(p-1)}} \int_{\Omega}e^{U_{\delta_j,\xi_j}}Z_{ij}PZ_{lk} \\ &+\sum_{i,\,j}|c_{ij}(\xi)|O\Big(\frac{1}{\gamma}+\|\phi_\xi\|_{\infty} \int_{\Omega}\big| \partial_{(\xi_k)_l}\big(e^{U_{\delta_j,\xi_j}}Z_{ij}\big)\big|\Big) =0, \end{align*} so that, using \eqref{3.23} and \eqref{4.2}, $$ -\frac{64c_{lk}(\xi)}{\delta_k\gamma\mu_k^{2/(p-1)}} \int_{\mathbb{R}^2}\frac{|y|^2}{(1+|y|^2)^4}dy +O\Big(\frac{1}{\delta_kp\gamma }+ \frac{1}{\gamma}+ \frac{1}{p^3\delta_k} \Big) \sum_{i=1}^2\sum_{j=1}^m\big|c_{ij}(\xi)\big|=0, $$ which implies $c_{lk}(\xi)=0$. \end{proof} \begin{proposition} \label{prop4.2} Let $\varepsilon>0$ be fixed. There exist $p_0>0$ and $\lambda_0>0$ such that for any $p>p_0$ and $\lambda>\lambda_0$, \begin{equation}\label{5.7} \begin{aligned} F_p^\lambda(\xi) &=\frac{4\pi mp}{\gamma^2} -\frac{32\pi^2}{\gamma^2}\varphi_m^\lambda(\xi_1,\ldots,\xi_m) +\frac{4\pi m}{\gamma^2} \\ &\quad +\frac{m}{2\gamma^2}\int_{\mathbb{R}^2} \Big(\frac{8}{(1+|y|^2)^2}U_{1,0}-\Delta\omega_1\Big) +o\big(\frac{1}{p^2}\big) \end{aligned} \end{equation} $C^1$-uniformly with respect to $\xi=(\xi_1,\ldots,\xi_m)\in\mathcal{O}_\varepsilon$, where $\varphi^\lambda_m(\xi)$ is defined by \eqref{5.8}. \end{proposition} \begin{proof} According to the proof of \cite[Proposition 5.3]{EMP}, it suffices to establish the expansion \eqref{5.7} in the $C^0$-sense. Multiplying the first equation in \eqref{4.1} by $U_\xi+\phi_\xi$ and integrating by parts, we obtain \begin{align*} \int_{\Omega}(U_\xi+\phi_\xi)^{p+1} &= \int_{\Omega}|\nabla(U_\xi+\phi_\xi)|^2+ \lambda\int_{\partial\Omega} b(x)(U_\xi+\phi_\xi)^{2} \\ &\quad +\sum_{i=1}^2\sum_{j=1}^m c_{ij}(\xi)\int_{\Omega}e^{U_{\delta_j,\xi_j}} Z_{ij}U_\xi. \end{align*} Since $U_\xi$ is a bounded function, by \eqref{4.2} we obtain $$ \int_{\Omega}(U_\xi+\phi_\xi)^{p+1}= \int_{\Omega}|\nabla(U_\xi+\phi_\xi)|^2+ \lambda\int_{\partial\Omega} b(x)(U_\xi+\phi_\xi)^{2} +O(\frac{1}{p^4}) $$ uniformly for any set of points $\xi=(\xi_1,\ldots,\xi_m)\in\mathcal{O}_\varepsilon$. Hence, by \eqref{5.2}-\eqref{5.3} we obtain \begin{equation}\label{5.9} \begin{aligned} F_p^\lambda(\xi) &=\big(\frac1{2}-\frac{1}{p+1}\big) \Big[\Big(\int_{\Omega}|\nabla U_\xi|^2+ \lambda\int_{\partial\Omega} b(x)U_\xi^2\Big)\\ &\quad +2\Big(\int_{\Omega}\nabla U_\xi\nabla\phi_\xi +\lambda\int_{\partial\Omega} b(x)U_\xi\phi_\xi\Big) \\ &\quad +\Big(\int_{\Omega}|\nabla\phi_\xi|^2+ \lambda\int_{\partial\Omega} b(x)\phi_\xi^2\Big) \Big]+O\big(\frac1{p^4}\big). \end{aligned} \end{equation} We expand the term $\int_{\Omega}|\nabla U_\xi|^2+ \lambda\int_{\partial\Omega}U_\xi^2$: in view of \eqref{2.21} we deduce \begin{align*} &\int_{\Omega}|\nabla U_\xi|^2+ \lambda\int_{\partial\Omega} b(x)U_\xi^2=\int_{\Omega}(-\Delta U_\xi)U_\xi\\ &=\sum_{j=1}^m\frac{1}{\gamma\mu_j^{2/(p-1)}} \int_{B(\xi_j,\varepsilon)} \Big[e^{U_{\delta_j,\xi_j}} -\frac{1}{p\delta_j^2}\Delta\omega_1\big(\frac{x-\xi_j}{\delta_j}\Big) -\frac{1}{p^2\delta_j^2}\Delta\omega_2\big(\frac{x-\xi_j}{\delta_j}\big) \Big]U_\xi \\ &\quad +O(\rho^2)\\ &=\sum_{j=1}^m\frac{1}{\gamma^2\mu_j^{4/(p-1)}} \int_{B(0,\varepsilon/\delta_j)} \Big[\frac{8}{(1+|y|^2)^2} -\frac{1}{p}\Delta\omega_1(y) -\frac{1}{p^2}\Delta\omega_2(y)\Big]\\ &\quad \times\Big[ p+U_{1,0}(y)+\frac1p\omega_1(y)+\frac1{p^2}\omega_2(y)+O(\rho|y|) +O\big(\frac{\rho}p\big)\Big]dy+O\big(\rho^2\big)\\ &=\sum_{j=1}^m\frac{1}{\gamma^2\mu_j^{4/(p-1)}}\Big[8\pi p + \int_{\mathbb{R}^2} \Big(\frac{8}{(1+|y|^2)^2}U_{1,0} -\Delta\omega_1\Big)+O(\frac1p) \Big]\\ &=\frac{8\pi mp}{\gamma^2}-\frac{32\pi}{\gamma^2} \sum_{j=1}^m\log\mu_j +\frac{m}{\gamma^2}\int_{\mathbb{R}^2} \Big(\frac{8}{(1+|y|^2)^2}U_{1,0}-\Delta\omega_1\Big)+O\big(\frac1{p^3}\big) \end{align*} since $\mu_j^{-\frac4{p-1}}=1-\frac4p\log\mu_j+O(\frac{1}{p^2})$. Recalling the expansion \eqref{2.17} of $\mu_j$, then we obtain \begin{equation}\label{5.10} \begin{aligned} \int_{\Omega}|\nabla U_\xi|^2+ \lambda\int_{\partial\Omega} b(x)U_\xi^2 &= \frac{8\pi mp}{\gamma^2} -\frac{64\pi^2}{\gamma^2}\varphi_m^\lambda(\xi) +\frac{24\pi m}{\gamma^2} \\ &\quad +\frac{m}{\gamma^2}\int_{\mathbb{R}^2} \Big(\frac{8}{(1+|y|^2)^2}U_{1,0}-\Delta\omega_1\Big) +O\big(\frac{1}{p^3}\big). \end{aligned} \end{equation} On the other hand, by \eqref{4.2}, we have \begin{equation}\label{5.11} \begin{aligned} &2\Big(\int_{\Omega}\nabla U_\xi\nabla\phi_\xi +\lambda\int_{\partial\Omega} b(x)U_\xi\phi_\xi\Big) +\Big(\int_{\Omega}|\nabla\phi_\xi|^2+ \lambda\int_{\partial\Omega} b(x)\phi_\xi^2\Big)\\ &=O\big(\frac1{p^{7/2}}\big). \end{aligned} \end{equation} Consequently, inserting \eqref{5.10}-\eqref{5.11} in \eqref{5.9}, we obtain \eqref{5.7}. \end{proof} \section{Proof of Theorem \ref{thm1.1}} \begin{definition} \label{def5.1}\rm Let $\mathcal {D}$ be an open set compactly contained in $\Omega^m$ with smooth boundary. We recall that $\varphi_m^\infty$ links in $\mathcal {D}$ at critical level $\mathcal {C}$ relative to $B$ and $B_0$ if $B$ and $B_0$ are closed subsets of $\overline{\mathcal {D}}$ with $B$ connected and $B_0\subset B$ such that the following conditions hold: let us set $\Gamma$ to be the class of all maps $\Phi\in C(B,\mathcal {D})$ with the property that there exists a function $\Psi\in C([0,1]\times B, \mathcal {D})$ such that $$ \Psi(0,\cdot)=Id_B, \quad \Psi(1,\cdot)=\Phi \quad \Psi(t,\cdot)|_{B_0}=Id_{B_0}\quad \text{for all } t\in[0,1]. $$ \end{definition} We assume \begin{equation}\label{6.1} \sup_{y\in B_0}\varphi_m^\infty(y)<\mathcal {C}\equiv\inf_{\Phi\in\Gamma}\sup_{y\in B}\varphi_m^\infty(\Phi(y)), \end{equation} and for all $y\in\partial\mathcal {D}$ such that $\varphi_m^\infty(y)=\mathcal {C}$, there exists a vector $\tau_y$ tangent to $\partial\mathcal {D}$ at $y$ such that \begin{equation}\label{6.2} \nabla\varphi_m^\infty(y)\cdot\tau_y\neq0. \end{equation} Under these conditions a critical point $\bar{y}\in\mathcal {D}$ of $\varphi_m^\infty$ with $\varphi_m^\infty(\bar{y})=\mathcal {C}$ exists, as a standard deformation argument involving the negative gradient flow of $\varphi_m^\infty$ shows. It is easy to check that the above conditions hold if $$ \inf_{x\in\mathcal {D}}\varphi_m^\infty(y)<\inf_{x\in\partial\mathcal {D}}\varphi_m^\infty(x), \quad\text{or}\quad \sup_{x\in\mathcal {D}}\varphi_m^\infty(x)>\sup_{x\in\partial\mathcal {D}}\varphi_m^\infty(x), $$ namely the case of (possibly degenerate) local minimum or maximum points of $\varphi_m^\infty$. We call $\mathcal{C}$ a nontrivial critical level of $\varphi_m^\infty$ in $\mathcal {D}$. \begin{proof}[Proof of Theorem \ref{thm1.1}] Since $\Omega$ is not simply connected, from the proof of \cite[Theorem 1]{DKM1} it follows that given any $m\geq1$, $\varphi_m^\infty$ has a nontrivial critical level $\mathcal{C}$ in some open set $\mathcal {D}$, compactly contained in $\Omega^m$. Let us consider the set $\mathcal {D}$, the associated critical value $\mathcal{C}$ and $\xi\in\mathcal {D}$. According to Proposition \ref{prop4.1}, the function $u_{p,\lambda}=U_\xi+\phi_\xi$ where $U_\xi$ is defined in \eqref{2.12} and $\phi_\xi$ is the unique solution to problem \eqref{4.1} given by Proposition \ref{prop3.3}, is a solution to problem \eqref{a1} if we adjust $\xi$ so that it is a critical point of $F_p^\lambda(\xi)$ defined by \eqref{5.3}. This is equivalent to finding a critical point of $$ \widetilde{F}_p^\lambda(\xi)=-\frac{\gamma^2}{32\pi^2} \Big[F_p^\lambda(\xi) -\frac{4\pi mp}{\gamma^2} -\frac{4\pi m}{\gamma^2}-\frac{m}{2\gamma^2} \int_{\mathbb{R}^2}\Big(\frac{8}{(1+|y|^2)^2}U_{1,0}-\Delta\omega_1\Big) \Big]. $$ On the other hand, from Proposition \ref{prop4.2}, for $\xi=(\xi_1,\ldots,\xi_m)\in\mathcal {D}\cap\mathcal{O}_\varepsilon$, we have \begin{equation}\label{6.3} \widetilde{F}_p^\lambda(\xi)=\varphi_m^\lambda(\xi) +o(1)\Theta_{p,\lambda}(\xi), \end{equation} where $\Theta_{p,\lambda}$ and $\nabla_\xi\Theta_{p,\lambda}$ are uniformly bounded in the considered region as $p$ and $\lambda$ go to $+\infty$. We claim that \begin{equation}\label{6.4} \varphi_m^\lambda(\xi)=\varphi_m^\infty(\xi) +O\big(1/\lambda\big) \quad\text{uniformly in $C^1(\mathcal{O}_\varepsilon)$ as } \lambda\to+\infty. \end{equation} From the definitions of $\varphi_m^\lambda$ and $\varphi_m^\infty$ it suffices to establish that for any set of points $\xi=(\xi_1,\ldots,\xi_m)\in\mathcal{O}_\varepsilon$, $$ H_\lambda(x,\xi_j) =H_\infty(x,\xi_j)+O\big(1/\lambda\big) $$ in $C(\overline{\Omega})$ and in $C_{\rm loc}^1(\Omega)$ as $\lambda\to+\infty$. Indeed, if we set $h(x)=H_\lambda(x,\xi_j)-H_\infty(x,\xi_j)$, then by \eqref{2.15}, \begin{gather*} -\Delta h=0\quad \text{in }\Omega,\\ \frac{\partial h}{\partial \nu}+\lambda b(x)h= -\frac{\partial G_\infty}{\partial\nu}(x,\xi_j) \quad\text{on } \partial\Omega. \end{gather*} Furthermore, by the maximum principle with Robin boundary condition and the definition of $\mathcal{O}_\varepsilon$ in \eqref{2.3}, we deduce $$ \max_{\overline{\Omega}}\big|h(x)\big| + \max_{\overline{\Omega}}\big|\operatorname{dist}(x,\partial\Omega)\nabla h(x)\big| \leq \frac{C}{\lambda}\big\|\frac{\partial G_\infty}{\partial\nu}(x,\xi_j) \big\|_{L^\infty(\partial\Omega)} =O\big(1/\lambda\big). $$ According to Definition \ref{def5.1}, we have that if $M>\mathcal{C}$, then assumptions \eqref{6.1}, \eqref{6.2} still hold for the function $\min\{M,\varphi_m^\infty(\xi)\}$ as well as for $\min\{M,\varphi_m^\infty(\xi) +o(1)\Theta_{p,\lambda}(\xi)+O(1/\lambda)\}$. By \eqref{6.3}-\eqref{6.4} it follows that the function $\min\{M,\widetilde{F}_p^\lambda(\xi)\}$ satisfies for all $p$ and $\lambda$ large assumptions \eqref{6.1}, \eqref{6.2} in $\mathcal {D}$ and therefore has a critical value $\mathcal{C}_{p,\lambda}0$ such that $|\xi_{i,p,\lambda}-\xi_{j,p,\lambda}|>2\varepsilon$, $\operatorname{dist}(\xi_{i,p,\lambda},\partial\Omega)>2\varepsilon$. This implies $C^1$-closeness of $\widetilde{F}_p^\lambda(\xi)$ and $\varphi_m^\infty$ at this level, hence $\nabla\varphi_m^\infty(\xi_{p,\lambda})\to0$ and thus $\nabla\varphi_m^\lambda(\xi_{p,\lambda})\to0$ as $p\to+\infty$ and $\lambda\to+\infty$. The function $u_{p,\lambda}(x)=U_{\xi_{p,\lambda}}(x)+\phi_{\xi_{p,\lambda}}(x)$ is therefore a solution with the qualitative properties as predicted in Theorem \ref{thm1.1}. \end{proof} \subsection*{Acknowledgments} This research is supported by the Fundamental Research Funds for the Central Universities (No. KYZ201541) and the Natural Science Foundation of Jiangsu Province (No. BK20150651). \begin{thebibliography}{99} \bibitem{AG} A. Adimurthi, M. Grossi; \emph{Asymptotic estimates for a two-dimensional problem with polynomial nonlinearity}, Proc. Amer. Math. Soc. 132 (2004), 1013--1019. \bibitem{BP} S. Baraket, F. Pacard; \emph{Construction of singular limits for a semilinear elliptic equation in dimension $2$}, Calc. Var. Partial Differential Equations 6 (1998), 1--38. \bibitem{CI} D. Chae, O. Imanuvilov; \emph{The existence of non-topological multivortex solutions in the relativistic self-dual Chern-Simons theory}, Comm. Math. Phys. 215 (2000), 119--142. \bibitem{CL} W. Chen, C. Li; \emph{Classification of solutions of some nonlinear elliptic equations}, Duke Math. J. 63 (1991), 615--622. \bibitem{D} D. Daners; \emph{Robin boundary value problems on arbitrary domains}, Trans. Amer. Math. Soc. 352 (2000) 4207--4236. \bibitem{DKM} J. D\'{a}vila, M. Kowalczyk, M. Montenegro; \emph{Critical points of the regular part of the harmonic Green function with Robin boundary condition}, J. Funct. Anal. 255 (2008) 1057--1101. \bibitem{DKM1} M. del Pino, M. Kowalczyk, M. Musso; \emph{Singular limits in Liouville-type equations}, Calc. Var. Partial Differential Equations 24 (2005), 47--81. \bibitem{DMO} R. Dillon, P. K. Maini, H. G. Othmer; \emph{Pattern formation in generalized Turing systems, I. Steady-state patterns in systems with mixed boundary conditions}, J. Math. Biol. 32 (1994), 345--393. \bibitem{EG} K. El Mehdi, M. Grossi; \emph{ Asymptotic estimates and qualitative properties of an elliptic problem in dimension two}, Adv. Nonlinear Stud. 4 (2004), 15--36. \bibitem{EGP} P. Esposito, M. Grossi, A. Pistoia; \emph{On the existence of blowing-up solutions for a mean field equation}, Ann. I. H. Poincar\'{e} Anal. Non Lin\'{e}aire 22 (2005), 227--257. \bibitem{EMP} P. Esposito, M. Musso, A. Pistoia; \emph{Concentrating solutions for a planar elliptic problem involving nonlinearities with large exponent}, J. Differential Equations 227 (2006), 29--68. \bibitem{GT1} M. Grossi, F. Takahashi; \emph{Nonexistence of multi-bubbles solutions to some elliptic equations on convex domains}, J. Funct. Anal. 259 (2010) 904--917. \bibitem{GT} D. Gilbarg, N. S. Trudinger; \emph{Elliptic Partial Differential Equations of Second Order}, Springer-Verlag, Berlin, 2001. \bibitem{MW1} M. Musso, J. Wei; \emph{Stationary solutions to a Keller-Segel chemotaxis system}, Asymptot. Anal. 49 (2006) 217--247. \bibitem{N} W.-M. Ni; \emph{Diffusion, cross-diffusion, and their spike-layer steady states}, Notices Amer. Math. Soc. 45 (1998) 9--18. \bibitem{RW1} X. Ren, J. Wei; \emph{On a two-dimensional elliptic problem with large exponent in nonlinearity}, Trans. Amer. Math. Soc. 343 (1994) 749--763. \bibitem{RW2} X. Ren, J. Wei; \emph{Single-point condensation and least-energy solutions}, Proc. Amer. Math. Soc. 124 (1996) 111--120. \bibitem{ZS} Y. Zhang, L. Shi; \emph{Concentrating solutions for a planar elliptic problem with large nonlinear exponent and Robin boundary condition}, Advances in Nonlinear Analysis, online, 2016. \end{thebibliography} \end{document}