#  Boundary clustered layer positive solutions for an elliptic Neumann problem
with large exponent

Yibin Zhang  College of Sciences, Nanjing Agricultural University, Nanjing
210095, China<EMAIL_ADDRESS>
###### Abstract.

Let $\mathcal{D}$ be a smooth bounded domain in $\mathbb{R}^{N}$ with $N\geq
3$, we study the existence and profile of positive solutions for the following
elliptic Nenumann problem

$\left\\{\begin{aligned} &-\Delta\upsilon+\upsilon=\upsilon^{p},\,\,\,\
\,\upsilon>0\,\,\,\ \,\textrm{in}\,\,\,\,\,\,\mathcal{D},\\\\[2.84526pt]
&\frac{\partial\upsilon}{\partial\nu}=0\,\,\,\,\,\ \,\ \,\qquad\,\ \,\,\,\,\
\,\,\quad\quad\,\,\
\textrm{on}\,\,\,\,\partial\mathcal{D},\end{aligned}\right.$

where $p>1$ is a large exponent and $\nu$ denotes the outer unit normal vector
to the boundary $\partial\mathcal{D}$. For suitable domains $\mathcal{D}$, by
a constructive way we prove that, for any integers $l$, $m$ with $0\leq l\leq
m$ and $m\geq 1$, if $p$ is large enough, such a problem has a family of
positive solutions with $l$ interior layers and $m-l$ boundary layers which
concentrate along $m$ distinct $(N-2)$-dimensional minimal submanifolds of
$\partial\mathcal{D}$, or collapse to the same $(N-2)$-dimensional minimal
submanifold of $\partial\mathcal{D}$ as $p\rightarrow+\infty$.

###### Key words and phrases:

Boundary clustered layer positive solutions; Large exponent; Anisotropic
Green’s function.

###### 2010 Mathematics Subject Classification:

Primary 35B25; Secondary 35B38, 35J25.

##  1\. Introduction

In this paper, we are interested in the following classical elliptic Neumann
problem

$\left\\{\begin{aligned} &-d^{2}\Delta\upsilon+\upsilon=\upsilon^{p},\,\,\,\
\,\upsilon>0\,\,\,\ \,\textrm{in}\,\,\,\,\,\,\mathcal{D},\\\\[2.84526pt]
&\frac{\partial\upsilon}{\partial\nu}=0\,\,\,\,\,\ \,\
\,\quad\qquad\,\,\,\,\,\,\ \,\,\quad\quad\,\,\
\textrm{on}\,\,\,\,\partial\mathcal{D},\end{aligned}\right.$ (1.1)

where $\mathcal{D}$ is a smooth bounded domain in $\mathbb{R}^{N}$, $N\geq 2$,
$d>0$ is a parameter, $p>1$ is a exponent and $\nu$ denotes the outer unit
normal vector to the boundary $\partial\mathcal{D}$.

Problem (1.1) has received considerable attention in the last three decades,
because it appears in many different mathematical models: for example, it
arises from the study of steady states for the logarithmic Keller-Segel system
in chemotaxis [32] and the shadow system of the Gierer-Meinhardt model in
biological pattern formation [22]. In particular, it has been shown that the
solutions of (1.1) exhibit a variety of interesting concentration phenomena as
either the exponent $p$ tends to some critical values or the parameter $d$
tends to zero.

Let us first define $p_{h+1}^{*}:=(N-h+2)/(N-h-2)$, $0\leq h\leq N-3$ as the
$(h+1)$-th critical exponent and set $p_{N-1}^{*}:=+\infty$. In the
subcritical case, i.e. $N=2$, or $N\geq 3$ and $p<p_{1}^{*}=(N+2)/(N-2)$,
compactness of Sobolev’s embedding ensures the existence of a positive least
energy solution of (1.1). For $d>0$ sufficiently small, Lin, Ni and Takagi in
[33, 44, 45] proved that this least energy solution has exactly one bounded,
very sharp spike located on the boundary and near the most curved part of the
boundary, i.e., the region where the mean curvature of the boundary attains
its maximum. Higher energy solutions of (1.1) with multiple boundary peaks as
well as multiple interior peaks have been established in [11, 12, 13, 19, 24,
28, 29, 30, 31, 35, 54, 56]. It turns out that multiple boundary spikes tend
to cluster around critical points of the mean curvature of the boundary, while
the location of multiple interior spikes is determined by the distance between
the peaks and the boundary. In particular, Gui and Wei in [29] proved that for
any integers $l$, $m$ with $0\leq l\leq m$ and $m\geq 1$, problem (1.1) has a
solution with exactly $l$ interior spikes and $m-l$ boundary spikes provided
that $d$ is small and $p$ is subcritical. Generally, such spiky solutions are
called solutions with $0$-dimensional concentration sets.

In the critical case, i.e. $p=p_{1}^{*}$, such type of concentration phenomena
occurs, but the situation is quite different. The lack of compactness of
Sobolev’s embedding makes it non-obvious to apply variational techniques to
obtain a nonconstant least energy solution of (1.1) when $d$ is sufficiently
small. The first existence result for a nonconstant least energy solution of
(1.1) in general domains and $d$ small was obtained by Adimurthi-Mancini [1]
and Wang [53]. The profile and asymptotic behavior of this least energy
solution has been clarified in the subsequent works [4, 43, 47]. As in the
subcritical case, the least energy solution has a unique maximum point or peak
that lies on the boundary and goes, as $d$ tends to zero, to a maximum point
of the mean curvature of the boundary. Unlike the subcritical case, the least
energy solution blows up, as $d$ tends to zero, at a maximum point of the mean
curvature of the boundary. Higher energy solutions with one or more separated
or clustered boundary blow-up points have been exhibited for instance in [2,
3, 23, 26, 34, 47, 48, 57], where these blow-up points are nothing but
degenerate or non-degenerate critical points of the mean curvature of the
boundary. A major difference with the subcritical case is that the condition
of positivity for the mean curvature at these critical points turns out to be
necessary for the boundary bubbling phenomenon to take place [27, 47]. Note
that in striking contrast with the subcritical case, there are no solutions
blowing up at only interior points when $d$ tends to zero, namely at least one
blow-up point has to lie on the boundary, even all of the blow-up points for
solutions with uniformly bounded energy have to lie only on the boundary as
established in [9, 20, 49]. Moreover, for the subcritical case, spiky
solutions with only interior peaks always exist, but these interior peaks must
stay in the domain as $d$ tends to zero [55]; while for the critical case,
interior peaks of solutions could not remain in the domain as $d$ tends to
zero [23, 25].

The almost first critical case, i.e. $p=p_{1}^{*}\pm\epsilon$ with $\epsilon$
positive and sufficiently small, has been widely studied. In the slightly
subcritical case, i.e. $p=p_{1}^{*}-\epsilon$, the existence of solutions with
simple or non-simple blow-up points located on the boundary near critical
points of the mean curvature of the boundary with negative value as
$\epsilon\rightarrow 0$ for fixed $d$ and $N\geq 4$ was established in [10,
50, 52]; while in the slightly supercritical case, i.e.
$p=p_{1}^{*}+\epsilon$, it was proved in [16, 50, 52] that there also exists a
solution with simple or non-simple blow-up points located on the boundary near
critical points of the mean curvature of the boundary with positive value as
$\epsilon\rightarrow 0$ for fixed $d$ and $N\geq 4$. In [51], it was proved
that if $N=3$ and $\epsilon$ tends to zero from below or above, then a
solution with one interior blow-up point may exist for finite $d$. Moreover,
if $N=2$, $d$ is fixed and the exponent $p$ goes to $+\infty$, Musso and Wei
in [41] proved that for any integers $l$, $m$ with $0\leq l\leq m$ and $m\geq
1$, problem (1.1) admits a solution with exactly $l$ interior bubbles and
$m-l$ boundary bubbles, whose location can be characterized by critical points
of a certain combination of Green’s function and its regular part.

It seems natural to ask if problem (1.1) has solutions that exhibit
concentration phenomena on $h$-dimensional subsets of $\overline{\mathcal{D}}$
for every $1\leq h\leq N-1$, as conjectured by Ni in [42]. In particular,
given an $h$-dimensional submanifold $\Gamma$ of $\partial\mathcal{D}$ and
assuming that either $h\geq N-2$ or $p\leq p_{h+1}^{*}+\epsilon$ with
$\epsilon$ positive and sufficiently small, one question is whether problem
(1.1) admits a solution that concentrates along $\Gamma$ as either $d$ or
$\epsilon$ tends to zero. For results in this direction, we first mention the
$(h+1)$-th subcritical case, i.e. $p<p_{h+1}^{*}$, where Malchiodi and
Montenegro [38, 39] proved that, given any $N\in\mathbb{N}$ and any $p>1$,
there exist solutions concentrating on the whole boundary if the sequence $d$
satisfies some gap condition, corresponding to $h=N-1$. Furthermore, the
result was extended in [36, 37] for general $h$, and the concentration set
$\Gamma$ is an embedded closed minimal submanifold of $\partial\mathcal{D}$
which is in addition nondegenerate  in the sense that its Jacobi operator is
invertible. Indeed, this phenomenon is rather subtle compared with pointwise
concentration: existence can only be achieved along a sequence of parameters
$d=d_{i}\rightarrow 0$. The sequence of parameters $d_{i}$ must be suitable
away from certain values where resonance occurs, and the topological type of
the solution changes: unlike the pointwise concentration, the Morse index of
these solutions is very large and grows as $d_{i}\rightarrow 0$. Del Pino,
Mahmoudi and Musso [15] extended this type of results to the $(h+1)$-th
critical case $p=p_{h+1}^{*}$ and proved that if $\Gamma$ is an embedded
closed minimal submanifold of $\partial\mathcal{D}$ with dimension $h\leq N-7$
(in particular, $N\geq 8$) which is nondegenerate, and a certain weighted
average of sectional curvature of $\partial\mathcal{D}$ is positive along
$\Gamma$, then problem (1.1) admits a solution for a suitable sequence of
parameters $d_{i}\rightarrow 0$ which blows up along $\Gamma$. Further for the
almost $(h+1)$-th critical case $p=p_{h+1}^{*}\pm\epsilon$ with
$0<\epsilon\rightarrow 0$ but $d$ fixed, the same conclusion may be
established by Deng, Mahmoudi and Musso [18] under the analogous assumption.
Meanwhile, for any integer $h=1,\ldots,N-3$, Manna and Pistoia [40] proved
that in some suitable domains $\mathcal{D}$, problem (1.1) has a solution
which blows up along an $h$-dimensional minimal submanifold $\Gamma$ of
$\partial\mathcal{D}$ as $p$ approaches from either below or above the
$(h+1)$th critical exponent $p_{h+1}^{*}$.

In the present paper we consider the almost $(N-1)$th critical case, i.e.
$h=N-2$, and give a positive answer when $d$ is fixed and $p$ goes to
$+\infty$. More precisely, we find some domains $\mathcal{D}$ such that if
$d=1$ and $p$ is large enough, then for any positive integer $m$, problem
(1.1) has a positive solution with $m$ distinct mixed interior and boundary
layers which concentrate along $m$ distinct $(N-2)$-dimensional minimal
submanifolds of $\partial\mathcal{D}$, or collapse to the same
$(N-2)$-dimensional minimal submanifold of $\partial\mathcal{D}$ as $p$ goes
to $+\infty$.

Let $\Omega$ be a smooth bounded domain $\Omega$ in $\mathbb{R}^{2}$ such that

$\displaystyle\overline{\Omega}\subset\\{(x_{1},x_{2})\in\mathbb{R}^{2}|\,\,x_{1}>0\\}\,\,\qquad\,\,\textrm{or}\,\,\,\qquad\,\overline{\Omega}\subset\\{(x_{1},x_{2})\in\mathbb{R}^{2}|\,\,x_{1},\,x_{2}>0\\}.$
(1.2)

Let $n=1$ or $n=2$ be fixed. Fix $k_{1},k_{n}\in\mathbb{N}$ with
$h:=k_{1}+k_{n}=N-2$ and set

$\displaystyle\mathcal{D}:=\big{\\{}(y_{1},y_{n},x^{\prime})\in\mathbb{R}^{k_{1}+1}\times\mathbb{R}^{k_{n}+1}\times\mathbb{R}^{2-n}|\,\,(|y_{1}|,|y_{n}|,x^{\prime})\in\Omega\big{\\}}.$
(1.3)

Then $\mathcal{D}$ is a smooth bounded domain $\mathbb{R}^{N}$ which is
$\Upsilon$-invariant for the action of the group
$\Upsilon:=\mathcal{O}(k_{1}+1)\times\mathcal{O}(k_{n}+1)$ on $R^{N}$ given by

$\displaystyle(g_{1},g_{n})(y_{1},y_{n},x^{\prime}):=(g_{1}y_{1},g_{n}y_{n},x^{\prime}).$

Here $\mathcal{O}(k_{i}+1)$ denotes the group of linear isometries of
$\mathbb{R}^{k_{i}+1}$. For $p>1$ large enough we shall look for
$\Upsilon$-invariant solutions of problem (1.1) with $d=1$, i.e. solutions
$\upsilon$ of the form

$\displaystyle\upsilon(y_{1},y_{n},x^{\prime})=u(|y_{1}|,|y_{n}|,x^{\prime}).$
(1.4)

Then a simple calculation shows that $\upsilon$ solves problem (1.1) with
$d=1$ if and only if $u$ satisfies

$\left\\{\begin{aligned} &-\Delta
u-\sum_{i=1}^{n}\frac{k_{i}}{\,x_{i}\,}\frac{\partial u}{\partial
x_{i}}+u=u^{p},\,\,\,\,\,u>0\,\,\,\,\,\,\textrm{in}\,\,\,\,\,\Omega,\\\
&\frac{\partial u}{\partial\nu}=0\,\,\,\,\,\ \,\
\,\qquad\qquad\qquad\qquad\,\,\,\,\,\ \,\,\quad\quad\,\,\
\textrm{on}\,\,\,\,\partial\Omega.\end{aligned}\right.$ (1.5)

Thus, we are led to study the more general anisotropic problem

$\left\\{\begin{aligned} &-\nabla(a(x)\nabla
u)+a(x)u=a(x)u^{p},\,\,\,\,u>0\,\,\,\,\,\textrm{in}\,\,\,\,\,\Omega,\\\\[2.84526pt]
&\frac{\partial
u}{\partial\nu}=0\,\,\qquad\quad\qquad\qquad\qquad\qquad\qquad\ \ \ \
\,\,\textrm{on}\,\,\,\partial\Omega,\end{aligned}\right.$ (1.6)

where $\Omega$ is a smooth bounded domain in $\mathbb{R}^{2}$, $a(x)$ is a
positive smooth function over $\overline{\Omega}$, $p>1$ is a large exponent
and $\nu$ denotes the outer unit normal vector to $\partial\Omega$. Note that
if

$\displaystyle a(x)=a(x_{1},x_{n},x^{\prime}):=x_{1}^{k_{1}}x_{n}^{k_{n}},$
(1.7)

then problem (1.6) can be rewritten as equation (1.5).

Our goal is to construct solutions $u_{p}$ to problem (1.6) with $m$ distinct
mixed interior and boundary spikes which concentrate at points
$\xi^{*}_{1},\ldots,\xi^{*}_{m}$ of $\partial\Omega$, or accumulate to the
same point $\xi_{*}$ of $\partial\Omega$ as $p$ goes to $+\infty$. They
correspond, via (1.4), to $\Upsilon$-invariant solutions $\upsilon_{p}$ of
problem (1.1)$\large|_{d=1}$ with $m$ distinct mixed interior and boundary
layers which concentrate along the $\Upsilon$-orbits
$\Xi(\xi^{*}_{1}),\ldots,\Xi(\xi^{*}_{m})$ of $\partial\mathcal{D}$, or
collapse to the same $\Upsilon$-orbit $\Xi(\xi_{*})$ of $\partial\mathcal{D}$
as $p$ goes to $+\infty$. Here

$\displaystyle\Xi(\zeta):=\big{\\{}(y_{1},y_{n},x^{\prime})\in\partial\mathcal{D}|\,\,(|y_{1}|,|y_{n}|,x^{\prime})=\zeta\in\partial\Omega\big{\\}}$

is a $(N-2)$-dimensional minimal submanifold of $\partial\mathcal{D}$
diffeomorphic to $\mathbb{S}^{k_{1}}\times\mathbb{S}^{k_{n}}$ (note that
$k_{1}+k_{n}=h=N-2$), where $\mathbb{S}^{k_{i}}$ is the unit sphere in
$\mathbb{R}^{k_{i}+1}$.

Let us define the linear differential operator

$\displaystyle\Delta_{a}u=\frac{1}{a(x)}\nabla(a(x)\nabla u)=\Delta
u+\nabla\log a(x)\nabla u$

and the Green’s function associated with the Neumann problem

$\left\\{\begin{aligned}
&-\Delta_{a}G(x,y)+G(x,y)=\delta_{y}(x),\,\,\,\,\,\,\,x\in\Omega,\\\
&\frac{\partial
G}{\partial\nu_{x}}(x,y)=0,\qquad\,\,\qquad\qquad\qquad\,\,x\in\partial\Omega,\end{aligned}\right.$
(1.8)

for every $y\in\overline{\Omega}$. The regular part of $G(x,y)$ is defined
depending on whether $y$ lies in the domain or on its boundary as

$\displaystyle H(x,y)=\left\\{\begin{aligned}
&G(x,y)+\frac{1}{2\pi}\log|x-y|,\,\quad\,y\in\Omega,\\\\[2.84526pt]
&G(x,y)+\frac{1}{\pi}\log|x-y|,\,\
\quad\,y\in\partial\Omega.\end{aligned}\right.$ (1.9)

Our first result concerns the existence of solutions of problem (1.6) whose
interior and boundary spikes are uniformly far away from each other and
interior spikes lie in the domain with distance to the boundary uniformly
approaching zero.

Theorem 1.1. Let $l\in\\{0,\ldots,m\\}$, $m\geq 1$ and assume that there exist
$m$ different points $\xi^{*}_{1},\ldots,\xi^{*}_{m}\in\partial\Omega$ such
that each point $\xi_{i}^{*}$ is either a strict local maximum or a strict
local minimum point of $a(x)$ on $\partial\Omega$ and satisfies that for all
$i=1,\ldots,l$, $\partial_{\nu}a(\xi_{i}^{*}):=\langle\nabla
a(\xi_{i}^{*}),\,\nu(\xi_{i}^{*})\rangle>0$. Then, there exists $p_{m}$ such
that for any $p>p_{m}$, there is a family of positive solutions $u_{p}$ for
problem (1.6) with $m-l$ different boundary spikes and $l$ different interior
spikes located at distance $O\left(1/p\right)$ from $\partial\Omega$. More
precisely,

$\displaystyle
u_{p}(x)=\sum\limits_{i=1}^{m}\frac{1}{\gamma\mu_{i}^{2/(p-1)}}\left[\,\log\frac{1}{(\delta_{i}^{2}+|x-\xi^{p}_{i}|^{2})^{2}}+c_{i}H(x,\xi^{p}_{i})\right]+o\left(1\right),$

where $o(1)\rightarrow 0$, as $p\rightarrow+\infty$, on each compact subset of
$\overline{\Omega}\setminus\\{\xi^{p}_{1},\ldots,\xi^{p}_{m}\\}$, the
parameters $\gamma$, $\delta_{i}$ and $\mu_{i}$ satisfy

$\displaystyle\gamma=p^{p/(p-1)}\varepsilon^{2/(p-1)},\,\,\quad\quad\quad\,\,\delta_{i}=\mu_{i}\varepsilon,\quad\quad\quad\,\,\varepsilon=e^{-p/4},\quad\quad\quad\,\,\frac{1}{C}<\mu_{i}<Cp^{2(m^{2}+1)},$

for some $C>0$,
$(\xi^{p}_{1},\ldots,\xi^{p}_{m})\in\Omega^{l}\times(\partial\Omega)^{m-l}$
satisfies

$\displaystyle\xi^{p}_{i}\rightarrow\xi^{*}_{i}\,\quad\,\textrm{for
all}\,\,\,i,\,\qquad\quad\,\textrm{and}\,\qquad\quad\,\text{dist}(\xi^{p}_{i},\partial\Omega)=O\left(1/p\right)\quad\,\,\forall\,\,i=1,\ldots,l,$

and $c_{i}=8\pi$ for $i=1,\ldots,l$, but $c_{i}=4\pi$ for $i=l+1,\ldots,m$. In
particular, for any $d>0$, as $p\rightarrow+\infty$,

$\displaystyle pu_{p}^{p+1}\rightharpoonup 8\pi
e\sum_{i=1}^{l}\delta_{\xi^{*}_{i}}+4\pi
e\sum_{i=l+1}^{m}\delta_{\xi^{*}_{i}}\,\quad\textrm{weakly in the sense of
measure in}\,\,\,\,\overline{\Omega},$ $\displaystyle u_{p}\rightarrow
0\,\quad\textrm{uniformly
in}\,\,\,\,\overline{\Omega}\setminus\bigcup_{i=1}^{m}B_{d}(\xi^{*}_{i}),$

and

$\displaystyle\sup\limits_{x\in\overline{\Omega}\cap
B_{d}(\xi^{*}_{i})}u_{p}(x)\rightarrow\sqrt{e}.$

The corresponding result for problem (1.1)$\large|_{d=1}$ can be stated as
follows.

Theorem 1.2. Let $l\in\\{0,\ldots,m\\}$, $m\geq 1$ and $\mathcal{D}$ be as in
(1.3). If the assumption of Theorem 1.1 holds, then, there exists $p_{m}$ such
that for any $p>p_{m}$, problem (1.1)$\large|_{d=1}$ has a positive solution
$\upsilon_{p}$ with $m-l$ boundary layers and $l$ interior layers which
concentrate along $m$ different $(N-2)$-dimensional minimal submanifolds of
$\partial\mathcal{D}$, namely the $\Upsilon$-orbit $\Xi(\xi^{*}_{i})$ of
$\xi^{*}_{i}$ for every $i=1,\ldots,m$, as $p\rightarrow+\infty$.

Our next result concerns the existence of solutions of problem (1.6) with
mixed interior and boundary spikes which accumulate to the same point of the
boundary.

Theorem 1.3. Let $l\in\\{0,\ldots,m\\}$, $m\geq 1$ and assume that
$\xi_{*}\in\partial\Omega$ is a strict local maximum point of $a(x)$ and
satisfies $\partial_{\nu}a(\xi_{*}):=\langle\nabla
a(\xi_{*}),\,\nu(\xi_{*})\rangle=0$. Then, there exists $p_{m}$ such that for
any $p>p_{m}$, there is a family of positive solutions $u_{p}$ for problem
(1.6) with $m-l$ different boundary spikes and $l$ different interior spikes
which accumulate to $\xi_{*}$ as $p\rightarrow+\infty$. More precisely,

$\displaystyle
u_{p}(x)=\sum\limits_{i=1}^{m}\frac{1}{\gamma\mu_{i}^{2/(p-1)}}\left[\,\log\frac{1}{(\delta_{i}^{2}+|x-\xi^{p}_{i}|^{2})^{2}}+c_{i}H(x,\xi^{p}_{i})\right]+o\left(1\right),$

where $o(1)\rightarrow 0$, as $p\rightarrow+\infty$, on each compact subset of
$\overline{\Omega}\setminus\\{\xi^{p}_{1},\ldots,\xi^{p}_{m}\\}$, the
parameters $\gamma$, $\delta_{i}$ and $\mu_{i}$ satisfy

$\displaystyle\gamma=p^{p/(p-1)}\varepsilon^{2/(p-1)},\,\,\quad\quad\quad\,\,\delta_{i}=\mu_{i}\varepsilon,\quad\quad\quad\,\,\varepsilon=e^{-p/4},\quad\quad\quad\,\,\frac{1}{C}<\mu_{i}<Cp^{2(m^{2}+1)},$

for some $C>0$,
$(\xi^{p}_{1},\ldots,\xi^{p}_{m})\in\Omega^{l}\times(\partial\Omega)^{m-l}$
satisfies

$\displaystyle\xi^{p}_{i}\rightarrow\xi_{*}\,\quad\,\textrm{for
all}\,\,\,i,\qquad|\xi^{p}_{i}-\xi^{p}_{k}|>p^{-2(m^{2}+1)}\quad\forall\,\,i\neq
k,\qquad\textrm{and}\qquad\text{dist}(\xi^{p}_{i},\partial\Omega)>p^{-2(m^{2}+1)}\quad\forall\,\,i=1,\ldots,l,$

and $c_{i}=8\pi$ for $i=1,\ldots,l$, but $c_{i}=4\pi$ for $i=l+1,\ldots,m$. In
particular, for any $d>0$, as $p\rightarrow+\infty$,

$\displaystyle pu_{p}^{p+1}\rightharpoonup 4\pi
e(m+l)\delta_{\xi_{*}}\,\quad\textrm{weakly in the sense of measure
in}\,\,\,\,\overline{\Omega},$ $\displaystyle u_{p}\rightarrow
0\,\quad\textrm{uniformly in}\,\,\,\,\overline{\Omega}\setminus
B_{d}(\xi_{*}),$

and

$\displaystyle\sup\limits_{x\in\overline{\Omega}\cap
B_{d/p^{4(m^{2}+1)}}(\xi_{i}^{p})}u_{p}(x)\rightarrow\sqrt{e}.$

The corresponding result for problem (1.1)$\large|_{d=1}$ can be stated as
follows.

Theorem 1.4. Let $m\geq 1$, $l\in\\{0,\ldots,m\\}$ and $\mathcal{D}$ be as in
(1.3). If the assumption of Theorem 1.3 holds, then, there exists $p_{m}$ such
that for any $p>p_{m}$, problem (1.1)$\large|_{d=1}$ has a positive solution
$\upsilon_{p}$ with $m-l$ boundary layers and $l$ interior layers which
collapse to the same $(N-2)$-dimensional minimal submanifold of
$\partial\mathcal{D}$, namely the $\Upsilon$-orbit $\Xi(\xi_{*})$ of
$\xi_{*}$, as $p\rightarrow+\infty$.

Let us remark that the assumption condition in Theorems $1.3$-$1.4$ contains
the following two cases:
(A1) $\xi_{*}\in\partial\Omega$ is a strict local maximum point of $a(x)$
restricted on $\partial\Omega$;
(A2) $\xi_{*}\in\partial\Omega$ is a strict local maximum point of $a(x)$
restricted in $\Omega$ and satisfies $\partial_{\nu}a(\xi_{*})=\langle\nabla
a(\xi_{*}),\,\nu(\xi_{*})\rangle=0$.
In fact, arguing exactly along the sketch of the proof of Theorem $1.3$, we
can easily find that if (A1) holds, then problem (1.6) has positive solutions
with arbitrarily many boundary spikes which accumulate to $\xi_{*}$ along
$\partial\Omega$; while if (A2) holds, then problem (1.6) has positive
solutions with arbitrarily many interior spikes which accumulate to $\xi_{*}$
along the inner normal direction of $\partial\Omega$. For the latter case, it
seems that this paper is the first one in the literature obtaining this type
of concentration phenomenon for positive solutions of some two-dimensional
anisotropic nonlinear elliptic Neumann problems, see [5] as an instance.

The general strategy for proving our main results relies on a very well known
Lyapunov-Schmidt reduction. In Section $2$ we provide an appropriate
approximation for a solution of problem (1.6) and give a basic estimate for
the scaling error term created by the choice of our approximation. Then we
rewrite problem (1.1) 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 the problem of
finding spike solutions of (1.6) to that of finding a critical point of a
finite-dimensional function. Section $6$ concerns with an asymptotic expansion
for the finite-dimensional function appeared in Section $5$. Finally, in
Section $7$ we provide the detailed proof of Theorems $1.1$ and $1.3$.

In this paper, the letter $C$ will always denote a generic positive constant
independent of $p$, which could be changed from one line to another. The
symbol $o(t)$ (respectively $O(t)$) will denote a quantity for which
$\frac{o(t)}{|t|}$ tends to zero (respectively, $\frac{O(t)}{|t|}$ stays
bounded ) as parameter $t$ goes to zero. Moreover, we will use the notation
$o(1)$ (respectively $O(1)$) to stand for a quantity which tends to zero
(respectively, which remains uniformly bounded) as $p\rightarrow+\infty$.

##  2\. An approximation for the solution

In this section we provide an appropriate approximation for a solution of
problem (1.6) and give a basic estimate for the scaling error term created by
the choice of our approximation. Since the function $H(x,y)$ defined in (1.9)
plays an essential role in our construction, we shall first state its
asymptotic behavior without proof, see [5] for details.

Consider the vector function $T(x)=(T_{1}(x),T_{2}(x))$ as the solution of

$\displaystyle\Delta_{x}T-T=\frac{x}{|x|^{2}}\,\,\quad\,\,\textrm{in}\,\,\,\,\mathbb{R}^{2},\,\,\quad\,\,T(x)\in
L^{\infty}_{loc}(\mathbb{R}^{2}).$ (2.1)

Then standard elliptic regularity theory implies that for any $1<q<2$,
$T(x)\in W^{2,q}_{loc}(\mathbb{R}^{2})\cap
C^{\infty}(\mathbb{R}^{2}\setminus\\{0\\})$, and the Sobolev embeddings yield
that $T(x)\in W^{1,1/\alpha}(B_{r}(0))\cap C^{\alpha}(\overline{B_{r}(0)})$
for any $r>0$ and $0<\alpha<1$.

Lemma 2.1. Let $T(x)$ be the function described in (2.1). There exists a
function $H_{1}(x,y)$ such that
(i) for every $x,y\in\overline{\Omega}$,

$\displaystyle H(x,y)=H_{1}(x,y)+\left\\{\begin{aligned}
&\frac{1}{2\pi}\nabla\log a(y)\cdot
T(x-y),\,\quad\,\,y\in\Omega,\\\\[2.84526pt] &\,\frac{1}{\pi}\,\nabla\log
a(y)\cdot T(x-y),\,\,\quad\,y\in\partial\Omega,\end{aligned}\right.$ (2.2)

(ii) the mapping $y\in\overline{\Omega}\mapsto H_{1}(\cdot,y)$ belongs to
$C^{1}\big{(}\Omega,C^{1}(\overline{\Omega})\big{)}\cap
C^{1}\big{(}\partial\Omega,C^{1}(\overline{\Omega})\big{)}$.
In this way, $y\in\overline{\Omega}\mapsto H(\cdot,y)\in
C\big{(}\Omega,C^{\alpha}(\overline{\Omega})\big{)}\cap
C\big{(}\partial\Omega,C^{\alpha}(\overline{\Omega})\big{)}$ and $H(x,y)\in
C^{\alpha}\big{(}\overline{\Omega}\times\Omega\big{)}\cap
C^{\alpha}\big{(}\overline{\Omega}\times\partial\Omega\big{)}\cap
C^{1}\big{(}\overline{\Omega}\times\Omega\setminus\\{x=y\\}\big{)}\cap
C^{1}\big{(}\overline{\Omega}\times\partial\Omega\setminus\\{x=y\\}\big{)}$
for any $\alpha\in(0,1)$, and the corresponding Robin’s function
$y\in\overline{\Omega}\mapsto H(y,y)$ belongs to $C^{1}(\Omega)\cap
C^{1}(\partial\Omega)$.

Let $d>0$ be a sufficiently small but fixed number such that for any
$y\in\Omega$ with $\text{dist}(y,\partial\Omega)<d$, we can define a
reflection of $y$ across $\partial\Omega$ along the outer normal direction,
$y^{*}\in\Omega^{c}$, and get that $|y-y^{*}|=2\text{dist}(y,\partial\Omega)$.
Set

$\displaystyle\Omega_{d}:=\left\\{y\in\Omega\big{|}\,\,\text{dist}(y,\partial\Omega)<d\,\right\\}.$

Lemma 2.2. There exists a mapping $y\in\Omega_{d}\mapsto\mathrm{z}(\cdot,y)\in
C\big{(}\Omega_{d},C^{\alpha}(\overline{\Omega})\big{)}\cap
L^{\infty}\big{(}\Omega_{d},C^{\alpha}(\overline{\Omega})\big{)}$ for any
$\alpha\in(0,1)$ such that for any $x\in\overline{\Omega}$ and
$y\in\Omega_{d}$,

$\displaystyle H(x,y)=\frac{1}{2\pi}\log\frac{1}{|x-y^{*}|}+\mathrm{z}(x,y).$
(2.3)

Even more, for any $x\in\overline{\Omega}$ and $y\in\Omega_{d}$,

$\displaystyle\mathrm{z}(x,y)=\frac{1}{2\pi}\nabla\log a(y)\cdot
T(x-y)-\frac{1}{2\pi}\nabla\log a(y^{*})\cdot
T(x-y^{*})+\tilde{\mathrm{z}}(x,y),$ (2.4)

where the mapping
$y\in\overline{\Omega_{d}}\mapsto\tilde{\mathrm{z}}(\cdot,y)$ belongs to
$C^{1}\big{(}\overline{\Omega_{d}},\,C^{1}(\overline{\Omega})\big{)}$.

Corollary 2.3. Under the assumptions in Lemma 2.2, the Robin’s function
$y\in\Omega\mapsto H(y,y)$ satisfies

$\displaystyle
H(y,y)=\frac{1}{2\pi}\log\frac{1}{|y-y^{*}|}+\mathrm{z}(y),\,\,\quad\forall\,\,y\in\Omega_{d},$
(2.5)

where $\mathrm{z}\in C^{1}\big{(}\overline{\Omega_{d}}\big{)}$ and

$\displaystyle\mathrm{z}(y)=\frac{1}{2\pi}\nabla\log a(y)\cdot
T(0)-\frac{1}{2\pi}\nabla\log a(y^{*})\cdot
T(y-y^{*})+\tilde{\mathrm{z}}(y,y),\,\,\quad\forall\,\,y\in\Omega_{d}.$ (2.6)

The key ingredient to describe the shape of the approximate solution of (1.6)
is based on the standard bubble

$\displaystyle
U_{\delta,\xi}(x)=\log\frac{8\delta^{2}}{(\delta^{2}+|x-\xi|^{2})^{2}},\quad\,\,\delta>0,\,\,\,\,\xi\in\mathbb{R}^{2}.$
(2.7)

It is well known in [8] that those are all the solutions of the equation

$\left\\{\begin{aligned} &-\Delta u=e^{u}\,\
\,\,\,\,\textrm{in}\,\,\,\,\mathbb{R}^{2},\\\
&\int_{\mathbb{R}^{2}}e^{u}<+\infty.\end{aligned}\right.$ (2.8)

The configuration space for $m$ concentration points
$\xi=(\xi_{1},\ldots,\xi_{m})$ we try to seek is the following

$\displaystyle\mathcal{O}_{p}:=\left\\{\,\xi=(\xi_{1},\ldots,\xi_{m})\in\Omega^{l}\times(\partial\Omega)^{m-l}\left|\,\min_{i,j=1,\ldots,m,\,i\neq
j}|\xi_{i}-\xi_{j}|>\frac{1}{p^{\kappa}},\,\,\,\,\,\,\,\min_{1\leq i\leq
l}\text{dist}(\xi_{i},\partial\Omega)>\frac{1}{p^{\kappa}}\right.\right\\},$
(2.9)

where $l\in\\{0,\ldots,m\\}$ and $\kappa$ is given by

$\displaystyle\kappa=2(m^{2}+1).$ (2.10)

Let $m\in\mathbb{N}^{*}$ and $\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$
be fixed. For numbers $\mu_{i}>0$, $i=1,\ldots,m$, yet to be chosen, we define

$\displaystyle\gamma=p^{\frac{p}{p-1}}e^{-\frac{p}{2(p-1)}},\,\,\qquad\qquad\,\,\varepsilon=e^{-\frac{1}{4}p},\,\,\qquad\qquad\
\,\,\delta_{i}=\varepsilon\mu_{i},$ (2.11)

and

$\displaystyle
U_{i}(x)=\frac{1}{\gamma\mu_{i}^{2/(p-1)}}\left[U_{\delta_{i},\xi_{i}}(x)+\frac{1}{p}\omega_{1}\left(\frac{x-\xi_{i}}{\delta_{i}}\right)+\frac{1}{p^{2}}\omega_{2}\left(\frac{x-\xi_{i}}{\delta_{i}}\right)\right].$
(2.12)

Here, $\omega_{j}$, $j=1,2$, are radial solutions of

$\displaystyle\Delta\omega_{j}+\frac{8}{(1+|z|^{2})^{2}}\omega_{j}=\frac{8}{(1+|z|^{2})^{2}}f_{j}(|z|)\quad\,\textrm{in}\,\,\,\,\mathbb{R}^{2},$
(2.13)

with

$\displaystyle
f_{1}=\frac{1}{2}U_{1,0}^{2},\,\quad\quad\,f_{2}=\omega_{1}U_{1,0}-\frac{1}{3}U_{1,0}^{3}-\frac{1}{2}\omega_{1}^{2}-\frac{1}{8}U_{1,0}^{4}+\frac{1}{2}\omega_{1}U_{1,0}^{2},$
(2.14)

having asymptotic expansions

$\displaystyle\omega_{j}(r)=\frac{C_{j}}{2}\log(1+r^{2})+O\left(\frac{1}{1+r}\right),\,\,\quad\,\,\partial_{r}\omega_{j}(r)=\frac{C_{j}r}{1+r^{2}}+O\left(\frac{1}{1+r^{2}}\right)\quad\,\,\,\textrm{as}\,\,\,r\rightarrow+\infty,\,\,\,r=|z|,$
(2.15)

where

$\displaystyle
C_{j}=8\int_{0}^{\infty}t\frac{t^{2}-1}{(t^{2}+1)^{3}}f_{j}(t)dt,$ (2.16)

in particular,

$\displaystyle\omega_{1}(z)=\frac{1}{2}U_{1,0}^{2}(z)+6\log(|z|^{2}+1)+\frac{2\log
8-10}{|z|^{2}+1}+\frac{|z|^{2}-1}{|z|^{2}+1}\times\left\\{-\frac{1}{2}\log^{2}8+2\log^{2}(|z|^{2}+1)\right.$
$\displaystyle\left.+4\int_{|z|^{2}}^{\infty}\frac{ds}{s+1}\log\frac{s+1}{s}-8\log|z|\log(|z|^{2}+1)\right\\},\,\,\qquad\qquad\qquad\qquad\qquad\qquad\qquad\,\,$
(2.17)

and

$\displaystyle C_{1}=12-4\log 8$ (2.18)

(see [7, 21]). We now approximate the solution of problem (1.6) by

$\displaystyle
U_{\xi}(x):=\sum_{i=1}^{m}PU_{i}(x)=\sum_{i=1}^{m}\big{[}U_{i}(x)+H_{i}(x)\big{]},$
(2.19)

where $H_{i}$ is a correction term defined as the solution of

$\displaystyle\left\\{\begin{aligned} &-\Delta_{a}H_{i}+H_{i}=\nabla\log
a(x)\nabla
U_{i}-U_{i}\,\,\,\,\,\,\textrm{in}\,\,\,\,\,\,\Omega,\\\\[2.84526pt]
&\frac{\partial H_{i}}{\partial\nu}=-\frac{\partial
U_{i}}{\partial\nu}\,\,\qquad\qquad\qquad\qquad\qquad\,\,\,\textrm{on}\,\,\,\,\partial\Omega.\end{aligned}\right.$
(2.20)

In order to understand the asymptotic behavior of the correction term $H_{i}$,
let us make first the convention that

$\displaystyle c_{i}=\left\\{\begin{aligned}
&8\pi,\qquad\textrm{if}\quad\,\xi_{i}\in\Omega,\\\\[2.84526pt]
&4\pi,\qquad\textrm{if}\quad\xi_{i}\in\partial\Omega.\end{aligned}\right.$
(2.21)

Furthermore, we have the following result, whose proof is postponed to the
Appendix.

Lemma 2.4. For any $0<\alpha<1$ and
$\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$, then we have

$\displaystyle
H_{i}(x)=\frac{1}{\gamma\mu_{i}^{2/(p-1)}}\left[\left(1-\frac{C_{1}}{4p}-\frac{C_{2}}{4p^{2}}\right)c_{i}H(x,\xi_{i})-\log(8\delta_{i}^{2})+\left(\frac{C_{1}}{p}+\frac{C_{2}}{p^{2}}\right)\log\delta_{i}+O\left(\delta_{i}^{\alpha/2}\right)\right]$
(2.22)

uniformly in $\overline{\Omega}$, where $H$ is the regular part of Green’s
function defined in (1.9).

From Lemma 2.4 we have that away from each point $\xi_{i}$, namely
$|x-\xi_{i}|\geq 1/p^{2\kappa}$ for any $i=1,\ldots,m$,

$\displaystyle
U_{\xi}(x)=\sum\limits_{i=1}^{m}\frac{1}{\gamma\mu_{i}^{2/(p-1)}}\left[\left(1-\frac{C_{1}}{4p}-\frac{C_{2}}{4p^{2}}\right)c_{i}G(x,\xi_{i})+\frac{1}{p}O\left(\delta_{i}^{\alpha/2}\right)\right].$
(2.23)

While for $|x-\xi_{i}|<1/p^{2\kappa}$ with some $i$, by (2.7), (2.12), (2.15),
(2.22) and the fact that $H(\cdot,\xi_{k})\in C^{\alpha}(\overline{\Omega})$
for any $\xi_{k}\in\overline{\Omega}$ and any $\alpha\in(0,1)$ we obtain

$\displaystyle PU_{i}(x)=$
$\displaystyle\frac{1}{\gamma\mu_{i}^{2/(p-1)}}\left[U_{1,0}\left(\frac{x-\xi_{i}}{\delta_{i}}\right)+\frac{1}{p}\omega_{1}\left(\frac{x-\xi_{i}}{\delta_{i}}\right)+\frac{1}{p^{2}}\omega_{2}\left(\frac{x-\xi_{i}}{\delta_{i}}\right)+\left(1-\frac{C_{1}}{4p}-\frac{C_{2}}{4p^{2}}\right)c_{i}H(\xi_{i},\xi_{i})\right.$
$\displaystyle\left.-\log(8\delta_{i}^{4})+\left(\frac{C_{1}}{p}+\frac{C_{2}}{p^{2}}\right)\log\delta_{i}+O\left(|x-\xi_{i}|^{\alpha}+\delta_{i}^{\alpha/2}\right)\right],$

and for any $k\neq i$,

$\displaystyle PU_{k}(x)=$
$\displaystyle\frac{1}{\gamma\mu_{k}^{2/(p-1)}}\left[\log\frac{8\delta_{k}^{2}}{(\delta_{k}^{2}+|x-\xi_{k}|^{2})^{2}}+\frac{1}{p}\omega_{1}\left(\frac{x-\xi_{k}}{\delta_{k}}\right)+\frac{1}{p^{2}}\omega_{2}\left(\frac{x-\xi_{k}}{\delta_{k}}\right)\right.$
$\displaystyle\left.+\left(1-\frac{C_{1}}{4p}-\frac{C_{2}}{4p^{2}}\right)c_{k}H(x,\xi_{k})-\log(8\delta_{k}^{2})+\left(\frac{C_{1}}{p}+\frac{C_{2}}{p^{2}}\right)\log\delta_{k}+O\left(\delta_{k}^{\alpha/2}\right)\right]$
$\displaystyle=$
$\displaystyle\frac{1}{\gamma\mu_{k}^{2/(p-1)}}\left[\left(1-\frac{C_{1}}{4p}-\frac{C_{2}}{4p^{2}}\right)c_{k}G(\xi_{i},\xi_{k})+O\left(|x-\xi_{i}|^{\alpha}+\delta_{k}^{\alpha/2}\right)\right].$

Hence for $|x-\xi_{i}|<1/p^{2\kappa}$,

$\displaystyle
U_{\xi}(x)=\frac{1}{\gamma\mu_{i}^{2/(p-1)}}\left[\,p+U_{1,0}\left(\frac{x-\xi_{i}}{\delta_{i}}\right)+\frac{1}{p}\omega_{1}\left(\frac{x-\xi_{i}}{\delta_{i}}\right)+\frac{1}{p^{2}}\omega_{2}\left(\frac{x-\xi_{i}}{\delta_{i}}\right)+O\left(|x-\xi_{i}|^{\alpha}+\sum_{k=1}^{m}\delta_{k}^{\alpha/2}\right)\right]$
(2.24)

is an appropriate approximation for a solution of problem (1.6) provided that
the concentration parameters $\mu_{i}$, $i=1,\ldots,m$, are the solution of
the nonlinear system

$\displaystyle\log\big{(}8\mu_{i}^{4}\big{)}=\left(1-\frac{C_{1}}{4p}-\frac{C_{2}}{4p^{2}}\right)c_{i}H(\xi_{i},\xi_{i})+\left(\frac{C_{1}}{p}+\frac{C_{2}}{p^{2}}\right)\log\delta_{i}+\left(1-\frac{C_{1}}{4p}-\frac{C_{2}}{4p^{2}}\right)\sum_{k\neq
i}^{m}\left(\frac{\mu_{i}}{\mu_{k}}\right)^{2/(p-1)}c_{k}G(\xi_{i},\xi_{k}).$
(2.25)

Indeed, the parameters $\mu=(\mu_{1},\ldots,\mu_{m})$ are well defined in
system (2.25), which is stated as follows and proved in the Appendix.

Lemma 2.5. For any points $\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$ and
any $p>1$ large enough, system (2.25) has a unique solution
$\mu=(\mu_{1},\ldots,\mu_{m})$ satisfying

$\displaystyle 1/C\leq\mu_{i}\leq Cp^{\kappa}\qquad\forall\,\,\,i=1,\ldots,m,$
(2.26)

for some $C>0$. Moreover, for any $i=1,\ldots,m$, one has

$\displaystyle\big{|}D_{\xi}\log\mu_{i}\big{|}\leq Cp^{\kappa},$ (2.27)

and

$\displaystyle\mu_{i}=e\large^{-\frac{3}{4}+\frac{1}{4}c_{i}H(\xi_{i},\xi_{i})+\frac{1}{4}\sum\large_{k=1,\,k\neq
i}^{m}c_{k}G(\xi_{i},\xi_{k})}\left[\,1+O\left(\frac{\log^{2}p}{p}\right)\right].$
(2.28)

Remark 2.6. Observe that for $|x-\xi_{i}|=\delta_{i}|z|<1/p^{2\kappa}$, by
(2.7), (2.11), (2.15) and (2.26),

$\displaystyle
p+U_{1,0}(z)+\frac{1}{p}\omega_{1}(z)+\frac{1}{p^{2}}\omega_{2}(z)\geq
p-2\log\left(1+|z|^{2}\right)+O\left(1\right)\geq 8\kappa\log
p+4\log\mu_{i}+O\left(1\right)>7\kappa\log p.$

Hence by (2.24), we can easily get that $0<U_{\xi}\leq 2\sqrt{e}$ in
$B_{1/p^{2\kappa}}(\xi_{i})$, and
$\sup_{B_{1/p^{2\kappa}}(\xi_{i})}U_{\xi}\rightarrow\sqrt{e}$ as
$p\rightarrow+\infty$. Moreover, by the maximum principle, we see that
$G(x,\xi_{i})>0$ over $\overline{\Omega}$ and thus by (2.23), $U_{\xi}$ is a
positive, uniformly bounded function over $\overline{\Omega}$. In conclusion,
$0<U_{\xi}\leq 2\sqrt{e}$ over $\overline{\Omega}$.

Let us perform the change of variables

$\displaystyle\upsilon(y)=\varepsilon^{2/(p-1)}u(\varepsilon
y),\,\,\,\,\quad\forall\,\,\,y\in\Omega_{p}:=e^{p/4}\Omega.$ (2.29)

Then by the definition of $\varepsilon$ in (2.11), $u(x)$ solves equation
(1.6) if and only if the function $\upsilon(y)$ satisfies

$\begin{array}[]{ll}\left\\{\begin{aligned} &-\Delta_{a(\varepsilon
y)}\upsilon+\varepsilon^{2}\upsilon=\upsilon^{p},\,\,\,\,\,\,\upsilon>0\,\,\,\,\,\,\textrm{in}\,\,\,\,\,\Omega_{p},\\\\[2.84526pt]
&\frac{\partial\upsilon}{\partial\nu}=0\qquad\qquad\qquad\qquad\qquad\quad\,\textrm{on}\,\,\,\partial\Omega_{p}.\end{aligned}\right.\end{array}$
(2.30)

We write $\xi_{i}^{\prime}=\xi_{i}/\varepsilon$, $i=1,\ldots,m$ and define the
initial approximate solution of (2.30) as

$\displaystyle V_{\xi^{\prime}}(y)=\varepsilon^{2/(p-1)}U_{\xi}(\varepsilon
y),$ (2.31)

with $\xi^{\prime}=(\xi_{1}^{\prime},\ldots,\xi_{m}^{\prime})$ and $U_{\xi}$
defined in (2.19). Let us set

$\displaystyle S_{p}(\upsilon)=-\Delta_{a(\varepsilon
y)}\upsilon+\varepsilon^{2}\upsilon-\upsilon^{p}_{+},\,\qquad\,\textrm{where}\,\,\,\,\upsilon_{+}=\max\\{\upsilon,0\\},$

and introduce the functional

$\displaystyle I_{p}(\upsilon)=\frac{1}{2}\int_{\Omega_{p}}a(\varepsilon
y)\left(|\nabla\upsilon|^{2}+\varepsilon^{2}\upsilon^{2}\right)dy-\frac{1}{p+1}\int_{\Omega_{p}}a(\varepsilon
y)\upsilon^{p+1}_{+}dy,\,\quad\,\upsilon\in H^{1}(\Omega_{p}),$ (2.32)

whose nontrivial critical points are solutions of problem (2.30). In fact, by
the maximum principle, problem (2.30) is equivalent to

$\displaystyle S_{p}(\upsilon)=0,\,\quad\,\upsilon_{+}\not\equiv
0\,\quad\,\textrm{in}\,\,\,\,\Omega_{p},\,\qquad\qquad\,\frac{\partial\upsilon}{\partial\nu}=0\,\quad\,\textrm{on}\,\,\,\partial\Omega_{p}.$

We will look for solutions of problem (2.30) in the form
$\upsilon=V_{\xi^{\prime}}+\phi$, where $\phi$ will represent a higher-order
correction in the expansion of $\upsilon$. Observe that

$\displaystyle
S_{p}(V_{\xi^{\prime}}+\phi)=\mathcal{L}(\phi)+R_{\xi^{\prime}}+N(\phi)=0,$
(2.33)

where

$\displaystyle\mathcal{L}(\phi)=-\Delta_{a(\varepsilon
y)}\phi+\varepsilon^{2}\phi-
W_{\xi^{\prime}}\phi\,\quad\,\textrm{with}\,\quad\,W_{\xi^{\prime}}=pV_{\xi^{\prime}}^{p-1},$
(2.34)

and

$\displaystyle R_{\xi^{\prime}}=-\Delta_{a(\varepsilon
y)}V_{\xi^{\prime}}+\varepsilon^{2}V_{\xi^{\prime}}-V_{\xi^{\prime}}^{p},\quad\quad\quad
N(\phi)=-\big{[}(V_{\xi^{\prime}}+\phi)_{+}^{p}-V_{\xi^{\prime}}^{p}-pV_{\xi^{\prime}}^{p-1}\phi\big{]}.$
(2.35)

In terms of $\phi$, problem (2.30) becomes

$\displaystyle\left\\{\begin{aligned}
&\mathcal{L}(\phi)=-\big{[}R_{\xi^{\prime}}+N(\phi)\big{]}\quad\textrm{in}\,\,\,\,\,\,\Omega_{p},\\\
&\frac{\partial\phi}{\partial\nu}=0\quad\qquad\qquad\qquad\,\,\,\,\textrm{on}\,\,\,\,\partial\Omega_{p}.\end{aligned}\right.$
(2.36)

For any $\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$ and $h\in
L^{\infty}(\Omega_{p})$, let us introduce a weighted $L^{\infty}$-norm defined
as

$\displaystyle\|h\|_{*}=\sup_{y\in\Omega_{p}}\left|\left(\sum\limits_{i=1}^{m}\frac{\mu_{i}^{\sigma}}{(\mu_{i}+|y-\xi^{\prime}_{i}|)^{2+\sigma}}+\varepsilon^{2}\right)^{-1}h(y)\right|,$
(2.37)

where $\sigma>0$ is small but fixed independent of $p$. With respect to the
$\|\cdot\|_{*}$-norm, the error term $R_{\xi^{\prime}}$ defined in (2.35) can
be estimated as follows.

Proposition 2.7. Let $m$ be a positive integer. There exist constants $C>0$
and $p_{m}>1$ such that for any
$\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$ and any $p>p_{m}$,

$\displaystyle\|R_{\xi^{\prime}}\|_{*}\leq\frac{C}{p^{4}}.$ (2.38)

###### Proof.

Observe that, by (2.19), (2.20) and (2.31),

$\displaystyle-\Delta_{a(\varepsilon
y)}V_{\xi^{\prime}}+\varepsilon^{2}V_{\xi^{\prime}}=\varepsilon^{2}\sum_{i=1}^{m}\varepsilon^{2/(p-1)}\left[-\Delta_{a}\big{(}U_{i}+H_{i}\big{)}+\big{(}U_{i}+H_{i}\big{)}\right]=-\varepsilon^{2p/(p-1)}\sum_{i=1}^{m}\Delta
U_{i}.$

Then by (2.7), (2.11), (2.12) and (2.13),

$\displaystyle-\Delta_{a(\varepsilon
y)}V_{\xi^{\prime}}+\varepsilon^{2}V_{\xi^{\prime}}=\sum_{i=1}^{m}\frac{1}{p^{p/(p-1)}\mu_{i}^{2p/(p-1)}}e^{U_{1,0}(z)}\left[1-\frac{1}{p}f_{1}(z)-\frac{1}{p^{2}}f_{2}(z)+\frac{1}{p}\omega_{1}(z)+\frac{1}{p^{2}}\omega_{2}(z)\right]$
(2.39)

with $z=(\varepsilon y-\xi_{i})/\delta_{i}$. By (2.7), (2.15) and (2.26) we
get, if $|\varepsilon y-\xi_{i}|=\delta_{i}|z|\geq 1/p^{2\kappa}$ for any
$i=1,\ldots,m$, then

$\displaystyle U_{1,0}(z)=-p+O\left(\log
p\right),\,\qquad\quad\,\omega_{j}(z)=\frac{1}{4}C_{j}p+O\left(\log
p\right),\,\,\quad\,j=1,\,2,$

and hence, by (2.14) and (2.39),

$\displaystyle-\Delta_{a(\varepsilon
y)}V_{\xi^{\prime}}+\varepsilon^{2}V_{\xi^{\prime}}=\sum_{i=1}^{m}\frac{e^{U_{1,0}(z)}}{p^{p/(p-1)}\mu_{i}^{2p/(p-1)}}O\left(p^{2}\right).$
(2.40)

On the other hand, in the same region, by (2.23) and (2.31) we get

$\displaystyle V_{\xi^{\prime}}^{p}=O\left(\frac{1}{p}\left(\frac{\log
p}{p}\right)^{p}\right),$ (2.41)

which, together with (2.11), (2.26) and (2.40), implies

$\displaystyle\left|\left(\sum\limits_{i=1}^{m}\frac{\mu_{i}^{\sigma}}{(\mu_{i}+|y-\xi^{\prime}_{i}|)^{2+\sigma}}+\varepsilon^{2}\right)^{-1}R_{\xi^{\prime}}(y)\right|\leq
C\left[p\sum_{i=1}^{m}\left|\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right|^{\sigma-2}+\frac{1}{p}\left(\frac{\sqrt{e}\log
p}{p}\right)^{p}\right]=o\left(e^{-p/4}\right).$ (2.42)

Let us fix the index $i\in\\{1,\ldots,m\\}$ and the region $|\varepsilon
y-\xi_{i}|=\delta_{i}|z|<1/p^{2\kappa}$. By (2.24), (2.31) and the relation

$\displaystyle\left(\frac{p\varepsilon^{2/(p-1)}}{\gamma\mu_{i}^{2/(p-1)}}\right)^{p}=\frac{1}{p^{p/(p-1)}\mu_{i}^{2p/(p-1)}},$
(2.43)

we get, for $|\varepsilon y-\xi_{i}|=\delta_{i}|z|<1/p^{2\kappa}$,

$\displaystyle
V^{p}_{\xi^{\prime}}(y)=\frac{1}{p^{p/(p-1)}\mu_{i}^{2p/(p-1)}}\left\\{1+\frac{U_{1,0}(z)}{p}+\frac{\omega_{1}(z)}{p^{2}}+\frac{1}{p^{3}}\left[\omega_{2}(z)+O\left(p^{2}\delta_{i}^{\alpha}|z|^{\alpha}+p^{2}\sum_{k=1}^{m}\delta_{k}^{\alpha/2}\right)\right]\right\\}^{p}.$
(2.44)

From Taylor expansions of the exponential and logarithmic function

$\displaystyle\left(1+\frac{a}{p}+\frac{b}{p^{2}}+\frac{c}{p^{3}}\right)^{p}=e^{a}\left[1+\frac{1}{p}\left(b-\frac{a^{2}}{2}\right)+\frac{1}{p^{2}}\left(c-ab+\frac{a^{3}}{3}+\frac{b^{2}}{2}-\frac{a^{2}b}{2}+\frac{a^{4}}{8}\right)\right.$
$\displaystyle\left.+O\left(\frac{\log^{6}(|z|+2)}{p^{3}}\right)\right],\,\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\,\,\,$
(2.45)

which holds for $|z|\leq Ce^{p/8}$ provided $-4\log(|z|+2)\leq a(z)\leq C$ and
$|b(z)|+|c(z)|\leq C\log(|z|+2)$, we have that for $|\varepsilon
y-\xi_{i}|=\delta_{i}|z|\leq\sqrt{\delta_{i}}/p^{2\kappa}$,

$\displaystyle
V^{p}_{\xi^{\prime}}(y)=\frac{1}{p^{p/(p-1)}\mu_{i}^{2p/(p-1)}}e^{U_{1,0}(z)}\left[\,1+\frac{1}{p}\left(\omega_{1}-\frac{1}{2}U_{1,0}^{2}\right)(z)+\frac{1}{p^{2}}\left(\omega_{2}-\omega_{1}U_{1,0}+\frac{1}{3}U_{1,0}^{3}+\frac{1}{2}\omega_{1}^{2}\right.\right.$
$\displaystyle\left.\left.-\frac{1}{2}\omega_{1}U_{1,0}^{2}+\frac{1}{8}U_{1,0}^{4}\right)(z)+O\left(\frac{\log^{6}(|z|+2)}{p^{3}}+\delta_{i}^{\alpha}|z|^{\alpha}+\sum_{k=1}^{m}\delta_{k}^{\alpha/2}\right)\right],\,\qquad\qquad\qquad\qquad$
(2.46)

which combined with (2.14) and (2.39) gives

$\displaystyle R_{\xi^{\prime}}(y)=-\Delta_{a(\varepsilon
y)}V_{\xi^{\prime}}+\varepsilon^{2}V_{\xi^{\prime}}-V_{\xi^{\prime}}^{p}=\frac{1}{p^{p/(p-1)}\mu_{i}^{2p/(p-1)}}e^{U_{1,0}(z)}O\left(\frac{\log^{6}(|z|+2)}{p^{3}}+\delta_{i}^{\alpha}|z|^{\alpha}+\sum_{k=1}^{m}\delta_{k}^{\alpha/4}\right).$

Hence, in this region we get

$\displaystyle\left|\left(\sum\limits_{i=1}^{m}\frac{\mu_{i}^{\sigma}}{(\mu_{i}+|y-\xi^{\prime}_{i}|)^{2+\sigma}}+\varepsilon^{2}\right)^{-1}R_{\xi^{\prime}}(y)\right|\leq\frac{C}{p^{4}\big{(}\big{|}\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\big{|}+1\big{)}^{2-\sigma}}\log^{6}\left(\left|\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right|+2\right)=O\left(\frac{1}{p^{4}}\right).$
(2.47)

Finally, in the remaining region $\sqrt{\delta_{i}}/p^{2\kappa}<|\varepsilon
y-\xi_{i}|=\delta_{i}|z|<1/p^{2\kappa}$, we have that, by (2.14) and (2.39),

$\displaystyle-\Delta_{a(\varepsilon
y)}V_{\xi^{\prime}}+\varepsilon^{2}V_{\xi^{\prime}}=\frac{e^{U_{1,0}(z)}}{p^{p/(p-1)}\mu_{i}^{2p/(p-1)}}O\left(p^{2}\right),$
(2.48)

and by (2.44),

$\displaystyle
V^{p}_{\xi^{\prime}}(y)=\frac{e^{U_{1,0}(z)}}{p^{p/(p-1)}\mu_{i}^{2p/(p-1)}}O\left(1\right),$
(2.49)

since $(1+\frac{s}{p})^{p}\leq e^{s}$. Thus, in this region,

$\displaystyle\left|\left(\sum\limits_{i=1}^{m}\frac{\mu_{i}^{\sigma}}{(\mu_{i}+|y-\xi^{\prime}_{i}|)^{2+\sigma}}+\varepsilon^{2}\right)^{-1}R_{\xi^{\prime}}(y)\right|\leq
Cp\left|\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right|^{\sigma-2}=o\left(e^{-p/8}\right).$
(2.50)

Combining (2.37), (2.42), (2.47) with (2.50), we conclude that estimate (2.38)
holds. ∎

##  3\. Analysis of the linearized operator

In this section, we prove bounded invertibility of the operator $\mathcal{L}$,
uniformly on $\xi\in\mathcal{O}_{p}$, by using the weighted $L^{\infty}$-norm
defined in (2.37). Let us recall that
$\mathcal{L}(\phi)=-\Delta_{a(\varepsilon y)}\phi+\varepsilon^{2}\phi-
W_{\xi^{\prime}}\phi$, where $W_{\xi^{\prime}}=pV_{\xi^{\prime}}^{p-1}$. As in
Proposition 2.7, we have the following expansions with respect to the
potential $W_{\xi^{\prime}}$.

Lemma 3.1. Let $m$ be a positive integer. There exist constants $D_{0}>0$ and
$p_{m}>1$ such that

$\displaystyle W_{\xi^{\prime}}(y)\leq
D_{0}\sum_{i=1}^{m}\frac{1}{\mu_{i}^{2}}e^{U_{1,0}\big{(}\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\big{)}}$
(3.1)

for any points $\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$ and any $p\geq
p_{m}$. Furthermore,

$\displaystyle
W_{\xi^{\prime}}(y)=\frac{8}{\mu_{i}^{2}(1+|z|^{2})^{2}}\left[1+\frac{1}{p}\left(\omega_{1}-U_{1,0}-\frac{1}{2}U_{1,0}^{2}\right)(z)+O\left(\frac{\log^{4}\big{(}|z|+2\big{)}}{p^{2}}\right)\right]$
(3.2)

for any $|\varepsilon y-\xi_{i}|\leq\sqrt{\delta_{i}}/p^{2\kappa}$, where
$z=(y-\xi^{\prime}_{i})/\mu_{i}$.

Proof. If $|\varepsilon y-\xi_{i}|=\delta_{i}|z|<1/p^{2\kappa}$ for some
$i=1,\ldots,m$, by (2.24), (2.31) and (2.43),

$\displaystyle W_{\xi^{\prime}}(y)=$ $\displaystyle
p\left(\frac{\varepsilon^{2/(p-1)}}{\gamma\mu_{i}^{2/(p-1)}}\right)^{p-1}\left[\,p+U_{1,0}(z)+\frac{1}{p}\omega_{1}(z)+\frac{1}{p^{2}}\omega_{2}(z)+O\left(\delta_{i}^{\alpha}|z|^{\alpha}+\sum_{k=1}^{m}\delta_{k}^{\alpha/2}\right)\right]^{p-1}$
$\displaystyle=$
$\displaystyle\frac{1}{\mu_{i}^{2}}\left[1+\frac{1}{p}U_{1,0}(z)+\frac{1}{p^{2}}\omega_{1}(z)+\frac{1}{p^{3}}\omega_{2}(z)+\frac{1}{p}O\left(\delta_{i}^{\alpha}|z|^{\alpha}+\sum_{k=1}^{m}\delta_{k}^{\alpha/2}\right)\right]^{p-1}.$

In this region, using the fact that $(1+a/p)^{p-1}\leq e^{(p-1)a/p}$ and
$U_{1,0}(z)\geq-p+O\big{(}\log p\big{)}$, we get

$\displaystyle
W_{\xi^{\prime}}(y)\leq\frac{C}{\mu_{i}^{2}}e^{U_{1,0}(z)}e^{-U_{1,0}(z)/p}=O\left(\frac{1}{\mu_{i}^{2}}e^{U_{1,0}(z)}\right).$

In particular, from a slight modification of formula (2), namely

$\displaystyle\left(1+\frac{a}{p}+\frac{b}{p^{2}}+\frac{c}{p^{3}}\right)^{p-1}=e^{a}\left[1+\frac{1}{p}\left(b-a-\frac{a^{2}}{2}\right)+O\left(\frac{\log^{4}(|z|+2)}{p^{2}}\right)\right],$

we conclude that if $|\varepsilon
y-\xi_{i}|=\delta_{i}|z|\leq\sqrt{\delta_{i}}/p^{2\kappa}$, then

$\displaystyle
W_{\xi^{\prime}}(y)=\frac{1}{\mu_{i}^{2}}e^{U_{1,0}(z)}\left[1+\frac{1}{p}\left(\omega_{1}-U_{1,0}-\frac{1}{2}U_{1,0}^{2}\right)(z)+O\left(\frac{\log^{4}(|z|+2)}{p^{2}}\right)\right].$

Additionally, if $|\varepsilon y-\xi_{i}|=\delta_{i}|z|\geq 1/p^{2\kappa}$ for
all $i$, then by (2.11), (2.23), (2.26) and (2.31),

$\displaystyle\,\,\quad\qquad\qquad\,\,W_{\xi^{\prime}}(y)=p\varepsilon^{2}U_{\xi}^{p-1}(\varepsilon
y)\leq p\varepsilon^{2}\left(\frac{C\log
p}{\gamma}\right)^{p-1}=O\left(\left(\frac{\log
p}{p}\right)^{p-1}\right)=o\left(\frac{1}{\mu_{i}^{2}}e^{U_{1,0}(z)}\right).\qquad\qquad\quad\,\,\,\,\square$

Remark 3.2. As for $W_{\xi^{\prime}}$, we mention that if $|\varepsilon
y-\xi_{i}|<1/p^{2\kappa}$ for some $i=1,\ldots,m$, then

$\displaystyle
p\left[V_{\xi^{\prime}}(y)+O\left(\frac{1}{p^{3}}\right)\right]^{p-2}\leq
Cp\left(\frac{p\varepsilon^{2/(p-1)}}{\gamma\mu_{i}^{2/(p-1)}}\right)^{p-2}e^{\frac{p-2}{p}U_{1,0}\big{(}\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\big{)}}=O\left(\frac{1}{\mu_{i}^{2}}e^{U_{1,0}\big{(}\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\big{)}}\right).$

Since this estimate is true if $|\varepsilon y-\xi_{i}|\geq 1/p^{2\kappa}$ for
all $i$, we get

$\displaystyle
p\left[V_{\xi^{\prime}}(y)+O\left(\frac{1}{p^{3}}\right)\right]^{p-2}\leq
C\sum\limits_{i=1}^{m}\frac{1}{\mu_{i}^{2}}e^{U_{1,0}\big{(}\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\big{)}}.$
(3.3)

Let

$\displaystyle
Z_{0}(z)=\frac{|z|^{2}-1}{|z|^{2}+1},\,\,\quad\qquad\quad\,\,Z_{j}(z)=\frac{z_{j}}{|z|^{2}+1},\,\,\,\,j=1,\,2.$
(3.4)

It is well known (see [6, 8]) that

  * •

any bounded solution to

$\displaystyle\Delta\phi+\frac{8}{(1+|z|^{2})^{2}}\phi=0\,\quad\textrm{in}\,\,\,\mathbb{R}^{2},$
(3.5)

is a linear combination of $Z_{j}$, $j=0,1,2$;

  * •

any bounded solution to

$\displaystyle\Delta\phi+\frac{8}{(1+|z|^{2})^{2}}\phi=0\,\quad\textrm{in}\,\,\,\mathbb{R}^{2}_{+},\,\qquad\qquad\,\frac{\partial\phi}{\partial\nu}=0\,\quad\textrm{on}\,\,\,\partial\mathbb{R}^{2}_{+},$
(3.6)

where $\mathbb{R}^{2}_{+}=\\{(z_{1},z_{2}):\,z_{2}>0\\}$, is a linear
combination of $Z_{j}$, $j=0,1$.

Now we consider the following linear problem: given $h\in
C(\overline{\Omega}_{p})$ and points
$\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$, we find a function $\phi\in
H^{2}(\Omega_{p})$ and scalars $c_{ij}\in\mathbb{R}$, $i=1,\ldots,m$,
$j=1,J_{i}$, such that

$\left\\{\begin{aligned} &\mathcal{L}(\phi)=-\Delta_{a(\varepsilon
y)}\phi+\varepsilon^{2}\phi-W_{\xi^{\prime}}\phi=h+\frac{1}{a(\varepsilon
y)}\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{J_{i}}c_{ij}\chi_{i}\,Z_{ij}\,\,\
\,\,\textrm{in}\,\,\,\,\,\,\Omega_{p},\\\
&\frac{\partial\phi}{\partial\nu}=0\,\,\,\,\,\,\,\,\ \ \ \ \ \ \ \ \ \ \ \ \
\,\,\qquad\qquad\qquad\quad\qquad\qquad\qquad\qquad\ \,\ \ \ \,\,\,\ \,\
\,\textrm{on}\,\,\,\,\partial\Omega_{p},\\\\[2.84526pt]
&\int_{\Omega_{p}}\chi_{i}\,Z_{ij}\phi=0\,\qquad\qquad\qquad\quad\qquad\qquad\qquad\forall\,\,i=1,\ldots,m,\,\,\,j=1,J_{i},\end{aligned}\right.$
(3.7)

where $J_{i}=2$ if $i=1,\ldots,l$ while $J_{i}=1$ if $i=l+1,\ldots,m$, and
$Z_{ij}$, $\chi_{i}$, are defined as follows.

Let $\chi:\mathbb{R}\rightarrow[0,1]$ be 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$.

For $i=1,\ldots,l$ (corresponding to interior spike case), we define

$\displaystyle\chi_{i}(y)=\chi\left(\frac{|y-\xi^{\prime}_{i}|}{\mu_{i}}\right),\,\qquad\,Z_{ij}(y)=\frac{1}{\mu_{i}}Z_{j}\left(\frac{y-\xi_{i}^{\prime}}{\mu_{i}}\right),\,\quad\,j=0,1,2.$
(3.8)

For $i=l+1,\ldots,m$ (corresponding to boundary spike case), we have to
strengthen the boundary first. More precisely, at the boundary point
$\xi_{i}\in\partial\Omega$, we define a rotation map
$A_{i}:\mathbb{R}^{2}\mapsto\mathbb{R}^{2}$ such that
$A_{i}\nu_{\Omega}(\xi_{i})=\nu_{\mathbb{R}_{+}^{2}}(0)$. Let
$\mathcal{G}(x_{1})$ be the defining function for the boundary
$A_{i}(\partial\Omega-\\{\xi_{i}\\})$ in a small neighborhood
$B_{\delta}(0,0)$ of the origin, that is, there exist $R_{1}>0$, $\delta>0$
small and a smooth function $\mathcal{G}:(-R_{1},R_{1})\mapsto\mathbb{R}$
satisfying $\mathcal{G}(0)=0$, $\mathcal{G}^{\prime}(0)=0$ and such that
$A_{i}(\Omega-\\{\xi_{i}\\})\cap
B_{\delta}(0,0)=\\{(x_{1},x_{2}):\,-R_{1}<x_{1}<R_{1},\,x_{2}>\mathcal{G}(x_{1})\\}\cap
B_{\delta}(0,0)$. Furthermore, we consider the flattening change of variables
$F_{i}:B_{\delta}(0,0)\cap\overline{A_{i}(\Omega-\\{\xi_{i}\\})}\mapsto\mathbb{R}^{2}$
be defined by $F_{i}=(F_{i1},F_{i2})$, where

$\displaystyle
F_{i1}=x_{1}+\frac{x_{2}-\mathcal{G}(x_{1})}{\,1+|\mathcal{G}^{\prime}(x_{1})|^{2}\,}\mathcal{G}^{\prime}(x_{1})\
\qquad\textrm{and}\ \qquad\,F_{i2}=x_{2}-\mathcal{G}(x_{1}).$ (3.9)

Then for any $i=l+1,\ldots,m$, we set

$\displaystyle
F_{i}^{p}(y)=e^{p/4}F_{i}\big{(}A_{i}(e^{-p/4}y-\xi_{i})\big{)}=\frac{1}{\varepsilon}F_{i}\big{(}A_{i}(\varepsilon
y-\xi_{i})\big{)},$ (3.10)

and define

$\displaystyle\chi_{i}(y)=\chi\left(\frac{1}{\mu_{i}}|F_{i}^{p}(y)|\right),\,\qquad\,Z_{ij}(y)=\frac{1}{\mu_{i}}Z_{j}\left(\frac{1}{\mu_{i}}F_{i}^{p}(y)\right),\,\quad\,j=0,1.$
(3.11)

It is important to note that $F^{p}_{i}$, $i=l+1,\ldots,m$, preserves the
Neumann boundary condition and

$\displaystyle\Delta_{a(\varepsilon
y)}Z_{i0}+\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}Z_{i0}=O\left(\frac{\varepsilon\mu_{i}}{(\mu_{i}+|y-\xi^{\prime}_{i}|)^{3}}\right),\,\,\,\,\
\,\forall\,\,i=1,\ldots,m.$ (3.12)

Proposition 3.3. Let $m$ be a positive integer. Then there exist constants
$C>0$ and $p_{m}>1$ such that for any $p>p_{m}$, any points
$\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$ and any $h\in
C(\overline{\Omega}_{p})$, there is a unique solution $\phi\in
H^{2}(\Omega_{p})$ of problem (3.7) for some coefficients
$c_{ij}\in\mathbb{R}$, $i=1,\ldots,m$, $j=1,J_{i}$, which satisfies

$\displaystyle\|\phi\|_{L^{\infty}(\Omega_{p})}\leq Cp\|h\|_{*}.$ (3.13)

The proof of this result will be split into four steps which we state and
prove next.

Step 1: Constructing a suitable barrier.

Lemma 3.4. There exist constants $R_{1}>0$ and $C>0$, independent of $p$, such
that for any sufficiently large $p$, any points
$\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$ and any $\sigma\in(0,1)$,
there is

$\displaystyle\psi:\,\,\Omega_{p}\setminus\bigcup_{i=1}^{m}B_{R_{1}\mu_{i}}(\xi^{\prime}_{i})\,\,\mapsto\mathbb{R}$

smooth and positive so that

$\displaystyle\mathcal{L}(\psi)=-\Delta_{a(\varepsilon
y)}\psi+\varepsilon^{2}\psi-W_{\xi^{\prime}}\psi$
$\displaystyle\geq\sum_{i=1}^{m}\frac{\mu_{i}^{\sigma}}{|y-\xi^{\prime}_{i}|^{2+\sigma}}+\varepsilon^{2}\,\,\,\quad\,\,\textrm{in}\,\,\,\,\,\,\,\Omega_{p}\setminus\bigcup_{i=1}^{m}B_{R_{1}\mu_{i}}(\xi^{\prime}_{i}),$
$\displaystyle\frac{\partial\psi}{\partial\nu}$ $\displaystyle\geq
0\,\,\,\qquad\qquad\qquad\qquad\quad\,\,\textrm{on}\,\,\,\,\partial\Omega_{p}\setminus\bigcup_{i=1}^{m}B_{R_{1}\mu_{i}}(\xi^{\prime}_{i}),$
$\displaystyle\psi$
$\displaystyle>0\,\,\,\qquad\qquad\qquad\qquad\quad\,\,\textrm{in}\,\,\,\,\,\,\,\Omega_{p}\setminus\bigcup_{i=1}^{m}B_{R_{1}\mu_{i}}(\xi^{\prime}_{i}),$
$\displaystyle\psi$ $\displaystyle\geq
1\,\,\,\qquad\qquad\qquad\qquad\quad\,\,\textrm{on}\,\,\,\,\Omega_{p}\cap\left(\bigcup_{i=1}^{m}\partial
B_{R_{1}\mu_{i}}(\xi^{\prime}_{i})\right).$

Moreover, $\psi$ is bounded uniformly:

$\displaystyle 1<\psi\leq
C\,\,\,\quad\textrm{in}\,\,\,\,\Omega_{p}\setminus\bigcup_{i=1}^{m}B_{R_{1}\mu_{i}}(\xi^{\prime}_{i}).$

###### Proof.

Let us take

$\displaystyle\psi=\sum_{i=1}^{m}\left(1-\frac{\mu_{i}^{\sigma}}{|y-\xi^{\prime}_{i}|^{\sigma}}\right)+C_{1}\Psi_{0}(y),$

where $\Psi_{0}$ is the unique solution of

$\displaystyle\left\\{\begin{aligned} &-\Delta_{a(\varepsilon
y)}\Psi_{0}+\varepsilon^{2}\Psi_{0}=\varepsilon^{2}\,\quad\,\,\textrm{in}\,\,\,\,\Omega_{p},\\\
&\frac{\partial\Psi_{0}}{\partial\nu}=\varepsilon\,\qquad\qquad\qquad\qquad\,\textrm{on}\,\,\,\,\partial\Omega_{p}.\end{aligned}\right.$

Observing that $\Psi_{0}$ is uniformly bounded in $\Omega_{p}$, it is directly
checked that, choosing the positive constant $C_{1}$ larger if necessary,
$\psi$ satisfies all the properties of the lemma for numbers $R_{1}$ and $p$
large enough. ∎

Step 2: Transferring a linear equation. Given $h\in
C^{0,\alpha}(\overline{\Omega}_{p})$ and
$\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$, we first study the linear
equation

$\left\\{\begin{aligned} &\mathcal{L}(\phi)=-\Delta_{a(\varepsilon
y)}\phi+\varepsilon^{2}\phi-W_{\xi^{\prime}}\phi=h\,\,\
\,\,\textrm{in}\,\,\,\,\,\,\Omega_{p},\\\
&\frac{\partial\phi}{\partial\nu}=0\,\,\,\,\,\qquad\qquad\qquad\qquad\qquad\qquad\
\,\textrm{on}\,\,\,\,\partial\Omega_{p}.\end{aligned}\right.$ (3.14)

For the solution of (3.14) satisfying the orthogonality conditions with
respect to $Z_{ij}$, $i=1,\ldots,m$, $j=0,1,J_{i}$, we prove the following a
priori estimate.

Lemma 3.5. There exist $R_{0}>0$ and $p_{m}>1$ such that for any $p>p_{m}$ and
any solution $\phi$ of (3.14) with the orthogonality conditions

$\displaystyle\int_{\Omega_{p}}\chi_{i}Z_{ij}\phi=0\,\,\,\,\,\,\,\,\,\forall\,\,i=1,\ldots,m,\,\,j=0,1,J_{i},$
(3.15)

we have

$\displaystyle\|\phi\|_{L^{\infty}(\Omega_{p})}\leq C\|h\|_{*},$ (3.16)

where $C>0$ is independent of $p$.

###### Proof.

Take $R_{0}=2R_{1}$, $R_{1}$ being the constant of Lemma 3.4. Since
$\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$ and
$\varepsilon\mu_{i}=o(1/p^{\kappa})$ for $p$ large enough, we find
$B_{R_{1}\mu_{i}}(\xi^{\prime}_{i})$ disjointed. Let $h$ be bounded and $\phi$
a solution to (3.14) satisfying (3.15). We first consider the following inner
norm of $\phi$:

$\displaystyle\|\phi\|_{i}=\sup_{y\in\overline{\Omega}_{p}\cap\left(\bigcup_{i=1}^{m}B_{R_{1}\mu_{i}}(\xi^{\prime}_{i})\right)}|\phi(y)|,$

and claim that there is a constant $C>0$ independent of $p$ such that

$\displaystyle\|\phi\|_{L^{\infty}(\Omega_{p})}\leq
C\left(\|\phi\|_{i}+\|h\|_{*}\right).$ (3.17)

Indeed, set

$\displaystyle\widetilde{\phi}(y)=C_{1}\left(\|\phi\|_{i}+\|h\|_{*}\right)\psi(y)\,\qquad\forall\,\,\,y\in\overline{\Omega}_{p}\setminus\bigcup_{i=1}^{m}B_{R_{1}\mu_{i}}(\xi^{\prime}_{i}),$

where $\psi$ is the positive, uniformly bounded barrier constructed by the
previous lemma and the constant $C_{1}>0$ is chosen larger if necessary,
independent of $p$. Then for
$y\in\Omega_{p}\setminus\bigcup_{i=1}^{m}B_{R_{1}\mu_{i}}(\xi^{\prime}_{i})$,

$\displaystyle\mathcal{L}(\widetilde{\phi}\pm\phi)(y)\geq
C_{1}\,\|h\|_{*}\left\\{\sum_{i=1}^{m}\frac{\mu_{i}^{\sigma}}{|y-\xi^{\prime}_{i}|^{2+\sigma}}+\varepsilon^{2}\right\\}\pm
h(y)\geq|h(y)|\pm h(y)\geq 0,$

for
$y\in\partial\Omega_{p}\setminus\bigcup_{i=1}^{m}B_{R_{1}\mu_{i}}(\xi^{\prime}_{i})$,

$\displaystyle\frac{\partial}{\partial\nu}(\widetilde{\phi}\pm\phi)(y)\geq 0,$

and for $y\in\Omega_{p}\cap\left(\bigcup_{i=1}^{m}\partial
B_{R_{1}\mu_{i}}(\xi^{\prime}_{i})\right)$,

$\displaystyle(\widetilde{\phi}\pm\phi)(y)>\|\phi\|_{i}\pm\phi(y)\geq|\phi(y)|\pm\phi(y)\geq
0.$

From the maximum principle (see [46]), it follows that
$-\widetilde{\phi}\leq\phi\leq\widetilde{\phi}$ on
$\overline{\Omega}_{p}\setminus\bigcup_{i=1}^{m}B_{R_{1}\mu_{i}}(\xi^{\prime}_{i})$,
which deduces estimate (3.17).

We prove the lemma by contradiction. Assume that there are sequences of
parameters $p_{n}\rightarrow+\infty$, points
$\xi^{n}=(\xi_{1}^{n},\ldots,\xi_{m}^{n})\in\mathcal{O}_{p_{n}}$, functions
$h_{n}$, and associated solutions $\phi_{n}$ of equation (3.14) with
orthogonality conditions (3.15) such that

$\displaystyle\|\phi_{n}\|_{L^{\infty}(\Omega_{p_{n}})}=1\,\,\quad\,\,\textrm{and}\,\,\quad\,\,\|h_{n}\|_{*}\rightarrow
0,\,\,\quad\,\,\textrm{as}\,\,\,\,n\rightarrow+\infty.$ (3.18)

For each $k\in\\{1,\ldots,l\\}$, we have $\xi_{k}^{n}\in\Omega$ and we
consider
$\widehat{\phi}^{n}_{k}(z)=\phi_{n}\big{(}\mu_{k}^{n}z+(\xi^{n}_{k})^{\prime}\big{)}$,
where $\mu^{n}=(\mu^{n}_{1},\ldots,\mu_{m}^{n})$,
$(\xi^{n}_{k})^{\prime}=\xi^{n}_{k}/\varepsilon_{n}$ and
$\varepsilon_{n}=\exp\left\\{-p_{n}/4\right\\}$. Note that

$\displaystyle
h_{n}(y)=\big{(}-\Delta_{a(\varepsilon_{n}y)}\phi_{n}+\varepsilon_{n}^{2}\phi-
W_{(\xi^{n})^{\prime}}\phi_{n}\big{)}\big{|}_{y=\mu_{k}^{n}z+(\xi^{n}_{k})^{\prime}}=(\mu_{k}^{n})^{-2}\left[-\Delta_{\widehat{a}_{n}}\widehat{\phi}_{k}^{n}+\varepsilon_{n}^{2}(\mu_{k}^{n})^{2}\widehat{\phi}_{k}^{n}-(\mu_{k}^{n})^{2}\widehat{W}^{n}\widehat{\phi}_{k}^{n}\right](z),$

where

$\displaystyle\widehat{a}_{n}(z)=a(\varepsilon_{n}\mu_{k}^{n}z+\xi^{n}_{k}),\qquad\qquad\widehat{W}^{n}(z)=W_{(\xi^{n})^{\prime}}(\mu_{k}^{n}z+(\xi^{n}_{k})^{\prime}).$

By the expansion of $W^{n}$ in (3.2) and elliptic regularity,
$\widehat{\phi}^{n}_{k}$ converges uniformly over compact sets to a bounded
solution $\widehat{\phi}^{\infty}_{k}$ of equation $(\ref{3.3})$, which
satisfies

$\displaystyle\int_{\mathbb{R}^{2}}\chi
Z_{j}\widehat{\phi}_{k}^{\infty}=0\quad\,\,\,\textrm{for}\,\,\,\,j=0,\,1,\,2.$
(3.19)

Thus $\widehat{\phi}^{\infty}_{k}$ is a linear combination of $Z_{j}$,
$j=0,1,2$. Notice that $\int_{\mathbb{R}^{2}}\chi Z_{j}Z_{t}=0$ for $j\neq t$
and $\int_{\mathbb{R}^{2}}\chi Z_{j}^{2}>0$. Hence (3.19) implies
$\widehat{\phi}_{k}^{\infty}\equiv 0$.

As for each $k\in\\{l+1,\ldots,m\\}$, we have $\xi_{k}^{n}\in\partial\Omega$
and we consider
$\widehat{\phi}^{n}_{k}(z)=\phi_{n}\big{(}(A^{n}_{k})^{-1}\mu_{k}^{n}z+(\xi^{n}_{k})^{\prime}\big{)}$,
where $A_{k}^{n}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ is a rotation map
such that
$A_{k}^{n}\nu_{\Omega_{p_{n}}}\big{(}(\xi_{k}^{n})^{\prime}\big{)}=\nu_{\mathbb{R}_{+}^{2}}\big{(}0\big{)}$.
Similar to the above argument, we can get that $\widehat{\phi}^{n}_{k}$
converges uniformly over compact sets to a bounded solution
$\widehat{\phi}^{\infty}_{k}$ of equation $(\ref{3.4})$, which satisfies

$\displaystyle\int_{\mathbb{R}_{+}^{2}}\chi
Z_{j}\widehat{\phi}_{k}^{\infty}=0\quad\,\,\,\textrm{for}\,\,\,\,j=0,\,1.$
(3.20)

Thus $\widehat{\phi}^{\infty}_{k}$ is a linear combination of $Z_{j}$,
$j=0,1$. Notice that $\int_{\mathbb{R}_{+}^{2}}\chi Z_{j}Z_{t}=0$ for $j\neq
t$ and $\int_{\mathbb{R}_{+}^{2}}\chi Z_{j}^{2}>0$. Hence (3.20) implies
$\widehat{\phi}_{k}^{\infty}=0$. Furthermore, we find that
$\lim_{n\rightarrow+\infty}\|\phi_{n}\|_{i}=0$. But (3.17) and (3.18) tell us
$\liminf_{n\rightarrow+\infty}\|\phi_{n}\|_{i}>0$, which is a contradiction. ∎

Step 3: Proving an a priori estimate for solutions to (3.14) that satisfy
orthogonality conditions with respect to $Z_{ij}$, $j=1,J_{i}$ only.

Lemma 3.6. For $p$ large enough, if $\phi$ solves (3.14) and satisfies

$\displaystyle\int_{\Omega_{p}}\chi_{i}Z_{ij}\phi=0\,\,\,\,\,\,\,\,\,\forall\,\,i=1,\ldots,m,\,\,j=1,J_{i},$
(3.21)

then

$\displaystyle\|\phi\|_{L^{\infty}(\Omega_{p})}\leq Cp\,\|h\|_{*},$ (3.22)

where $C>0$ is independent of $p$.

###### Proof.

According to the results in Lemma 3.4 of [17] and Lemma 3.2 of [41], for
simplicity of argument we only consider the validity of estimate (3.22) when
the $m$ concentration points $\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$
satisfy the relation $|\xi_{i}-\xi_{k}|\leq 2d$ for any $i,k=1,\ldots,m$,
$i\neq k$ and for any $d>0$ sufficiently small, fixed and independent of $p$.
Let $R>R_{0}+1$ be a large but fixed number. Denote for $i=1,\ldots,m$,

$\displaystyle\widehat{Z}_{i0}(y)=Z_{i0}(y)-\frac{1}{\mu_{i}}+a_{i0}G(\varepsilon
y,\xi_{i}),$ (3.23)

where

$\displaystyle
a_{i0}=\frac{1}{\mu_{i}\big{[}H(\xi_{i},\xi_{i})-\frac{4}{c_{i}}\log(\varepsilon\mu_{i}R)\big{]}}.$
(3.24)

Note that by estimate (2.26), expansions (2.3) and (2.5), and definitions
(2.11), (3.4), (3.8) and (3.11),

$\displaystyle\frac{1}{8}p\leq-\log(\varepsilon\mu_{i}R)\leq\frac{1}{2}p,$
(3.25)

and

$\displaystyle\widehat{Z}_{i0}(y)=O\left(\frac{\,G(\varepsilon
y,\xi_{i})\,}{p\mu_{i}}\right).$ (3.26)

Let $\eta_{1}$ and $\eta_{2}$ be radial smooth cut-off functions in
$\mathbb{R}^{2}$ such that

$\displaystyle 0\leq\eta_{1}\leq 1;\,\,\,\ \,\,\,|\nabla\eta_{1}|\leq C\,\,\
\textrm{in}\,\,\,\mathbb{R}^{2};\,\,\,\ \,\,\,\eta_{1}\equiv 1\,\,\
\textrm{in}\,\,\,B_{R}(0);\,\,\,\,\ \,\,\,\,\eta_{1}\equiv 0\,\,\
\textrm{in}\,\,\,\mathbb{R}^{2}\setminus B_{R+1}(0);$ $\displaystyle
0\leq\eta_{2}\leq 1;\,\,\,\ \,\,\,|\nabla\eta_{2}|\leq C\,\,\
\textrm{in}\,\,\,\mathbb{R}^{2};\,\,\,\ \,\,\,\eta_{2}\equiv 1\,\,\
\textrm{in}\,\,\,B_{3d}(0);\,\,\,\,\,\,\,\,\,\eta_{2}\equiv 0\,\,\
\textrm{in}\,\,\,\mathbb{R}^{2}\setminus B_{6d}(0).$

We set, for $i\in\\{1,\ldots,l\\}$,

$\displaystyle\eta_{i1}(y)=\eta_{1}\left(\frac{1}{\mu_{i}}\big{|}y-\xi_{i}^{\prime}\big{|}\right),\,\,\quad\quad\,\,\eta_{i2}(y)=\eta_{2}\left(\varepsilon\big{|}y-\xi^{\prime}_{i}\big{|}\right),$
(3.27)

and for $i\in\\{l+1,\ldots,m\\}$,

$\displaystyle\eta_{i1}(y)=\eta_{1}\left(\frac{1}{\mu_{i}}\big{|}F_{i}^{p}(y)\big{|}\right),\,\,\quad\quad\,\,\eta_{i2}(y)=\eta_{2}\left(\varepsilon\big{|}F_{i}^{p}(y)\big{|}\right).$
(3.28)

Now we define the test function

$\displaystyle\widetilde{Z}_{i0}(y)=\eta_{i1}Z_{i0}+(1-\eta_{i1})\eta_{i2}\widehat{Z}_{i0}.$
(3.29)

Given $\phi$ satisfying (3.14) and (3.21), let

$\displaystyle\widetilde{\phi}=\phi+\sum\limits_{i=1}^{m}d_{i}\widetilde{Z}_{i0}+\sum_{i=1}^{m}\sum\limits_{j=1}^{J_{i}}e_{ij}\chi_{i}Z_{ij}.$
(3.30)

We will first prove the existence of $d_{i}$ and $e_{ij}$ such that
$\widetilde{\phi}$ satisfies the orthogonality condition

$\displaystyle\int_{\Omega_{p}}\chi_{i}Z_{ij}\widetilde{\phi}=0\,\,\quad\,\,\forall\,\,i=1,\ldots,m,\,\,j=0,1,J_{i}.$
(3.31)

Multiplying definition (3.30) by $\chi_{i}Z_{ij}$, $i=1,\ldots,m$,
$j=0,1,J_{i}$ and using orthogonality conditions (3.21) and (3.31) and the
fact that $\chi_{i}\chi_{k}\equiv 0$ if $i\neq k$, we get

$\displaystyle
d_{i}\int_{\Omega_{p}}\chi_{i}Z_{i0}\widetilde{Z}_{i0}+\sum_{k\neq
i}^{m}d_{k}\int_{\Omega_{p}}\chi_{i}Z_{i0}\widetilde{Z}_{k0}+\sum\limits_{t=1}^{J_{i}}e_{it}\int_{\Omega_{p}}\chi_{i}^{2}Z_{i0}Z_{it}=-\int_{\Omega_{p}}\chi_{i}Z_{i0}\phi,$
(3.32) $\displaystyle
d_{i}\int_{\Omega_{p}}\chi_{i}Z_{ij}\widetilde{Z}_{i0}+\sum_{k\neq
i}^{m}d_{k}\int_{\Omega_{p}}\chi_{i}Z_{ij}\widetilde{Z}_{k0}+\sum_{t=1}^{J_{i}}e_{it}\int_{\Omega_{p}}\chi^{2}_{i}Z_{ij}Z_{it}=0,\,\quad\
\ \,\ \,j=1,\,J_{i}.$ (3.33)

Remark that for any $i=1,\ldots,l$, $\widetilde{Z}_{i0}$ coincides with
$Z_{i0}$ in $B_{R\mu_{i}}(\xi^{\prime}_{i})$, while for any $i=l+1,\ldots,m$,
$\widetilde{Z}_{i0}$ coincides with $Z_{i0}$ in the region
$\\{y\in\Omega_{p}:|F_{i}^{p}(y)|\leq R\mu_{i}\\}$. Moreover, from definitions
(3.9)-(3.10) we can write $z=\frac{1}{\mu_{i}}F^{p}_{i}(y)$ and its inverse
$y=\big{(}\frac{1}{\mu_{i}}F^{p}_{i}\big{)}^{-1}(z)=\xi^{\prime}_{i}+\frac{1}{\varepsilon}A_{i}^{-1}F^{-1}_{i}(\varepsilon\mu_{i}z)$
such that
$\det\big{(}\nabla\big{(}\frac{1}{\mu_{i}}F_{i}^{p}\big{)}^{-1}(z)\big{)}=\mu_{i}^{2}+O(\varepsilon\mu_{i}^{3}|z|)$
holds in the upper half-ball
$B^{+}_{R\mu_{i}}(0):=B_{R\mu_{i}}(0)\cap\mathbb{R}^{2}_{+}$. Then for any
$i=1,\ldots,l,\,$ $j=1,2\,$ and $\,t=1,2$,

$\displaystyle\int_{\Omega_{p}}\chi_{i}Z_{i0}\widetilde{Z}_{i0}=\int_{\mathbb{R}^{2}}\chi
Z^{2}_{0}=C_{0}>0,\qquad\qquad\qquad\qquad\int_{\Omega_{p}}\chi_{i}^{2}Z_{i0}Z_{it}=\int_{\mathbb{R}^{2}}\chi^{2}Z_{0}Z_{t}=0,$
$\displaystyle\int_{\Omega_{p}}\chi_{i}Z_{ij}\widetilde{Z}_{i0}=\int_{\mathbb{R}^{2}}\chi
Z_{j}Z_{0}=0,\qquad\qquad\qquad\qquad\int_{\Omega_{p}}\chi^{2}_{i}Z_{ij}Z_{it}=\int_{\mathbb{R}^{2}}\chi^{2}Z_{j}Z_{t}=C_{j}\delta_{jt},$

where $\delta_{jt}$ denotes the Kronecker’s symbol, but for any
$i=l+1,\ldots,m$ and $j=t=J_{i}=1$,

$\displaystyle\int_{\Omega_{p}}\chi_{i}Z_{i0}\widetilde{Z}_{i0}=\int_{\mathbb{R}_{+}^{2}}\chi
Z^{2}_{0}[1+O\big{(}\varepsilon\mu_{i}|z|\big{)}]=\frac{C_{0}}{2}+O\left(\varepsilon\mu_{i}\right),\qquad\int_{\Omega_{p}}\chi_{i}^{2}Z_{i0}Z_{i1}=\int_{\mathbb{R}_{+}^{2}}\chi^{2}Z_{0}Z_{1}[1+O\big{(}\varepsilon\mu_{i}|z|\big{)}]=O\left(\varepsilon\mu_{i}\right),$
$\displaystyle\int_{\Omega_{p}}\chi_{i}Z_{i1}\widetilde{Z}_{i0}=\int_{\mathbb{R}_{+}^{2}}\chi
Z_{1}Z_{0}[1+O\big{(}\varepsilon\mu_{i}|z|\big{)}]=O\left(\varepsilon\mu_{i}\right),\,\,\qquad\,\,\int_{\Omega_{p}}\chi^{2}_{i}Z_{i1}^{2}=\int_{\mathbb{R}_{+}^{2}}\chi^{2}Z_{1}^{2}[1+O\big{(}\varepsilon\mu_{i}|z|\big{)}]=\frac{C_{1}}{2}+O\left(\varepsilon\mu_{i}\right).$

Moreover, from (3.26) and (3.29) it follows that for any $i=1,\ldots,m$ and
$j=0,1,J_{i}$,

$\displaystyle\int_{\Omega_{p}}\chi_{i}Z_{ij}\widetilde{Z}_{k0}=O\left(\frac{\mu_{i}\log
p}{p\mu_{k}}\right),\,\quad\,\forall\,\,k\neq i.$

Thus by (3.33),

$\displaystyle
e_{ij}=\left(-d_{i}\int_{\Omega_{p}}\chi_{i}Z_{ij}\widetilde{Z}_{i0}-\sum_{k\neq
i}^{m}d_{k}\int_{\Omega_{p}}\chi_{i}Z_{ij}\widetilde{Z}_{k0}\right)\left/\int_{\Omega_{p}}\chi^{2}_{i}Z^{2}_{ij},\right.\qquad
i=1,\ldots,m,\,\,\,j=1,J_{i}.$ (3.34)

Furthermore,

$\displaystyle|e_{ij}|\leq\left\\{\begin{aligned} &C\sum_{k\neq
i}^{m}|d_{k}|\frac{\mu_{i}\log
p}{p\mu_{k}},\,\quad\qquad\qquad\quad\,\,\,\,\forall\,\,\,i=1,\ldots,l,\,\ \
\,\,\,j=1,2,\\\ &C\varepsilon\mu_{i}|d_{i}|+C\sum_{k\neq
i}^{m}|d_{k}|\frac{\mu_{i}\log
p}{p\mu_{k}},\quad\quad\forall\,\,\,i=l+1,\ldots,m,\,\,j=1.\end{aligned}\right.$
(3.35)

We need just to show that $d_{i}$ is well defined. From (3.32) we can easily
get that for any $i=1,\ldots,l$,

$\displaystyle d_{i}C_{0}+\sum_{k\neq i}^{m}d_{k}O\left(\frac{\mu_{i}\log
p}{p\mu_{k}}\right)=-\int_{\Omega_{p}}\chi_{i}Z_{i0}\phi,$ (3.36)

and for any $i=l+1,\ldots,m$,

$\displaystyle\frac{1}{2}d_{i}C_{0}\big{[}1+O\big{(}\varepsilon\mu_{i}\big{)}\big{]}+\sum_{k\neq
i}^{m}d_{k}O\left(\frac{\mu_{i}\log
p}{p\mu_{k}}\right)+e_{i1}O\big{(}\varepsilon\mu_{i}\big{)}=-\int_{\Omega_{p}}\chi_{i}Z_{i0}\phi,$
(3.37)

where $e_{i1}$ is defined in (3.34) and satisfies

$\displaystyle e_{i1}=d_{i}O\left(\varepsilon\mu_{i}\right)+\sum_{k\neq
i}^{m}d_{k}O\left(\frac{\mu_{i}\log p}{p\mu_{k}}\right).$

We denote $\mathcal{A}$ the coefficient matrix of equations (3.36)-(3.37). By
the above estimates, it is clear that $\mathcal{M}^{-1}\mathcal{A}\mathcal{M}$
is diagonally dominant and thus invertible, where
$\mathcal{M}=\text{diag}(\mu_{1},\ldots,\mu_{m})$. Hence $\mathcal{A}$ is also
invertible and $(d_{1},\ldots,d_{m})$ is well defined.

Estimate (3.22) is a direct consequence of the following two claims.

Claim 1. Let $\mathcal{L}=-\Delta_{a(\varepsilon
y)}+\varepsilon^{2}-W_{\xi^{\prime}}$, then for any $i=1,\ldots,m$ and
$j=1,J_{i}$,

$\displaystyle\big{\|}\mathcal{L}(\chi_{i}Z_{ij})\big{\|}_{*}\leq\frac{C}{\mu_{i}},\,\quad\quad\quad\,\,\,\quad\,\big{\|}\mathcal{L}(\widetilde{Z}_{i0})\big{\|}_{*}\leq
C\frac{\log p}{p\mu_{i}}.$ (3.38)

Claim 2. For any $i=1,\ldots,m$ and $j=1,J_{i}$,

$\displaystyle|d_{i}|\leq
Cp\mu_{i}\|h\|_{*},\,\quad\qquad\quad\,\,\quad|e_{ij}|\leq C\mu_{i}\log
p\,\|h\|_{*}.$ (3.39)

In fact, the definition of $\widetilde{\phi}$ in (3.30) tells us

$\displaystyle\left\\{\begin{aligned}
&\mathcal{L}(\widetilde{\phi})=h+\sum\limits_{i=1}^{m}d_{i}\mathcal{L}(\widetilde{Z}_{i0})+\sum_{i=1}^{m}\sum_{j=1}^{J_{i}}e_{ij}\mathcal{L}(\chi_{i}Z_{ij})\,\quad\,\textrm{in}\,\,\,\
\,\,\Omega_{p},\\\
&\frac{\partial\widetilde{\phi}}{\partial\nu}=0\,\qquad\qquad\qquad\qquad\quad\qquad\,\,\qquad\qquad\qquad\,\,\,\textrm{on}\,\,\,\,\,\partial\Omega_{p}.\end{aligned}\right.$
(3.40)

Then by Lemma 3.5, we obtain

$\displaystyle\|\widetilde{\phi}\|_{L^{\infty}(\Omega_{p})}\leq
C\left\\{\|h\|_{*}+\sum\limits_{i=1}^{m}|d_{i}|\big{\|}\mathcal{L}(\widetilde{Z}_{i0})\big{\|}_{*}+\sum_{i=1}^{m}\sum_{j=1}^{J_{i}}|e_{ij}|\big{\|}\mathcal{L}(\chi_{i}Z_{ij})\big{\|}_{*}\right\\}\leq
C\log p\,\|h\|_{*}.$ (3.41)

Using the definition of $\widetilde{\phi}$ again and the fact that

$\displaystyle\big{\|}\widetilde{Z}_{i0}\big{\|}_{L^{\infty}(\Omega_{p})}\leq\frac{C}{\mu_{i}}\quad\quad\,\textrm{and}\,\quad\quad\big{\|}\chi_{i}Z_{ij}\big{\|}_{L^{\infty}(\Omega_{p})}\leq\frac{C}{\mu_{i}},\,\,\quad\,\forall\,\,\,i=1,\ldots,m,\,\,j=1,J_{i},$
(3.42)

estimate (3.22) then follows from estimate (3.41) and Claim 2.

Proof of Claim 1. Observe that

$\displaystyle\mathcal{L}(\chi_{i}Z_{ij})=\,\chi_{i}\mathcal{L}(Z_{ij})-2\nabla\chi_{i}\nabla
Z_{ij}-Z_{ij}\big{[}\Delta\chi_{i}+\varepsilon\nabla\log a(\varepsilon
y)\nabla\chi_{i}\big{]}.$ (3.43)

For any $i=1,\ldots,l$ and $j=1,2$, we write $z_{i}:=y-\xi^{\prime}_{i}$ and
note that in the region $|z_{i}|=|y-\xi^{\prime}_{i}|\leq\mu_{i}(R_{0}+1)$, by
(3.2), (3.4), (3.5) and (3.8),

$\displaystyle\mathcal{L}(Z_{ij})=$ $\displaystyle\big{(}-\Delta-
W_{\xi^{\prime}}\big{)}\left[\frac{1}{\mu_{i}}Z_{j}\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)\right]-\varepsilon\nabla\log
a(\varepsilon
y)\nabla\left[\frac{1}{\mu_{i}}Z_{j}\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)\right]+\frac{\varepsilon^{2}}{\mu_{i}}Z_{j}\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)$
$\displaystyle=$ $\displaystyle
O\left(\frac{1}{p\mu_{i}}\cdot\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}\right)+O\left(\frac{\varepsilon}{\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2}}\right)+O\left(\frac{\varepsilon^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{1/2}}\right),$

and then, by (3.43),

$\displaystyle\mathcal{L}(\chi_{i}Z_{ij})=O\left(\frac{\mu_{i}p^{-1}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}\right)+O\left(\frac{\mu_{i}^{-1}}{\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2}}\right)+O\left(\frac{\mu_{i}^{-2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{1/2}}\right).$
(3.44)

As for any $i=l+1,\ldots,m$, owing to $F_{i}^{p}(\xi^{\prime}_{i})=(0,0)$ and
$\nabla F_{i}^{p}(\xi^{\prime}_{i})=A_{i}$, we know

$\displaystyle\nabla_{y}=A_{i}\nabla_{z_{i}}+O(\varepsilon|z_{i}|)\nabla_{z_{i}}\,\quad\,\,\,\quad\,\textrm{and}\,\quad\,\,\,\quad\,-\Delta_{y}=-\Delta_{z_{i}}+O(\varepsilon|z_{i}|)\nabla_{z_{i}}^{2}+O(\varepsilon)\nabla_{z_{i}},$
(3.45)

where

$\displaystyle
z_{i}:=F_{i}^{p}(y)=\frac{1}{\varepsilon}F_{i}\big{(}A_{i}(\varepsilon
y-\xi_{i})\big{)}=A_{i}(y-\xi_{i}^{\prime})\big{\\{}1+O\big{(}\varepsilon
A_{i}(y-\xi_{i}^{\prime})\big{)}\big{\\}}.$ (3.46)

In the region $|z_{i}|=|F_{i}^{p}(y)|\leq\mu_{i}(R_{0}+1)$, by (3.2), (3.4),
(3.6), (3.11), (3.45) and (3.46),

$\displaystyle\mathcal{L}(Z_{i1})=$
$\displaystyle\big{(}-\Delta_{z_{i}}-W_{\xi^{\prime}}\big{)}\left[\frac{1}{\mu_{i}}Z_{1}\left(\frac{z_{i}}{\mu_{i}}\right)\right]+\big{(}\varepsilon^{2}+O(\varepsilon|z_{i}|)\nabla_{z_{i}}^{2}+O(\varepsilon)\nabla_{z_{i}}\big{)}\left[\frac{1}{\mu_{i}}Z_{1}\left(\frac{z_{i}}{\mu_{i}}\right)\right]$
$\displaystyle=$ $\displaystyle
O\left(\frac{1}{p\mu_{i}}\cdot\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|z_{i}|^{2})^{2}}\right)+O\left(\frac{\varepsilon}{\mu_{i}^{2}+|z_{i}|^{2}}\right)+O\left(\frac{\varepsilon^{2}}{(\mu_{i}^{2}+|z_{i}|^{2})^{1/2}}\right).$

Thus by (3.43),

$\displaystyle\mathcal{L}(\chi_{i}Z_{i1})=O\left(\frac{\mu_{i}p^{-1}}{(\mu_{i}^{2}+|F_{i}^{p}(y)|^{2})^{2}}\right)+O\left(\frac{\mu_{i}^{-1}}{\mu_{i}^{2}+|F_{i}^{p}(y)|^{2}}\right)+O\left(\frac{\mu_{i}^{-2}}{(\mu_{i}^{2}+|F_{i}^{p}(y)|^{2})^{1/2}}\right).$
(3.47)

Hence by (3.44), (3.46), (3.47) and the definition of $\|\cdot\|_{*}$ in
(2.37), we obtain
$\big{\|}\mathcal{L}(\chi_{i}Z_{ij})\big{\|}_{*}=O\left(1/\mu_{i}\right)$ for
all $i=1,\ldots,m$ and $j=1,J_{i}$.

We now prove the second inequality in (3.38). In fact,

$\displaystyle\mathcal{L}(\widetilde{Z}_{i0})=\eta_{i1}\mathcal{L}(Z_{i0}-\widehat{Z}_{i0})+\eta_{i2}\mathcal{L}(\widehat{Z}_{i0})-(Z_{i0}-\widehat{Z}_{i0})\Delta_{a(\varepsilon
y)}\eta_{i1}-2\nabla\eta_{i1}\nabla(Z_{i0}-\widehat{Z}_{i0})$
$\displaystyle\qquad\qquad\,-2\nabla\eta_{i2}\nabla\widehat{Z}_{i0}-\widehat{Z}_{i0}\Delta_{a(\varepsilon
y)}\eta_{i2}.$ (3.48)

Recalling that $z_{i}=y-\xi^{\prime}_{i}$ for any $i=1,\ldots,l$, but
$z_{i}=F_{i}^{p}(y)$ for any $i=l+1,\ldots,m$, we now consider the four
regions

$\displaystyle\Omega_{1}=\left\\{y\in\Omega_{p}\big{|}\,|z_{i}|\leq\mu_{i}R\right\\},\,\quad\,\quad\quad\quad\quad\quad\quad\quad\,\Omega_{2}=\left\\{y\in\Omega_{p}\big{|}\,\mu_{i}R<|z_{i}|\leq\mu_{i}(R+1)\right\\},$
$\displaystyle\Omega_{3}=\left\\{y\in\Omega_{p}\left|\,\mu_{i}(R+1)<|z_{i}|\leq
3d/\varepsilon\right.\right\\},\,\quad\quad\quad\quad\,\,\,\,\,\,\,\Omega_{4}=\left\\{y\in\Omega_{p}\left|\,3d/\varepsilon<|z_{i}|\leq
6d/\varepsilon\right.\right\\}.$

Notice first that, by (3.4), (3.8), (3.11) and (3.46),

$\displaystyle\left|Z_{i0}-\frac{1}{\mu_{i}}\right|=\frac{2\mu_{i}}{\mu_{i}^{2}+|z_{i}|^{2}}=O\left(\frac{\mu_{i}}{(\mu_{i}+|y-\xi_{i}^{\prime}|)^{2}}\right),$
(3.49)

and for $\mu_{i}R<|z_{i}|\leq 6d/\varepsilon$, by (3.23)-(3.24),

$\displaystyle Z_{i0}-\widehat{Z}_{i0}=\frac{1}{\mu_{i}}-a_{i0}G(\varepsilon
y,\xi_{i})=\frac{1}{\mu_{i}\big{[}H(\xi_{i},\xi_{i})-\frac{4}{c_{i}}\log(\varepsilon\mu_{i}R)\big{]}}\left[\frac{4}{c_{i}}\log\frac{|y-\xi_{i}^{\prime}|}{\mu_{i}R}+O\big{(}\varepsilon^{\alpha}|y-\xi_{i}^{\prime}|^{\alpha}\big{)}\right].$
(3.50)

In $\Omega_{1}$, by (3),

$\displaystyle\mathcal{L}(\widetilde{Z}_{i0})=\mathcal{L}(Z_{i0})=\left[-\Delta_{a(\varepsilon
y)}Z_{i0}-\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}Z_{i0}\right]+\left[\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}-W_{\xi^{\prime}}\right]Z_{i0}+\varepsilon^{2}Z_{i0}.$

Note that, by (3.2),

$\displaystyle\left[\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}-W_{\xi^{\prime}}\right]Z_{i0}=O\left(\frac{\mu_{i}p^{-1}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}\right),\qquad\forall\,\,y\in\Omega_{1}\cup\Omega_{2}.$
(3.51)

Hence in $\Omega_{1}$, by (3.12),

$\displaystyle\mathcal{L}(\widetilde{Z}_{i0})=\mathcal{L}(Z_{i0})=O\left(\frac{1}{p\mu_{i}^{3}}\right).$
(3.52)

In $\Omega_{2}$, by (3.23) and (3),

$\displaystyle\mathcal{L}(\widetilde{Z}_{i0})=\left[-\Delta_{a(\varepsilon
y)}Z_{i0}-\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}Z_{i0}\right]+\left[\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}-W_{\xi^{\prime}}\right]Z_{i0}+\varepsilon^{2}\left(Z_{i0}-\frac{1}{\mu_{i}}\right)$
$\displaystyle+W_{\xi^{\prime}}(1-\eta_{i1})(Z_{i0}-\widehat{Z}_{i0})+\frac{\varepsilon^{2}}{\mu_{i}}\eta_{i1}-2\nabla\eta_{i1}\nabla(Z_{i0}-\widehat{Z}_{i0})-(Z_{i0}-\widehat{Z}_{i0})\Delta_{a(\varepsilon
y)}\eta_{i1}.\qquad\quad$

Note that in $\Omega_{2}$, by (3.25) and (3.50),

$\displaystyle|Z_{i0}-\widehat{Z}_{i0}|=O\left(\frac{1}{p\mu_{i}R}\right)\,\,\quad\quad\,\textrm{and}\,\,\quad\quad\,|\nabla\big{(}Z_{i0}-\widehat{Z}_{i0}\big{)}|=O\left(\frac{1}{p\mu_{i}^{2}R}\right).$
(3.53)

Moreover $|\nabla\eta_{i1}|=O\left(\mu_{i}^{-1}\right)$,
$|\Delta_{a(\varepsilon y)}\eta_{i1}|=O\left(\mu_{i}^{-2}\right)$. Hence in
$\Omega_{2}$, by (3.2), (3.12), (3.49), (3.51) and (3.53),

$\displaystyle\mathcal{L}(\widetilde{Z}_{i0})=O\left(\frac{1}{p\mu_{i}^{3}R}\right).$
(3.54)

In $\Omega_{3}$, by (3.12), (3.23), (3) and (3.49),

$\displaystyle\mathcal{L}(\widetilde{Z}_{i0})=$
$\displaystyle\mathcal{L}(\widehat{Z}_{i0})=\mathcal{L}(Z_{i0})-\mathcal{L}(Z_{i0}-\widehat{Z}_{i0})$
$\displaystyle=$
$\displaystyle\left[\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}-W_{\xi^{\prime}}\right]Z_{i0}+W_{\xi^{\prime}}\left[\frac{1}{\mu_{i}}-a_{i0}G(\varepsilon
y,\xi_{i})\right]+O\left(\frac{\varepsilon\mu_{i}}{(\mu_{i}+|y-\xi^{\prime}_{i}|)^{3}}\right)+O\left(\frac{\varepsilon^{2}\mu_{i}}{(\mu_{i}+|y-\xi_{i}^{\prime}|)^{2}}\right)$
$\displaystyle\equiv$ $\displaystyle
A_{1}+A_{2}+O\left(\frac{\varepsilon\mu_{i}}{(\mu_{i}+|y-\xi^{\prime}_{i}|)^{3}}\right)+O\left(\frac{\varepsilon^{2}\mu_{i}}{(\mu_{i}+|y-\xi_{i}^{\prime}|)^{2}}\right).$

To estimate the first two terms, we need to decompose $\Omega_{3}$ to some
subregions:

$\displaystyle\Omega_{3,i}=\left\\{y\in\Omega_{3}\left|\,\,\mu_{i}(R+1)<|z_{i}|\leq\mu_{i}^{1/2}\big{/}\big{(}3\varepsilon^{1/2}p^{2\kappa}\big{)}\,\right.\right\\},$
$\displaystyle\Omega_{3,k}=$
$\displaystyle\left\\{y\in\Omega_{3}\left|\,\,|z_{k}|\leq\mu_{k}^{1/2}\big{/}\big{(}3\varepsilon^{1/2}p^{2\kappa}\big{)}\,\right.\right\\},\,\,\,\,k\neq
i,\quad\,\,\textrm{and}\quad\,\,\widetilde{\Omega}_{3}=\Omega_{3}\setminus\bigcup_{t=1}^{m}\Omega_{3,t}.$

From (3.1), (3.2) and (3.46) we get

$\displaystyle A_{1}=\left\\{\begin{aligned}
&\frac{8\mu^{2}_{i}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}\left[\frac{1}{p}\left(\omega_{1}-U_{1,0}-\frac{1}{2}U_{1,0}^{2}\right)\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)+O\left(\frac{\log^{4}\big{(}\big{|}\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\big{|}+2\big{)}}{p^{2}}\right)\right]O\left(\frac{1}{\mu_{i}}\right)\,\,\,\textrm{in}\,\,\,\Omega_{3,i},\\\\[5.69054pt]
&O\left(\varepsilon^{2}p^{8\kappa}/\mu_{i}\right)\quad\qquad\qquad\qquad\qquad\qquad\,\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\,\textrm{in}\,\,\,\,\widetilde{\Omega}_{3}.\end{aligned}\right.$

Moreover, by (3.25) and (3.50),

$\displaystyle
A_{2}=\frac{8\mu^{2}_{i}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}O\left(\frac{\log|y-\xi_{i}^{\prime}|-\log\mu_{i}R+\varepsilon^{\alpha}|y-\xi_{i}^{\prime}|^{\alpha}}{p\mu_{i}}\right)\,\,\,\,\,\,\textrm{in}\,\,\,\Omega_{3,i}\cup\widetilde{\Omega}_{3}.$

Then in $\Omega_{3,i}\cup\widetilde{\Omega}_{3}$, by (2.7) and (2.15),

$\displaystyle\mathcal{L}(\widetilde{Z}_{i0})=O\left(\frac{1}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{(3+\sigma)/2}}\cdot\frac{\mu_{i}^{\sigma}}{p}\right)+O\left(\frac{\log|y-\xi^{\prime}_{i}|-\log\mu_{i}R}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}\cdot\frac{\mu_{i}}{p}\right)$
(3.55)

with $\sigma>0$ small but fixed, independent of $p$. In $\Omega_{3,k}$ with
$k\neq i$, by (3.2), (3.12), (3.26) and (3.49),

$\displaystyle\mathcal{L}(\widetilde{Z}_{i0})=\mathcal{L}(\widehat{Z}_{i0})=\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}Z_{i0}-\left[\Delta_{a(\varepsilon
y)}Z_{i0}+\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}Z_{i0}\right]+\varepsilon^{2}\left(Z_{i0}-\frac{1}{\mu_{i}}\right)-W_{\xi^{\prime}}\widehat{Z}_{i0}$
$\displaystyle=O\left(\frac{8\mu^{2}_{k}}{(\mu^{2}_{k}+|y-\xi^{\prime}_{k}|^{2})^{2}}\cdot\frac{\log
p}{p\mu_{i}}\right).\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\qquad\qquad\qquad\qquad\,\,\,\,\,$
(3.56)

In $\Omega_{4}$, by (3),

$\displaystyle\mathcal{L}(\widetilde{Z}_{i0})=$
$\displaystyle\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}\eta_{i2}Z_{i0}-\eta_{i2}\left[\Delta_{a(\varepsilon
y)}Z_{i0}+\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}Z_{i0}\right]+\varepsilon^{2}\eta_{i2}\left(Z_{i0}-\frac{1}{\mu_{i}}\right)$
$\displaystyle-\eta_{i2}W_{\xi^{\prime}}\widehat{Z}_{i0}-2\nabla\eta_{i2}\nabla\widehat{Z}_{i0}-\widehat{Z}_{i0}\Delta_{a(\varepsilon
y)}\eta_{i2}.$

By (2.26) and (3.1) we find $W_{\xi^{\prime}}=O(\varepsilon^{4-\sigma})$ in
$\Omega_{4}$. In addition, $|\nabla\eta_{i2}|=O\left(\varepsilon/d\right)$,
$|\Delta_{a(\varepsilon y)}\eta_{i2}|=O\left(\varepsilon^{2}/d^{2}\right)$,

$\displaystyle|\widehat{Z}_{i0}|=O\left(\frac{|\log
d|}{p\mu_{i}}\right)\,\,\quad\quad\,\textrm{and}\,\,\quad\quad\,|\nabla\widehat{Z}_{i0}|=O\left(\frac{\varepsilon}{pd\mu_{i}}\right).$
(3.57)

Hence by (3.12) and (3.49),

$\displaystyle\mathcal{L}(\widetilde{Z}_{i0})=O\left(\frac{\varepsilon^{2}|\log
d|}{pd^{2}\mu_{i}}\right).$ (3.58)

Combining (2.37), (3.52), (3.54), (3.55), (3) and (3.58), we arrive at

$\displaystyle\big{\|}\mathcal{L}(\widetilde{Z}_{i0})\big{\|}_{*}=O\left(\frac{\log
p}{p\mu_{i}}\right),\,\,\quad\,\forall\,\,i=1,\ldots,m.$

Proof of Claim 2. Multiplying equation (3.40) by $a(\varepsilon
y)\widetilde{Z}_{i0}$, integrating by parts and using relations (3.41)-(3.42),
we get

$\displaystyle\sum_{k=1}^{m}d_{k}$
$\displaystyle\int_{\Omega_{p}}a(\varepsilon
y)\widetilde{Z}_{k0}\mathcal{L}(\widetilde{Z}_{i0})$ $\displaystyle=$
$\displaystyle-\int_{\Omega_{p}}a(\varepsilon
y)h\widetilde{Z}_{i0}+\int_{\Omega_{p}}a(\varepsilon
y)\widetilde{\phi}\mathcal{L}(\widetilde{Z}_{i0})-\sum_{k=1}^{m}\sum_{t=1}^{J_{k}}e_{kt}\int_{\Omega_{p}}a(\varepsilon
y)\chi_{k}Z_{kt}\mathcal{L}(\widetilde{Z}_{i0})$ $\displaystyle\leq$
$\displaystyle
C\frac{\|h\|_{*}}{\mu_{i}}+C\big{\|}\mathcal{L}(\widetilde{Z}_{i0})\big{\|}_{*}\left(\|\widetilde{\phi}\|_{L^{\infty}(\Omega_{p})}+\sum_{k=1}^{m}\sum_{t=1}^{J_{k}}\frac{1}{\mu_{k}}|e_{kt}|\right)$
$\displaystyle\leq$ $\displaystyle
C\frac{\|h\|_{*}}{\mu_{i}}+C\big{\|}\mathcal{L}(\widetilde{Z}_{i0})\big{\|}_{*}\left[\|h\|_{*}+\sum\limits_{k=1}^{m}|d_{k}|\big{\|}\mathcal{L}(\widetilde{Z}_{k0})\big{\|}_{*}+\sum_{k=1}^{m}\sum_{t=1}^{J_{k}}|e_{kt}|\left(\frac{1}{\mu_{k}}+\big{\|}\mathcal{L}(\chi_{k}Z_{kt})\big{\|}_{*}\right)\right],$

where we have applied the following inequality

$\displaystyle\int_{\Omega_{p}}\frac{\mu_{i}^{\sigma}}{(|y-\xi^{\prime}_{i}|+\mu_{i})^{2+\sigma}}\leq
C,\,\,\,\quad\,\,\,\forall\,\,\,i=1,\ldots,m.$

But estimates (3.35) and (3.38) imply that for any $i=1,\ldots,m$,

$\displaystyle|d_{i}|\left|\int_{\Omega_{p}}a(\varepsilon
y)\widetilde{Z}_{i0}\mathcal{L}(\widetilde{Z}_{i0})\right|\leq
C\frac{\|h\|_{*}}{\mu_{i}}+C\sum_{k=1}^{m}\frac{|d_{k}|\log^{2}p}{p^{2}\mu_{i}\mu_{k}}+\sum_{k\neq
i}^{m}\left|d_{k}\int_{\Omega_{p}}a(\varepsilon
y)\widetilde{Z}_{k0}\mathcal{L}(\widetilde{Z}_{i0})\right|.$ (3.59)

From (3.23), (3.29) and (3) we decompose

$\displaystyle\int_{\Omega_{p}}a(\varepsilon
y)\widetilde{Z}_{i0}\mathcal{L}(\widetilde{Z}_{i0})=J+K,$ (3.60)

where

$\displaystyle J=$ $\displaystyle\int_{\Omega_{p}}a(\varepsilon
y)\widetilde{Z}_{i0}\left[\eta_{i1}\mathcal{L}(Z_{i0}-\widehat{Z}_{i0})+\eta_{i2}\mathcal{L}(\widehat{Z}_{i0})\right]$
$\displaystyle=$ $\displaystyle\int_{\Omega_{p}}a(\varepsilon
y)\eta_{i2}^{2}\left\\{Z_{i0}-\big{(}1-\eta_{i1}\big{)}\left[\frac{1}{\mu_{i}}-a_{i0}G(\varepsilon
y,\xi_{i})\right]\right\\}\times\left\\{W_{\xi^{\prime}}\big{(}1-\eta_{i1}\big{)}\left[\frac{1}{\mu_{i}}-a_{i0}G(\varepsilon
y,\xi_{i})\right]\right.$ $\displaystyle\left.-\left[\Delta_{a(\varepsilon
y)}Z_{i0}+\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}Z_{i0}\right]+\left[\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}-W_{\xi^{\prime}}\right]Z_{i0}+\varepsilon^{2}\left(Z_{i0}-\frac{1}{\mu_{i}}\right)+\frac{\varepsilon^{2}}{\mu_{i}}\eta_{i1}\right\\},$

and

$\displaystyle K=\int_{\Omega_{p}}a(\varepsilon
y)\widetilde{Z}_{i0}\left[-(Z_{i0}-\widehat{Z}_{i0})\Delta_{a(\varepsilon
y)}\eta_{i1}-2\nabla\eta_{i1}\nabla(Z_{i0}-\widehat{Z}_{i0})-2\nabla\eta_{i2}\nabla\widehat{Z}_{i0}-\widehat{Z}_{i0}\Delta_{a(\varepsilon
y)}\eta_{i2}\right].$

Let us analyze first the behavior of $J$. By (3.1), (3.2), (3.8), (3.11),
(3.25) and (3.50),

$\displaystyle\int_{|z_{i}|\leq\mu_{i}^{1/2}\big{/}\big{(}3\varepsilon^{1/2}p^{2\kappa}\big{)}}a(\varepsilon
y)\eta_{i2}^{2}\left\\{Z_{i0}-(1-\eta_{i1})\left[\frac{1}{\mu_{i}}-a_{i0}G(\varepsilon
y,\xi_{i})\right]\right\\}\times
W_{\xi^{\prime}}\big{(}1-\eta_{i1}\big{)}\left[\frac{1}{\mu_{i}}-a_{i0}G(\varepsilon
y,\xi_{i})\right]dy$ $\displaystyle=$ $\displaystyle
O\left(\frac{1}{p\mu_{i}^{2}R}\right).$

By (3.8), (3.11), (3.12), (3.25) and (3.50),

$\displaystyle\int_{|z_{i}|\leq\mu_{i}^{1/2}\big{/}\big{(}3\varepsilon^{1/2}p^{2\kappa}\big{)}}\,a(\varepsilon
y)\eta_{i2}^{2}\left\\{Z_{i0}-(1-\eta_{i1})\left[\frac{1}{\mu_{i}}-a_{i0}G(\varepsilon
y,\xi_{i})\right]\right\\}\times\left[\Delta_{a(\varepsilon
y)}Z_{i0}+\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}Z_{i0}\right]dy$
$\displaystyle=$ $\displaystyle O\left(\frac{\varepsilon}{\mu_{i}}\right).$

By (3.8), (3.11), (3.25), (3.49) and (3.50),

$\displaystyle\int_{|z_{i}|\leq\mu_{i}^{1/2}\big{/}\big{(}3\varepsilon^{1/2}p^{2\kappa}\big{)}}\,a(\varepsilon
y)\eta_{i2}^{2}\left\\{Z_{i0}-(1-\eta_{i1})\left[\frac{1}{\mu_{i}}-a_{i0}G(\varepsilon
y,\xi_{i})\right]\right\\}\times\left[\varepsilon^{2}\left(Z_{i0}-\frac{1}{\mu_{i}}\right)+\frac{\varepsilon^{2}}{\mu_{i}}\eta_{i1}\right]dy$
$\displaystyle=$ $\displaystyle O\left(p\varepsilon^{2}\right).$

By (3.2), (3.8), (3.11), (3.25) and (3.50),

$\displaystyle\int_{|z_{i}|\leq\mu_{i}^{1/2}\big{/}\big{(}3\varepsilon^{1/2}p^{2\kappa}\big{)}}\,a(\varepsilon
y)\eta_{i2}^{2}(1-\eta_{i1})\left[\frac{1}{\mu_{i}}-a_{i0}G(\varepsilon
y,\xi_{i})\right]\left[\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}-W_{\xi^{\prime}}\right]Z_{i0}dy$
$\displaystyle=$
$\displaystyle\int_{\mu_{i}R<|z_{i}|\leq\mu_{i}^{1/2}\big{/}\big{(}3\varepsilon^{1/2}p^{2\kappa}\big{)}}\,\frac{8\mu^{2}_{i}a(\varepsilon
y)}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}\left[\,\frac{1}{p}\left(\omega_{1}-U_{1,0}-\frac{1}{2}U_{1,0}^{2}\right)\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)+O\left(\frac{1}{p^{2}}\log^{4}\left|\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right|\right)\right]$
$\displaystyle\,\,\,\,\times
O\left(\frac{\log|y-\xi_{i}^{\prime}|-\log(\mu_{i}R)+\varepsilon^{\alpha}|y-\xi_{i}^{\prime}|^{\alpha}}{p\mu_{i}^{2}}\right)dy$
$\displaystyle=$ $\displaystyle\,O\left(\frac{1}{p^{2}\mu_{i}^{2}R}\right).$

While by (3.1), (3.2), (3.26), (3.29) and (3),

$\displaystyle\int_{|z_{i}|>\mu_{i}^{1/2}\big{/}\big{(}3\varepsilon^{1/2}p^{2\kappa}\big{)}}\,a(\varepsilon
y)\widetilde{Z}_{i0}\left[\eta_{i1}\mathcal{L}(Z_{i0}-\widehat{Z}_{i0})+\eta_{i2}\mathcal{L}(\widehat{Z}_{i0})\right]dy$
$\displaystyle=$ $\displaystyle\sum_{k\neq
i}^{m}\int_{\Omega_{3,k}}a(\varepsilon
y)\eta_{i2}^{2}\widehat{Z}_{i0}\mathcal{L}(\widehat{Z}_{i0})dy+\int_{\widetilde{\Omega}_{3}\cup\Omega_{4}}a(\varepsilon
y)\eta_{i2}^{2}\widehat{Z}_{i0}\mathcal{L}(\widehat{Z}_{i0})dy$
$\displaystyle=$ $\displaystyle\sum_{k\neq
i}^{m}O\left(\int_{0}^{\mu_{k}^{1/2}\big{/}\big{(}3\varepsilon^{1/2}p^{2\kappa}\big{)}}\frac{8\mu^{2}_{k}}{(\mu_{k}^{2}+r^{2})^{2}}\frac{\log^{2}p}{p^{2}\mu_{i}^{2}}rdr\right)+O\left(\int_{\mu_{i}^{1/2}\big{/}\big{(}3\varepsilon^{1/2}p^{2\kappa}\big{)}}^{6d/\varepsilon}\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+r^{2})^{2}}\frac{\log\varepsilon
r}{p\mu_{i}^{2}}rdr\right)$ $\displaystyle=$
$\displaystyle\,O\left(\frac{\log^{2}p}{p^{2}\mu^{2}_{i}}\right).$

Thus

$\displaystyle
J=\int_{|z_{i}|\leq\mu_{i}^{1/2}\big{/}\big{(}3\varepsilon^{1/2}p^{2\kappa}\big{)}}\,a(\varepsilon
y)\eta_{i2}^{2}Z_{i0}^{2}\left[\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}-W_{\xi^{\prime}}\right]dy+O\left(\frac{1}{p\mu_{i}^{2}R}\right).$

Note that in a straightforward but tedious way, by (2) we can compute

$\displaystyle\int_{\mathbb{R}^{2}}e^{U_{1,0}(z)}Z_{0}^{2}(z)\left(\omega_{1}-U_{1,0}-\frac{1}{2}U_{1,0}^{2}\right)(z)dz=-8\pi.$

Hence by (2.21), (3.2), (3.8), (3.11) and (3.46),

$\displaystyle
J=\frac{c_{i}a(\xi_{i})}{p\mu_{i}^{2}}\left[1+O\left(\frac{1}{R}\right)\right],\quad\forall\,\,i=1,\ldots,m.$
(3.61)

Next, we estimate $K$. Integrating by parts the first term and the last term
of $K$ respectively, we get

$\displaystyle K=$ $\displaystyle-\int_{\Omega_{2}}a(\varepsilon
y)Z_{i0}\nabla\eta_{i1}\nabla(Z_{i0}-\widehat{Z}_{i0})+\int_{\Omega_{2}}a(\varepsilon
y)\big{(}Z_{i0}-\widehat{Z}_{i0}\big{)}\nabla\eta_{i1}\nabla(Z_{i0}-\widehat{Z}_{i0})$
$\displaystyle+\int_{\Omega_{2}}a(\varepsilon
y)(Z_{i0}-\widehat{Z}_{i0})^{2}|\nabla\eta_{i1}|^{2}+\int_{\Omega_{2}}a(\varepsilon
y)(Z_{i0}-\widehat{Z}_{i0})\nabla\eta_{i1}\nabla\widehat{Z}_{i0}+\int_{\Omega_{4}}a(\varepsilon
y)|\widehat{Z}_{i0}|^{2}|\nabla\eta_{i2}|^{2}$ $\displaystyle\equiv$
$\displaystyle K_{21}+K_{22}+K_{23}+K_{24}+K_{4}.$

From (2.21), (3.4), (3.8), (3.11), (3.23), (3.24), (3.45), (3.46) and (3.50)
we can conclude that

$\displaystyle K_{21}=$
$\displaystyle-\frac{a_{i0}}{\mu_{i}^{2}}\int_{\\{\mu_{i}R<|z_{i}|\leq\mu_{i}(R+1)\\}}\frac{1}{|y-\xi^{\prime}_{i}|}a(\varepsilon
y)Z_{0}\left(\frac{z_{i}}{\mu_{i}}\right)\eta_{1}^{\prime}\left(\frac{|z_{i}|}{\mu_{i}}\right)\left[\frac{4}{c_{i}}+o(1)\right]dy$
$\displaystyle=$
$\displaystyle-\frac{c_{i}a_{i0}}{4\mu_{i}}\int_{R}^{R+1}a(\xi_{i})\eta_{1}^{\prime}(r)\left[\frac{4}{c_{i}}+O\left(\frac{1}{r^{2}}\right)\right]dr$
$\displaystyle=$
$\displaystyle\frac{c_{i}a(\xi_{i})}{p\mu_{i}^{2}}\left[1+O\left(\frac{1}{R^{2}}\right)\right].$

Moreover, by (3.4), (3.8), (3.11), (3.45) and (3.53) we find
$|\nabla\widehat{Z}_{i0}|=O\big{(}\frac{1}{\mu_{i}^{2}R^{3}}\big{)}$ in
$\Omega_{2}$. Furthermore,

$\displaystyle
K_{22}=O\left(\frac{1}{p^{2}\mu_{i}^{2}R}\right),\,\,\quad\quad\quad\,\,K_{23}=O\left(\frac{1}{p^{2}\mu_{i}^{2}R}\right)\,\,\quad\quad\quad\,\,\textrm{and}\,\,\quad\quad\quad\,\,K_{24}=O\left(\frac{1}{p\mu^{2}_{i}R^{3}}\right).$

By (3.57),

$\displaystyle K_{4}=O\left(\frac{|\log d|^{2}}{p^{2}\mu_{i}^{2}}\right).$

Combining all these estimates, we conclude

$\displaystyle\int_{\Omega_{p}}a(\varepsilon
y)\widetilde{Z}_{i0}\mathcal{L}(\widetilde{Z}_{i0})=J+K=\frac{2c_{i}a(\xi_{i})}{p\mu_{i}^{2}}\left[1+O\left(\frac{1}{R}\right)\right],\quad\forall\,\,i=1,\ldots,m.$
(3.62)

According to (3.59), we need just to consider $\int_{\Omega_{p}}a(\varepsilon
y)\widetilde{Z}_{k0}\mathcal{L}(\widetilde{Z}_{i0})$ when $k\neq i$. Using the
previous estimates of $\mathcal{L}(\widetilde{Z}_{i0})$ and
$\widetilde{Z}_{k0}$, we can easily prove that

$\displaystyle\int_{\Omega_{1}}a(\varepsilon
y)\widetilde{Z}_{k0}\mathcal{L}(\widetilde{Z}_{i0})=O\left(\frac{R^{2}\log
p}{p^{2}\mu_{i}\mu_{k}}\right),\,\qquad\quad\quad\,\,\int_{\Omega_{2}}a(\varepsilon
y)\widetilde{Z}_{k0}\mathcal{L}(\widetilde{Z}_{i0})=O\left(\frac{\log
p}{p^{2}\mu_{i}\mu_{k}}\right),$ $\displaystyle\int_{\Omega_{4}}a(\varepsilon
y)\widetilde{Z}_{k0}\mathcal{L}(\widetilde{Z}_{i0})=O\left(\frac{|\log
d|^{2}}{p^{2}\mu_{i}\mu_{k}}\right),\,\,\,\,\quad\qquad\,\,\,\,\int_{\Omega_{3,i}\cup\widetilde{\Omega}_{3}}a(\varepsilon
y)\widetilde{Z}_{k0}\mathcal{L}(\widetilde{Z}_{i0})=O\left(\frac{\log
p}{p^{2}\mu_{i}\mu_{k}}\right),$

and

$\displaystyle\int_{\Omega_{3,l}}a(\varepsilon
y)\widetilde{Z}_{k0}\mathcal{L}(\widetilde{Z}_{i0})=O\left(\frac{\log^{2}p}{p^{2}\mu_{i}\mu_{k}}\right)\quad\quad\textrm{for
all}\,\,\,l\neq i\,\,\,\textrm{and}\,\,\,l\neq k.$

It remains to calculate the integral over $\Omega_{3,k}$. From (3.29) and an
integration by parts we get

$\displaystyle\int_{\Omega_{3,k}}a(\varepsilon
y)\widetilde{Z}_{k0}\mathcal{L}(\widetilde{Z}_{i0})=\int_{\Omega_{3,k}}a(\varepsilon
y)\widetilde{Z}_{i0}\mathcal{L}(\widetilde{Z}_{k0})-\int_{\partial\Omega_{3,k}}a(\varepsilon
y)\widehat{Z}_{k0}\frac{\partial\widehat{Z}_{i0}}{\partial\nu}+\int_{\partial\Omega_{3,k}}a(\varepsilon
y)\widehat{Z}_{i0}\frac{\partial\widehat{Z}_{k0}}{\partial\nu}.$

As above, we know

$\displaystyle\int_{\Omega_{3,k}}a(\varepsilon
y)\widetilde{Z}_{i0}\mathcal{L}(\widetilde{Z}_{k0})=O\left(\frac{\log
p}{p^{2}\mu_{i}\mu_{k}}\right).$

On $\partial\Omega_{3,k}$, by (2.9) and (3.26) we have

$\displaystyle\,\,\,\,\,\,\,\widehat{Z}_{k0}=O\left(\frac{1}{\mu_{k}}\right),\,\,\qquad\,\qquad\qquad\,\qquad\,\,\widehat{Z}_{i0}=O\left(\frac{\log
p}{p\mu_{i}}\right),$

and

$\displaystyle|\nabla\widehat{Z}_{k0}|=O\left(\frac{\varepsilon^{1/2}p^{2\kappa-1}}{\mu_{k}^{3/2}}\right),\,\qquad\,\,\,\qquad\,|\nabla\widehat{Z}_{i0}|=O\left(\frac{\,\varepsilon
p^{\kappa-1}\,}{\mu_{i}}\right).$

So

$\displaystyle\int_{\Omega_{3,k}}a(\varepsilon
y)\widetilde{Z}_{k0}\mathcal{L}(\widetilde{Z}_{i0})=O\left(\frac{\log
p}{p^{2}\mu_{i}\mu_{k}}\right).$

By the above estimates, we obtain

$\displaystyle\int_{\Omega_{p}}a(\varepsilon
y)\widetilde{Z}_{k0}\mathcal{L}(\widetilde{Z}_{i0})=O\left(\frac{\log^{2}p}{p^{2}\mu_{i}\mu_{k}}\right),\,\,\quad\,\,\textrm{if}\,\,\,\,i\neq
k.$ (3.63)

As a consequence, replacing (3.62) and (3.63) to (3.59) we get

$\displaystyle\frac{|d_{i}|}{\mu_{i}}\leq
Cp\|h\|_{*}+C\sum_{k=1}^{m}\frac{|d_{k}|}{\mu_{k}}\frac{\log^{2}p}{p}.$ (3.64)

Using linear algebra arguments, we then prove Claim 2 for $d_{i}$ and complete
the proof by (3.35). ∎

Step 4: Proof of Proposition 3.3. We first establish the validity of the a
priori estimate

$\displaystyle\|\phi\|_{L^{\infty}(\Omega_{p})}\leq Cp\|h\|_{*}$ (3.65)

for any $\phi$, $c_{ij}$ solutions of problem (3.7) and any $h\in
C^{0,\alpha}(\overline{\Omega}_{p})$. Step 3 gives

$\displaystyle\|\phi\|_{L^{\infty}(\Omega_{p})}\leq
Cp\left(\|h\|_{*}+\sum_{i=1}^{m}\sum_{j=1}^{J_{i}}|c_{ij}|\cdot\|\chi_{i}Z_{ij}\|_{*}\right)\leq
Cp\left(\|h\|_{*}+\sum_{i=1}^{m}\sum_{j=1}^{J_{i}}\mu_{i}|c_{ij}|\right).$

As before, arguing by contradiction of (3.65), we can proceed as in Step 2 and
suppose further that

$\displaystyle\|\phi_{n}\|_{L^{\infty}(\Omega_{p_{n}})}=1,\,\,\
\,\quad\,\,p_{n}\|h_{n}\|_{*}\rightarrow 0,\,\,\
\,\quad\,\,p_{n}\sum_{i=1}^{m}\sum_{j=1}^{J_{i}}\mu_{i}^{n}|c_{ij}^{n}|\geq\delta>0\,\,\,\quad\,\,\,\textrm{as}\,\,\,n\rightarrow+\infty.$
(3.66)

We omit the dependence on $n$. It suffices to estimate the values of the
constants $c_{ij}$. For this aim, let us consider the cut-off function
$\eta_{i2}$ given by (3.27)-(3.28). For any $i=1,\ldots,m$ and $j=1,J_{i}$,
multiplying (3.7) by $a(\varepsilon y)\eta_{i2}Z_{ij}$ and integrating by
parts we find

$\displaystyle\int_{\Omega_{p}}a(\varepsilon
y)\phi\mathcal{L}(\eta_{i2}Z_{ij})=\int_{\Omega_{p}}a(\varepsilon
y)h\eta_{i2}Z_{ij}+\sum_{k=1}^{m}\sum_{t=1}^{J_{k}}c_{kt}\int_{\Omega_{p}}\chi_{k}Z_{kt}\eta_{i2}Z_{ij}.$
(3.67)

Notice that

$\displaystyle\mathcal{L}(\eta_{i2}Z_{ij})=$
$\displaystyle\eta_{i2}\mathcal{L}(Z_{ij})-Z_{ij}\Delta_{a(\varepsilon
y)}\eta_{i2}-2\nabla\eta_{i2}\nabla Z_{ij}$ $\displaystyle=$
$\displaystyle\eta_{i2}\left[-\Delta_{a(\varepsilon
y)}Z_{ij}-\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}Z_{ij}\right]+\left[\frac{8\mu_{i}^{2}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}-W_{\xi^{\prime}}\right]\eta_{i2}Z_{ij}-\varepsilon^{2}\eta_{i2}Z_{ij}+O\left(\varepsilon^{3}/d^{3}\right)$
$\displaystyle\equiv$ $\displaystyle
B_{1}+B_{2}-\varepsilon^{2}\eta_{i2}Z_{ij}+O\left(\varepsilon^{3}/d^{3}\right).$

By (3.4), (3.8), (3.11), (3.45) and (3.46) we can compute

$\displaystyle
B_{1}=O\left(\frac{\varepsilon}{(|y-\xi^{\prime}_{i}|+\mu_{i})^{2}}\right).$

To estimate $B_{2}$, we decompose $\text{supp}(\eta_{i2})$ into several
regions:

$\displaystyle\widehat{\Omega}_{1k}=\left\\{y\in\Omega_{p}\big{|}\,|z_{k}|\leq\frac{\mu_{k}^{1/2}}{3\varepsilon^{1/2}p^{2\kappa}}\right\\},\,\,\
\,\,\forall\,\,k=1,\ldots,m,$
$\displaystyle\widehat{\Omega}_{2}=\left\\{y\in\Omega_{p}\big{|}\,|z_{i}|\leq\frac{6d}{\varepsilon},\,\quad\,|z_{k}|>\frac{\mu_{k}^{1/2}}{3\varepsilon^{1/2}p^{2\kappa}},\,\,\,\,\,\,\,\,\,k=1,\ldots,m\right\\},$

where $z_{k}=y-\xi^{\prime}_{k}$ for $k=1,\ldots,l$, but $z_{k}=F_{k}^{p}(y)$
for $k=l+1,\ldots,m$. Note that, by (2.9) and (3.46),

$\displaystyle|y-\xi^{\prime}_{i}|\geq|\xi_{i}^{\prime}-\xi_{k}^{\prime}|-|y-\xi_{k}^{\prime}|\geq|\xi_{i}^{\prime}-\xi_{k}^{\prime}|-\frac{\mu_{k}^{1/2}}{\varepsilon^{1/2}p^{2\kappa}}\geq\frac{1}{\varepsilon
p^{\kappa}}\left[1-\frac{(\varepsilon\mu_{k})^{1/2}}{p^{\kappa}}\right],$
(3.68)

uniformly in $\widehat{\Omega}_{1k}$, $k\neq i$. In $\widehat{\Omega}_{1i}$,
by (3.2), (3.8) and (3.11) we have that for any $i=1,\ldots,l$ and $j=1,2$,

$\displaystyle
B_{2}=-\frac{8\mu^{2}_{i}(y-\xi^{\prime}_{i})_{j}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{3}}\left\\{\frac{1}{p}\left(\omega_{1}-U_{1,0}-\frac{1}{2}U_{1,0}^{2}\right)\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)+O\left(\frac{\log^{4}\big{(}\big{|}\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\big{|}+2\big{)}}{p^{2}}\right)\right\\},$

and for any $i=l+1,\ldots,m$ and $j=J_{i}=1$,

$\displaystyle
B_{2}=-\frac{8\mu^{2}_{i}}{(\mu_{i}^{2}+|y-\xi^{\prime}_{i}|^{2})^{2}}\frac{(F_{i}^{p}(y))_{1}}{\mu_{i}^{2}+|F_{i}^{p}(y)|^{2}}\left\\{\frac{1}{p}\left(\omega_{1}-U_{1,0}-\frac{1}{2}U_{1,0}^{2}\right)\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)+O\left(\frac{\log^{4}\big{(}\big{|}\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\big{|}+2\big{)}}{p^{2}}\right)\right\\}.$

In $\widehat{\Omega}_{1k}$, $k\neq i$, by (3.2), (3.8), (3.11), (3.46) and
(3.68),

$\displaystyle B_{2}=O\left(\frac{8\varepsilon
p^{\kappa}\mu^{2}_{k}}{(\mu_{k}^{2}+|y-\xi^{\prime}_{k}|^{2})^{2}}\right).$

In $\widehat{\Omega}_{2}$, by (3.1),

$\displaystyle
B_{2}=\sum_{k=1}^{m}O\left(\frac{\mu_{k}^{2}}{|y-\xi^{\prime}_{k}|^{4}}\cdot\frac{\varepsilon^{1/2}p^{2\kappa}}{\mu_{i}^{1/2}}\right).$

Hence by (2.7), (2.15) and (3.46),

$\displaystyle\int_{\Omega_{p}}a(\varepsilon
y)\phi\mathcal{L}(\eta_{i2}Z_{ij})=-\frac{1}{p\mu_{i}}a(\xi_{i})E_{j}(\widehat{\phi}_{i})+O\left(\frac{1}{p^{2}\mu_{i}}\|\phi\|_{L^{\infty}(\Omega_{p})}\right),$
(3.69)

where for any $i=1,\ldots,l$ and $j=1,2$,
$\widehat{\phi}_{i}(z)=\phi\big{(}\xi_{i}^{\prime}+\mu_{i}z\big{)}$ and

$\displaystyle
E_{j}(\widehat{\phi}_{i})=\int_{B_{\frac{1}{3\varepsilon^{1/2}\mu_{i}^{1/2}p^{2\kappa}}}\left(0\right)}\frac{8z_{j}}{(|z|^{2}+1)^{3}}\left(\omega_{1}-U_{1,0}-\frac{1}{2}U_{1,0}^{2}\right)(|z|)\widehat{\phi}_{i}(z)dz,$

but for any $i=l+1,\ldots,m$ and $j=1$,
$\widehat{\phi}_{i}(z)=\phi\big{(}(F_{i}^{p})^{-1}(\mu_{i}z)\big{)}$ and

$\displaystyle E_{j}(\widehat{\phi}_{i})=\int_{\mathbb{R}_{+}^{2}\bigcap
B_{\frac{1}{3\varepsilon^{1/2}\mu_{i}^{1/2}p^{2\kappa}}}\left(0\right)}\frac{8z_{j}}{(|z|^{2}+1)^{3}}\left(\omega_{1}-U_{1,0}-\frac{1}{2}U_{1,0}^{2}\right)(|z|)\widehat{\phi}_{i}(z)dz.$

On the other hand, since $\|\eta_{i2}Z_{ij}\|_{L^{\infty}(\Omega_{p})}\leq
C\mu_{i}^{-1}$, we obtain

$\displaystyle\int_{\Omega_{p}}a(\varepsilon
y)h\eta_{i2}Z_{ij}=O\left(\frac{1}{\mu_{i}}\|h\|_{*}\right).$ (3.70)

Moreover, if $1\leq k=i\leq l$, by (3.4) and (3.8),

$\displaystyle\int_{\Omega_{p}}\chi_{k}Z_{kt}\eta_{i2}Z_{ij}=\int_{\mathbb{R}^{2}}\chi
Z_{t}Z_{j}dz=D_{t}\delta_{tj},$ (3.71)

and if $l+1\leq k=i\leq m$, by (3.4), (3.10) and (3.11),

$\displaystyle\int_{\Omega_{p}}\chi_{k}Z_{k1}\eta_{i2}Z_{i1}=\int_{\mathbb{R}_{+}^{2}}\chi
Z^{2}_{1}\big{[}1+O\left(\varepsilon\mu_{i}|z|\right)\big{]}dz=\frac{1}{2}D_{1}\big{[}1+O\left(\varepsilon\mu_{i}\right)\big{]},$
(3.72)

and if $k\neq i$, by (3.68),

$\displaystyle\int_{\Omega_{p}}\chi_{k}Z_{kt}\eta_{i2}Z_{ij}=O\left(\mu_{k}\varepsilon
p^{\kappa}\right).$ (3.73)

Inserting estimates (3.69)-(3.73) into (3.67), we deduce that for any
$i=1,\ldots,m$ and $j=1,J_{i}$,

$\displaystyle D_{j}c_{ij}+O\left(\sum\limits_{k\neq
i}^{m}\sum\limits_{t=1}^{J_{k}}\mu_{k}\varepsilon
p^{\kappa}|c_{kt}|\right)=O\left(\frac{1}{\mu_{i}}\|h\|_{*}+\frac{1}{p\mu_{i}}\|\phi\|_{L^{\infty}(\Omega_{p})}\right).$

Furthermore,

$\displaystyle\sum_{i=1}^{m}\sum_{j=1}^{J_{i}}\mu_{i}|c_{ij}|=O\left(\|h\|_{*}+\frac{1}{p}\|\phi\|_{L^{\infty}(\Omega_{p})}\right).$
(3.74)

Since
$\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{J_{i}}\mu_{i}|c_{ij}|=o\left(1\right)$,
as in contradiction arguments of Step 2, we conclude that for any
$i=1,\ldots,l$,

$\displaystyle\widehat{\phi}_{i}\rightarrow
C_{i}\frac{|z|^{2}-1}{|z|^{2}+1}\,\,\ \,\,\,\,\,\textrm{uniformly
in}\,\,\,C_{loc}^{0}(\mathbb{R}^{2}),$

but for any $i=l+1,\ldots,m$,

$\displaystyle\widehat{\phi}_{i}\rightarrow
C_{i}\frac{|z|^{2}-1}{|z|^{2}+1}\,\,\ \,\,\,\,\,\textrm{uniformly
in}\,\,\,C_{loc}^{0}(\mathbb{R}_{+}^{2}),$

with some constant $C_{i}\in\mathbb{R}$. Hence in (3.69), we have a better
estimate, since by Lebesgue’s theorem we can derive that for any
$i=1,\ldots,l$ and $j=1,2$,

$\displaystyle E_{j}(\widehat{\phi}_{i})\longrightarrow
C_{i}\int_{\mathbb{R}^{2}}\frac{8z_{j}}{(|z|^{2}+1)^{3}}\frac{|z|^{2}-1}{|z|^{2}+1}\left(\omega_{1}-U_{1,0}-\frac{1}{2}U_{1,0}^{2}\right)(|z|)dz=0,$

and for any $i=l+1,\ldots,m$ and $j=1$,

$\displaystyle E_{j}(\widehat{\phi}_{i})\longrightarrow
C_{i}\int_{\mathbb{R}^{2}_{+}}\frac{8z_{j}}{(|z|^{2}+1)^{3}}\frac{|z|^{2}-1}{|z|^{2}+1}\left(\omega_{1}-U_{1,0}-\frac{1}{2}U_{1,0}^{2}\right)(|z|)dz=0.$

Therefore,

$\displaystyle\sum_{i=1}^{m}\sum_{j=1}^{J_{i}}\mu_{i}|c_{ij}|=o(\frac{1}{p})+O(\|h\|_{*}),$

which is impossible because of (3.66). So estimate (3.65) is established and
then by (3.74), we obtain

$\displaystyle|c_{ij}|\leq C\frac{1}{\mu_{i}}\|h\|_{*}.$ (3.75)

Now consider the Hilbert space

$\displaystyle H_{\xi}=\left\\{\phi\in
H^{1}(\Omega_{p})\left|\,\,\int_{\Omega_{p}}\chi_{i}Z_{ij}\phi=0\quad\textrm{for
any}\,\,\,i=1,\ldots,m,\,\,j=1,J_{i};\quad\frac{\partial\phi}{\partial\nu}=0\quad\textrm{on}\,\,\,\partial\Omega_{p}\right.\right\\}$

with the norm $\|\phi\|_{H^{1}(\Omega_{p})}^{2}=\int_{\Omega_{p}}a(\varepsilon
y)\big{(}|\nabla\phi|^{2}+\varepsilon^{2}\phi^{2}\big{)}$. Equation (3.7) is
equivalent to find $\phi\in H_{\xi}$ such that

$\displaystyle\int_{\Omega_{p}}a(\varepsilon
y)\big{(}\nabla\phi\nabla\psi+\varepsilon^{2}\phi\psi\big{)}-\int_{\Omega_{p}}a(\varepsilon
y)W_{\xi^{\prime}}\phi\psi=\int_{\Omega_{p}}a(\varepsilon
y)h\psi\,\,\quad\,\,\forall\,\,\psi\in H_{\xi}.$

By Fredholm’s alternative this is equivalent to the uniqueness of solutions to
this problem, which is guaranteed by estimate (3.65). Finally, for $p\geq
p_{m}$ fixed, by density of $C^{0,\alpha}(\overline{\Omega}_{p})$ in
$(C(\overline{\Omega}_{p}),\,\|\cdot\|_{L^{\infty}(\Omega_{p})})$, we can
approximate $h\in C(\overline{\Omega}_{p})$ by smooth functions and, by (3.65)
and elliptic regularity theory, we find that for any $h\in
C(\overline{\Omega}_{p})$, problem (3.7) admits a unique solution which
belongs to $H^{2}(\Omega_{p})$ and satisfies the a priori estimate (3.13). The
proof is complete. $\square$

Remark 3.7. Given $h\in C(\overline{\Omega}_{p})$ with $\|h\|_{*}<\infty$, let
$\phi$ be the solution to (3.7) given by Proposition 3.3. Testing the first
equation of (3.7) against $a(\varepsilon y)\phi$, we get

$\displaystyle\|\phi\|_{H^{1}(\Omega_{p})}^{2}=\int_{\Omega_{p}}a(\varepsilon
y)W_{\xi^{\prime}}\phi^{2}+\int_{\Omega_{p}}a(\varepsilon y)h\phi.$

Furthermore, by (3.1) we obtain

$\displaystyle\|\phi\|_{H^{1}(\Omega_{p})}\leq
C\big{(}\|h\|_{*}+\|\phi\|_{L^{\infty}(\Omega_{p})}\big{)}.$

##  4\. The nonlinear problem

In order to solve problem (2.36) we first consider the auxiliary nonlinear
problem: for any points $\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$, we
find a function $\phi$ and scalars $c_{ij}$, $i=1,\ldots,m$, $j=1,J_{i}$ such
that

$\left\\{\begin{aligned}
&\mathcal{L}(\phi)=-\big{[}R_{\xi^{\prime}}+N(\phi)\big{]}+\frac{1}{a(\varepsilon
y)}\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{J_{i}}c_{ij}\chi_{i}\,Z_{ij}\,\,\
\,\textrm{in}\,\,\,\,\,\,\Omega_{p},\\\
&\frac{\partial\phi}{\partial\nu}=0\,\,\,\,\,\,\,\qquad\qquad\qquad\quad\qquad\qquad\qquad\qquad\qquad\textrm{on}\,\,\,\,\partial\Omega_{p},\\\\[2.84526pt]
&\int_{\Omega_{p}}\chi_{i}\,Z_{ij}\phi=0\qquad\qquad\qquad\quad\,\forall\,\,i=1,\ldots,m,\,\,\,j=1,J_{i},\end{aligned}\right.$
(4.1)

where $W_{\xi^{\prime}}=pV_{\xi^{\prime}}^{p-1}$ satisfies (3.1)-(3.2), and
$R_{\xi^{\prime}}$, $N(\phi)$ are given by (2.35).

Proposition 4.1. Let $m$ be a positive integer. Then there exist constants
$C>0$ and $p_{m}>1$ such that for any $p>p_{m}$ and any points
$\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$, problem (4.1) admits a
unique solution $\phi_{\xi^{\prime}}$ for some coefficients
$c_{ij}(\xi^{\prime})$, $i=1,\ldots,m$, $j=1,J_{i}$, such that

$\displaystyle\|\phi_{\xi^{\prime}}\|_{L^{\infty}(\Omega_{p})}\leq\frac{C}{p^{3}},\,\,\quad\,\quad\,\,\sum_{i=1}^{m}\sum_{j=1}^{J_{i}}\mu_{i}|c_{ij}(\xi^{\prime})|\leq\frac{C}{p^{4}}\quad\,\quad\textrm{and}\quad\,\quad\|\phi_{\xi^{\prime}}\|_{H^{1}(\Omega_{p})}\leq\frac{C}{p^{3}}.$
(4.2)

Furthermore, the map $\xi^{\prime}\mapsto\phi_{\xi^{\prime}}$ is a
$C^{1}$-function in $C(\overline{\Omega}_{p})$ and $H^{1}(\Omega_{p})$.

###### Proof.

Proposition $3.3$, Remarks $3.2$ and $3.7$ allow us to apply the Contraction
Mapping Principle and the Implicit Function Theorem to find a solution for
problem (4.1) satisfying (4.2) and the corresponding regularity. Since it is a
standard procedure, we omit the detailed proof here, see Lemma 4.1 in [21] for
a proof. ∎

Remark 4.2. The function $V_{\xi^{\prime}}+\phi_{\xi^{\prime}}$, where
$\phi_{\xi^{\prime}}$ is given by Proposition 4.1, is positive in
$\overline{\Omega}_{p}$. In fact, we observe first that
$p^{2}\phi_{\xi^{\prime}}\rightarrow 0$ uniformly in
$C(\overline{\Omega}_{p})$. Furthermore, from Remark 2.6 and the definition of
$V_{\xi^{\prime}}$ in (2.31) we have that, in each region
$|y-\xi^{\prime}_{i}|<1/(\varepsilon p^{2\kappa})$,
$V_{\xi^{\prime}}+\phi_{\xi^{\prime}}$ is positive. Outside these regions, by
(2.23) and (2.31) we may get the same result.

##  5\. Variational reduction

After problem (4.1) has been solved, we find a solution of problem (2.36) and
hence to the original problem (1.6) if we find $\xi^{\prime}$ such that the
coefficient $c_{ij}(\xi^{\prime})$ in (4.1) satisfies

$\displaystyle c_{ij}(\xi^{\prime})=0\,\quad\,\,\textrm{for
all}\,\,\,i=1,\ldots,m,\,\,\,j=1,J_{i}.$ (5.1)

Equation (1.6) is the Euler-Lagrange equation of the energy functional $J_{p}$
given by

$\displaystyle J_{p}(u)=\frac{1}{2}\int_{\Omega}a(x)(|\nabla
u|^{2}+u^{2})dx-\frac{1}{p+1}\int_{\Omega}a(x)u_{+}^{p+1}dx,\,\,\quad\,\,\,\,\,u\in
H^{1}(\Omega).$ (5.2)

We define

$\displaystyle F_{p}(\xi)=J_{p}(U_{\xi}+\widetilde{\phi}_{\xi})\,\ \,\,\ \
\,\forall\,\,\xi\in\mathcal{O}_{p},$ (5.3)

where $U_{\xi}$ is the function defined in (2.19) and

$\displaystyle\widetilde{\phi}_{\xi}(x)=\varepsilon^{-2/(p-1)}\phi_{\xi^{\prime}}(\varepsilon^{-1}x),\,\quad\,x\in\Omega,$
(5.4)

with $\phi_{\xi^{\prime}}$ the unique solution to problem (4.1) given by
Proposition 4.1. Critical points of $F_{p}$ correspond to solutions of (5.1)
for large $p$, as the following results states.

Proposition 5.1. The function $F_{p}:\mathcal{O}_{p}\mapsto\mathbb{R}$ is of
class $C^{1}$. Moreover, for all $p$ sufficiently large, if
$D_{\xi}F_{p}(\xi)=0$, then $\xi^{\prime}=\xi/\varepsilon$ satisfies (5.1).

###### Proof.

A direct consequence of the results obtained in Proposition 4.1 and the
definition of function $U_{\xi}$ is the fact that $F_{p}(\xi)$ is a
$C^{1}$-function of $\xi$ in $\mathcal{O}_{p}$ since the map
$\xi\mapsto\widetilde{\phi}_{\xi}$ is a $C^{1}$-map into $H^{1}(\Omega)$.
Recall that

$\displaystyle I_{p}(\upsilon)=\frac{1}{2}\int_{\Omega_{p}}a(\varepsilon
y)\left(|\nabla\upsilon|^{2}+\varepsilon^{2}\upsilon^{2}\right)dy-\frac{1}{p+1}\int_{\Omega_{p}}a(\varepsilon
y)\upsilon^{p+1}_{+}dy,\,\quad\,\upsilon\in H^{1}(\Omega_{p}),$

Then, making a change of variable, we get

$\displaystyle
F_{p}(\xi)=J_{p}\big{(}U_{\xi}+\widetilde{\phi}_{\xi}\big{)}=\varepsilon^{-4/(p-1)}I_{p}\big{(}V_{\xi^{\prime}}+\phi_{\xi^{\prime}}\big{)}.$
(5.5)

Since $\phi_{\xi^{\prime}}$ solves equation (4.1) and $D_{\xi}F_{p}(\xi)=0$,
we have that for any $k=1,\ldots,m$ and $t=1,J_{k}$,

$\displaystyle
0=I^{\prime}_{p}\big{(}V_{\xi^{\prime}}+\phi_{\xi^{\prime}}\big{)}\partial_{(\xi^{\prime}_{k})_{t}}\big{(}V_{\xi^{\prime}}+\phi_{\xi^{\prime}}\big{)}$
$\displaystyle\,\,\,\,\,=\int_{\Omega_{p}}a(x)\left[-\Delta_{a(\varepsilon
y)}\big{(}V_{\xi^{\prime}}+\phi_{\xi^{\prime}}\big{)}+\varepsilon^{2}\big{(}V_{\xi^{\prime}}+\phi_{\xi^{\prime}}\big{)}-\big{(}V_{\xi^{\prime}}+\phi_{\xi^{\prime}}\big{)}^{p}\right]\partial_{(\xi^{\prime}_{k})_{t}}\big{(}V_{\xi^{\prime}}+\phi_{\xi^{\prime}}\big{)}$
$\displaystyle\,\,\,\,\,=\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{J_{i}}c_{ij}(\xi^{\prime})\int_{\Omega_{p}}\chi_{i}Z_{ij}\partial_{(\xi^{\prime}_{k})_{t}}\big{(}V_{\xi^{\prime}}+\phi_{\xi^{\prime}}\big{)}$
$\displaystyle\,\,\,\,\,=\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{J_{i}}c_{ij}(\xi^{\prime})\int_{\Omega_{p}}\chi_{i}Z_{ij}\partial_{(\xi^{\prime}_{k})_{t}}V_{\xi^{\prime}}-\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{J_{i}}c_{ij}(\xi^{\prime})\int_{\Omega_{p}}\phi_{\xi^{\prime}}\partial_{(\xi^{\prime}_{k})_{t}}\big{(}\chi_{i}Z_{ij}\big{)},$
(5.6)

because $\int_{\Omega_{p}}\chi_{i}Z_{ij}\phi_{\xi^{\prime}}=0$. Notice first
that by (3.4), (3.8), (3.11) and (3.46), a direct computation shows

$\displaystyle\big{|}\partial_{(\xi^{\prime}_{k})_{t}}\big{(}\chi_{i}Z_{ij}\big{)}\big{|}=O\left(\frac{1}{\mu_{i}}\big{|}\partial_{(\xi^{\prime}_{k})_{t}}\log\mu_{i}\big{|}+\frac{1}{\mu_{i}^{2}}\delta_{ki}\right).$
(5.7)

On the other hand, since
$D_{\xi^{\prime}}V_{\xi^{\prime}}(y)=\varepsilon^{2/(p-1)}D_{\xi^{\prime}}U_{\xi}(\varepsilon
y)$, by (2.11), (2.12), (2.19), (2.23) and (2.24) we obtain

$\displaystyle\partial_{(\xi^{\prime}_{k})_{t}}V_{\xi^{\prime}}(y)=$
$\displaystyle\,\sum_{i=1}^{m}\frac{\varepsilon^{2/(p-1)}}{\,\,\gamma\mu_{i}^{2/(p-1)}\,\,}\left\\{\partial_{(\xi^{\prime}_{k})_{t}}\left[U_{\delta_{i},\xi_{i}}(\varepsilon
y)+\frac{1}{p}\omega_{1}\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)+\frac{1}{p^{2}}\omega_{2}\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)+\gamma\mu_{i}^{2/(p-1)}H_{i}(\varepsilon
y)\right]\right.$
$\displaystyle\left.-\frac{\,2\,\partial_{(\xi^{\prime}_{k})_{t}}\log\mu_{i}\,}{p-1}\big{[}\,p+O\left(1\right)\big{]}\right\\}.$

By (2.7) and (3.4),

$\displaystyle\partial_{(\xi^{\prime}_{k})_{t}}U_{\delta_{i},\xi_{i}}(\varepsilon
y)=$
$\displaystyle\frac{4}{\mu_{i}}\delta_{ki}Z_{t}\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)+2Z_{0}\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)\partial_{(\xi^{\prime}_{k})_{t}}\log\mu_{i}$
$\displaystyle=$
$\displaystyle\frac{4}{\mu_{i}}\delta_{ki}Z_{t}\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)+2\left(1-\frac{2\mu_{i}^{2}}{|y-\xi^{\prime}_{i}|^{2}+\mu_{i}^{2}}\right)\partial_{(\xi^{\prime}_{k})_{t}}\log\mu_{i},$

while for $j=1,2$, by (2.15) and (3.4),

$\displaystyle\partial_{(\xi^{\prime}_{k})_{t}}\omega_{j}\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)=$
$\displaystyle-\omega_{j}^{\prime}\left(\frac{|y-\xi^{\prime}_{i}|}{\mu_{i}}\right)\left\\{\frac{(y-\xi^{\prime}_{i})_{t}}{|y-\xi^{\prime}_{i}|}\frac{1}{\mu_{i}}\delta_{ki}+\frac{|y-\xi^{\prime}_{i}|}{\mu_{i}}\partial_{(\xi^{\prime}_{k})_{t}}\log\mu_{i}\right\\}$
$\displaystyle=$
$\displaystyle-\frac{1}{\mu_{i}}\delta_{ki}\left[C_{j}Z_{t}\left(\frac{y-\xi^{\prime}_{i}}{\mu_{i}}\right)+O\left(\frac{\mu_{i}^{2}}{|y-\xi^{\prime}_{i}|^{2}+\mu_{i}^{2}}\right)\right]-C_{j}\left[1+O\left(\frac{\mu_{i}^{2}}{|y-\xi^{\prime}_{i}|^{2}+\mu_{i}^{2}}\right)\right]\partial_{(\xi^{\prime}_{k})_{t}}\log\mu_{i}.$

Additionally, as in the proof of Lemma 2.1, we can prove that

$\displaystyle\partial_{(\xi^{\prime}_{k})_{t}}\left[\,\gamma\mu_{i}^{2/(p-1)}H_{i}(\varepsilon
y)\right]=$
$\displaystyle\,\delta_{ki}\left(1-\frac{C_{1}}{4p}-\frac{C_{2}}{4p^{2}}\right)c_{i}\partial_{(\xi^{\prime}_{k})_{t}}H(\varepsilon
y,\xi_{i})+\left(-2+\frac{C_{1}}{p}+\frac{C_{2}}{p^{2}}\right)\partial_{(\xi^{\prime}_{k})_{t}}\log\mu_{i}$
$\displaystyle+O\left(\frac{1}{p}\big{|}\partial_{(\xi^{\prime}_{k})_{t}}\log\mu_{i}\big{|}\right).$

Thus by (2.11) and (2.26),

$\displaystyle\partial_{(\xi^{\prime}_{k})_{t}}V_{\xi^{\prime}}(y)=\frac{1}{p^{p/(p-1)}\mu_{k}^{(p+1)/(p-1)}}\left\\{\left(1-\frac{C_{1}}{4p}-\frac{C_{2}}{4p^{2}}\right)\left[4Z_{t}\left(\frac{y-\xi^{\prime}_{k}}{\mu_{k}}\right)+c_{k}\mu_{k}\partial_{(\xi^{\prime}_{k})_{t}}H(\varepsilon
y,\xi_{k})\right]\right.$
$\displaystyle\left.+O\left(\frac{1}{p}\frac{\mu_{k}^{2}}{|y-\xi^{\prime}_{k}|^{2}+\mu_{k}^{2}}\right)\right\\}+O\left(\frac{1}{p}\sum_{i=1}^{m}\big{|}\partial_{(\xi^{\prime}_{k})_{t}}\log\mu_{i}\big{|}\right).\qquad\quad\qquad\qquad\qquad\,\,\,\,\,\quad$
(5.8)

Hence for each $k=1,\ldots,m$ and $t=1,J_{k}$, (5) can be written as

$\displaystyle\sum_{i,j}\frac{c_{ij}(\xi^{\prime})}{p^{p/(p-1)}\mu_{k}^{2/(p-1)}}\left\\{\delta_{ki}\left[\frac{c_{i}}{2\pi}\int_{\mathbb{R}^{2}}\chi(|z|)Z_{j}(z)Z_{t}(z)dz+O\left(\frac{1}{p}\right)\right]+(1-\delta_{ki})O\left(\frac{\mu_{i}}{|\xi^{\prime}_{i}-\xi^{\prime}_{k}|}\right)\right\\}$
$\displaystyle+\sum_{ij}|c_{ij}(\xi^{\prime})|\left\\{O\left(\frac{1}{p}\mu_{i}\sum_{s=1}^{m}\big{|}\partial_{(\xi^{\prime}_{k})_{t}}\log\mu_{s}\big{|}\right)+\|\phi_{\xi^{\prime}}\|_{L^{\infty}(\Omega_{p})}\,O\left(\mu_{i}\big{|}\partial_{(\xi^{\prime}_{k})_{t}}\log\mu_{i}\big{|}+\delta_{ki}\right)\right\\}=0,$

and then, by (2.9), (2.26), (2.27) and (4.2),

$\displaystyle\frac{c_{k}c_{kt}(\xi^{\prime})}{p^{p/(p-1)}\mu_{k}^{2/(p-1)}}\int_{0}^{R_{0}+1}\chi(r)\frac{r^{3}}{(\,1+r^{2})^{2}}dr+\sum_{i=1}^{m}\sum_{j=1}^{J_{i}}\big{|}c_{ij}(\xi^{\prime})\big{|}O\left(\frac{\delta_{ki}}{p^{2}\mu_{k}^{2/(p-1)}}+\varepsilon\mu_{i}p^{\kappa-1}\right)=0,$

which implies $c_{kt}(\xi^{\prime})=0$ for each $k=1,\ldots,m$ and
$t=1,J_{k}$. ∎

##  6\. Expansion of the energy

Proposition 6.1. Let $m$ be a positive integer. With the choice of $\mu_{i}$’s
given by (2.25), there exists $p_{m}>1$ such that for any $p>p_{m}$ and any
points $\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$, the following
expansion holds

$\displaystyle
F_{p}(\xi)=\frac{e}{2p}\sum_{i=1}^{m}c_{i}a(\xi_{i})\left[1-\frac{2\log
p}{p}+\frac{\mathcal{K}+2}{p}-\frac{1}{p}c_{i}H(\xi_{i},\xi_{i})-\frac{1}{p}\sum_{k\neq
i}^{m}c_{k}G(\xi_{i},\xi_{k})\right]+O\left(\frac{\log^{2}p}{p^{3}}\right),$
(6.1)

where

$\displaystyle\mathcal{K}=\frac{1}{8\pi}\int_{\mathbb{R}^{2}}\left[\frac{8}{(1+|z|^{2})^{2}}U_{1,0}(z)-\Delta\omega_{1}(z)\right]dz.$

###### Proof.

Multiplying the first equation in (4.1) by $a(\varepsilon
y)(V_{\xi^{\prime}}+\phi_{\xi^{\prime}})$ and integrating by parts, we obtain

$\displaystyle\int_{\Omega_{p}}a(\varepsilon
y)\left[|\nabla(V_{\xi^{\prime}}+\phi_{\xi^{\prime}})|^{2}+\varepsilon^{2}(V_{\xi^{\prime}}+\phi_{\xi^{\prime}})^{2}\right]dy$
$\displaystyle=\int_{\Omega_{p}}a(\varepsilon
y)(V_{\xi^{\prime}}+\phi_{\xi^{\prime}})^{p+1}dy+\sum_{i=1}^{m}\sum_{j=1}^{J_{i}}c_{ij}(\xi^{\prime})\int_{\Omega_{p}}\chi_{i}Z_{ij}V_{\xi^{\prime}}dy$
$\displaystyle=\int_{\Omega_{p}}a(\varepsilon
y)(V_{\xi^{\prime}}+\phi_{\xi^{\prime}})^{p+1}dy+O\left(\frac{1}{p^{4}}\right)$

because $V_{\xi^{\prime}}=O\left(1\right)$ and
$|c_{ij}(\xi^{\prime})|=O\left(p^{-4}\mu_{i}^{-1}\right)$. Hence, by (2.31),
(4.2) and (5.5) we can write

$\displaystyle
F_{p}(\xi)=\left(\frac{1}{2}-\frac{1}{p+1}\right)\varepsilon^{-4/(p-1)}\int_{\Omega_{p}}a(\varepsilon
y)\left[|\nabla(V_{\xi^{\prime}}+\phi_{\xi^{\prime}})|^{2}+\varepsilon^{2}(V_{\xi^{\prime}}+\phi_{\xi^{\prime}})^{2}\right]dy+O\left(\frac{1}{p^{5}}\right)\,\,\,\,\,\quad$
$\displaystyle=\left(\frac{1}{2}-\frac{1}{p+1}\right)\varepsilon^{-4/(p-1)}\left\\{\int_{\Omega_{p}}a(\varepsilon
y)\big{(}|\nabla
V_{\xi^{\prime}}|^{2}+\varepsilon^{2}V_{\xi^{\prime}}^{2}\big{)}dy\right.\qquad\qquad\qquad\quad\,\,\qquad\qquad\,\,\,\quad\,$
$\displaystyle\left.+2\int_{\Omega_{p}}a(\varepsilon y)\big{(}\nabla
V_{\xi^{\prime}}\nabla\phi_{\xi^{\prime}}+\varepsilon^{2}V_{\xi^{\prime}}\phi_{\xi^{\prime}}\big{)}dy+\int_{\Omega_{p}}a(\varepsilon
y)\big{(}|\nabla\phi_{\xi^{\prime}}|^{2}+\varepsilon^{2}\phi_{\xi^{\prime}}^{2}\big{)}dy\right\\}+O\left(\frac{1}{p^{5}}\right)\quad$
$\displaystyle=\left(\frac{1}{2}-\frac{1}{p+1}\right)\int_{\Omega}a(x)\big{(}|\nabla
U_{\xi}|^{2}+U_{\xi}^{2}\big{)}dx+O\left(\frac{1}{p^{3}}\left|\int_{\Omega}a(x)\big{(}|\nabla
U_{\xi}|^{2}+U_{\xi}^{2}\big{)}dx\right|^{1/2}+\frac{1}{p^{5}}\right)$

uniformly for any points $\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$.
Now, in view of (2.11), (2.12), (2.19), (2.20), (2.21), (2.23) and (2.24) we
have

$\displaystyle\int_{\Omega}a(x)$ $\displaystyle\big{(}|\nabla
U_{\xi}|^{2}+U_{\xi}^{2}\big{)}dx=\int_{\Omega}a(x)(-\Delta_{a}U_{\xi}+U_{\xi})U_{\xi}dx$
$\displaystyle=\sum_{i=1}^{m}\frac{1}{\gamma\mu_{i}^{2/(p-1)}}\int_{\Omega\bigcap
B_{\frac{1}{p^{2\kappa}}}(\xi_{i})}a(x)\left[e^{U_{\delta_{i},\xi_{i}}}-\frac{1}{p\delta_{i}^{2}}\Delta\omega_{1}\left(\frac{x-\xi_{i}}{\delta_{i}}\right)-\frac{1}{p^{2}\delta_{i}^{2}}\Delta\omega_{2}\left(\frac{x-\xi_{i}}{\delta_{i}}\right)\right]U_{\xi}dx+o\left(\delta_{i}\right)$
$\displaystyle=\sum_{i=1}^{m}\frac{1}{\gamma^{2}\mu_{i}^{4/(p-1)}}\int_{\big{(}\frac{1}{\delta_{i}}(\Omega-\\{\xi_{i}\\})\big{)}\bigcap
B_{\frac{1}{\delta_{i}p^{2\kappa}}}(0)}a(\xi_{i}+\delta_{i}z)\left[\frac{8}{(1+|z|^{2})^{2}}-\frac{1}{p}\Delta\omega_{1}(z)-\frac{1}{p^{2}}\Delta\omega_{2}(z)\right]$
$\displaystyle\,\,\,\qquad\,\times\left[p+U_{1,0}(z)+\frac{1}{p}\omega_{1}(z)+\frac{1}{p^{2}}\omega_{2}(z)+O\left(\delta_{i}^{\alpha}|z|^{\alpha}+\sum_{k=1}^{m}\delta_{k}^{\alpha/2}\right)\right]dz+o\left(\delta_{i}\right)$
$\displaystyle=\sum_{i=1}^{m}\frac{c_{i}a(\xi_{i})}{\gamma^{2}\mu_{i}^{4/(p-1)}}\left[\,p+\mathcal{K}+O\left(\frac{1}{p}\right)\right]+o\left(\delta_{i}\right).$

Recalling that $\gamma=p^{p/(p-1)}e^{-p/(2p-2)}$, we get

$\displaystyle\frac{1}{\gamma^{2}}=\frac{e}{p^{2}}\left[1-\frac{2\log
p}{p}+\frac{1}{p}+O\left(\frac{\log^{2}p}{p^{2}}\right)\right],\,\qquad\qquad\,\frac{1}{\,\mu_{i}^{4/(p-1)}}=1-\frac{4\log\mu_{i}}{p}+O\left(\frac{\log^{2}\mu_{i}}{p^{2}}\right),$

and then, by (2.26),

$\displaystyle\frac{1}{\,\gamma^{2}\mu_{i}^{4/(p-1)}\,}=\frac{e}{p^{2}}\left[1-\frac{2\log
p}{p}-\frac{4\log\mu_{i}}{p}+\frac{1}{p}+O\left(\frac{\log^{2}p}{p^{2}}\right)\right].$

Therefore,

$\displaystyle
F_{p}(\xi)=\frac{e}{2p}\sum_{i=1}^{m}c_{i}a(\xi_{i})\left\\{1-\frac{2\log
p}{p}-\frac{4\log\mu_{i}}{p}+\frac{\mathcal{K}-1}{p}\right\\}+O\left(\frac{\log^{2}p}{p^{3}}\right),$

which, together with the expansion of $\mu_{i}$ in (2.28), easily implies that
(6.1) holds. ∎

##  7\. Proofs of theorems

Proof of Theorem 1.1. According to Proposition 5.1,
$U_{\xi}+\widetilde{\phi}_{\xi}$ is a solution of problem (1.6) if we adjust
$\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$ so that it is a critical
point of $F_{p}$ defined in (5.3). For this aim, we only choose points
$\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$ in the form of the
parametrization

$\displaystyle\xi_{i}\equiv\xi_{i}(\mathbf{s},\mathbf{t})=s_{i}-\frac{t_{i}}{p}\nu(s_{i}),\quad
i=1,\ldots,l,\qquad\quad\textrm{but}\quad\qquad\xi_{i}\equiv\xi_{i}(\mathbf{s},\mathbf{t})=s_{i},\quad
i=l+1,\ldots,m,$ (7.1)

where $\mathbf{s}=(s_{1},\ldots,s_{m})\in(\partial\Omega)^{m}$ and
$\mathbf{t}=(t_{1},\ldots,t_{l})\in\mathbb{R}^{l}_{+}$ lie in the
configuration space

$\displaystyle\Lambda_{d}=\left\\{(\mathbf{s},\,\mathbf{t})\in(\partial\Omega)^{m}\times\mathbb{R}^{l}_{+}\,\big{|}\,|s_{i}-s_{k}|>2d,\quad\,d<t_{\tilde{i}}<1/d,\quad\forall\,\,\,i,k=1,\ldots,m,\,\,\,\tilde{i}=1,\ldots,l,\,\,\,i\neq
k\right\\},$

for any $d>0$ sufficiently small, fixed and independent of $p$. Thus we can
easily prove that if $(\mathbf{s},\mathbf{t})$ is a critical point of the
reduced energy
$\widehat{F}_{p}\big{(}\mathbf{s},\mathbf{t}\big{)}:=F_{p}\big{(}\xi(\mathbf{s},\mathbf{t})\big{)}$
in $\Lambda_{d}$, then the function
$U_{\xi(\mathbf{s},\mathbf{t})}+\widetilde{\phi}_{\xi(\mathbf{s},\mathbf{t})}$
is a solution of problem (1.6) with the qualitative properties predicted by
Theorem 1.1. Therefore, we need first to compute the expansion of the reduced
energy $\widehat{F}_{p}$ with the aid of Lemmas 2.1-2.2, Corollary 2.3 and
Proposition 6.1.

Using the smooth property of $a(x)$ over $\overline{\Omega}$, we can perform a
Taylor expansion around each boundary point $s_{i}$ along the inner normal
vector $-\nu(s_{i})$ to derive that

$\displaystyle
a(\xi_{i})=a(s_{i})-\frac{t_{i}}{p}\partial_{\nu}a(s_{i})+O\left(\frac{t_{i}^{2}}{p^{2}}\right),\quad\forall\,\,i=1,\ldots,l.$
(7.2)

From the expansions of Robin’s function in (2.5)-(2.6) and the regularity of
the vector function $T(x)$, we obtain

$\displaystyle
H(\xi_{i},\xi_{i})=-\frac{1}{2\pi}\log\left(\frac{2t_{i}}{p}\right)+\tilde{\mathrm{z}}(s_{i},s_{i})+O\left(\frac{t_{i}^{\alpha}}{p^{\alpha}}\right),\quad\forall\,\,i=1,\ldots,l.$
(7.3)

On the other hand, from the expansions of the regular part of Green’s function
in (2.3)-(2.4), we find that if $i,k=1,\ldots,l$ with $i\neq k$,

$\displaystyle
G(\xi_{i},\xi_{k})=-\frac{1}{2\pi}\log\left|s_{i}-s_{k}-\frac{t_{i}}{p}\nu(s_{i})+\frac{t_{k}}{p}\nu(s_{k})\right|-\frac{1}{2\pi}\log\left|s_{i}-s_{k}-\frac{t_{i}}{p}\nu(s_{i})-\frac{t_{k}}{p}\nu(s_{k})\right|\quad\qquad\quad$
$\displaystyle+\frac{t_{k}}{\pi p}\left\langle\nabla\log a(s_{k}),\,\nabla
T(s_{i}-s_{k})\cdot\nu(s_{k})\right\rangle-\frac{t_{k}}{\pi
p}\left\langle\big{(}\nabla\times(\nabla\log
a)\big{)}(s_{k})\cdot\nu(s_{k}),\,T(s_{i}-s_{k})\right\rangle$
$\displaystyle+\tilde{\mathrm{z}}(s_{i},s_{k})-\frac{1}{p}\left\langle\nabla_{(s_{i},s_{k})}\tilde{\mathrm{z}}(s_{i},s_{k}),\,\big{(}t_{i}\nu(s_{i}),\,t_{k}\nu(s_{k})\big{)}\right\rangle+O\left(\frac{t_{i}^{2}+t_{k}^{2}}{p^{2}}\right),\qquad\qquad\qquad\quad\,\,\,\,$
(7.4)

while if $i=l+1,\ldots,m$ and $k=1,\ldots,l$,

$\displaystyle
G(\xi_{i},\xi_{k})=-\frac{1}{2\pi}\log\left|s_{i}-s_{k}+\frac{t_{k}}{p}\nu(s_{k})\right|-\frac{1}{2\pi}\log\left|s_{i}-s_{k}-\frac{t_{k}}{p}\nu(s_{k})\right|\qquad\qquad\qquad\qquad\quad\qquad\quad\,\,\,$
$\displaystyle+\frac{t_{k}}{\pi p}\left\langle\nabla\log a(s_{k}),\,\nabla
T(s_{i}-s_{k})\cdot\nu(s_{k})\right\rangle-\frac{t_{k}}{\pi
p}\left\langle\big{(}\nabla\times(\nabla\log
a)\big{)}(s_{k})\cdot\nu(s_{k}),\,T(s_{i}-s_{k})\right\rangle$
$\displaystyle+\tilde{\mathrm{z}}(s_{i},s_{k})-\frac{t_{k}}{p}\left\langle\nabla_{s_{k}}\tilde{\mathrm{z}}(s_{i},s_{k}),\,\nu(s_{k})\right\rangle+O\left(\frac{t_{k}^{2}}{p^{2}}\right).\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\,$
(7.5)

Substituting (2.21), (7.2)-(7) into (6.1) and using the fact that
$a(\xi_{i})G(\xi_{i},\xi_{k})=a(\xi_{k})G(\xi_{k},\xi_{i})$ for all
$i,k=1,\ldots,m$ with $i\neq k$, we conclude that $\widehat{F}_{p}$ becomes

$\displaystyle\widehat{F}_{p}\big{(}\mathbf{s},\mathbf{t}\big{)}=\frac{e}{2p}\left\\{8\pi\sum_{i=1}^{l}\left[a(s_{i})+\frac{1}{p}\large\big{(}4a(s_{i})\log
t_{i}-t_{i}\partial_{\nu}a(s_{i})\large\big{)}\right]+4\pi\sum_{k=l+1}^{m}a(s_{k})+\Gamma_{p}(\mathbf{s})\right\\}+O\left(\frac{\,|\mathbf{t}|^{\alpha}}{p^{2+\alpha}}+\frac{\,\log^{2}p}{p^{3}}\right)$
(7.6)

$C^{1}$-uniformly in $\Lambda_{d}$, where $\Gamma_{p}(\mathbf{s})$ is a smooth
function of $\mathbf{s}$ satisfying that $\Gamma_{p}(\mathbf{s})$ and
$|\nabla\Gamma_{p}(\mathbf{s})|$ uniformly converge to zero as
$p\rightarrow+\infty$.

We seek a critical point of $\widehat{F}_{p}$ in $\Lambda_{d}$. Let
$\partial_{T(s_{i})}$ be the tangential derivative which is defined on
$s_{i}\in\partial\Omega$. Set

$\displaystyle A(s_{i},t_{i})=a(s_{i})+\frac{1}{p}\large\big{(}4a(s_{i})\log
t_{i}-t_{i}\partial_{\nu}a(s_{i})\large\big{)},\,\,\quad\,i=1,\ldots,l.$

Then

$\displaystyle\partial_{T(s_{i})}A(s_{i},t_{i})=\left(1+\frac{4\log
t_{i}}{p}\right)\partial_{T}a(s_{i})-\frac{t_{i}}{p}\partial_{T}\partial_{\nu}a(s_{i})\large,\,\qquad\,\partial_{t_{i}}A(s_{i},t_{i})=\frac{1}{p}\left(\frac{4}{t_{i}}a(s_{i})-\partial_{\nu}a(s_{i})\right).$

In view of $\partial_{\nu}a(\xi^{*}_{i})>0$ with $i=1,\ldots,l$, we have that
for any sufficiently small $d>0$ and for any $s_{i}\in
B_{d}(\xi^{*}_{i})\cap\partial\Omega$, there exists a unique positive
$t_{i}=t_{i}(s_{i})=4a(s_{i})/\partial_{\nu}a(s_{i})$ such that
$\partial_{t_{i}}A(s_{i},t_{i})=0$ and
$\partial^{2}_{t_{i}t_{i}}A(s_{i},t_{i})<0$. Set
$t_{i}^{*}=t_{i}(\xi_{i}^{*})$, $i=1,\ldots,l$. Since
$\xi^{*}_{1},\ldots,\xi^{*}_{m}$ are $m$ different strict local maximum or
strict local minimum points of $a(x)$ on $\partial\Omega$, we find that for
any $i=1,\ldots,l$, any sufficiently small $d>0$ and any sufficiently large
$p>1$, the Brouwer degrees
$\deg\big{(}\big{(}\partial_{T(s_{i})}A,\,\partial_{t_{i}}A\big{)},\,\big{(}B_{d}(\xi^{*}_{i})\cap\partial\Omega\big{)}\times\big{(}t_{i}^{*}-d,\,t_{i}^{*}+d\big{)},\,0\big{)}$
and
$\deg\big{(}\big{(}\partial_{T(s_{l+1})}a,\ldots,\partial_{T(s_{m})}a\big{)},\,\prod_{k=l+1}^{m}\big{(}B_{d}(\xi^{*}_{k})\cap\partial\Omega\big{)},\,0\big{)}$
are well defined. By the definition and homotopy invariance of topological
degree we deduce

$\displaystyle\deg\big{(}\big{(}\partial_{T(s_{i})}A,$
$\displaystyle\,\partial_{t_{i}}A\big{)},\,\big{(}B_{d}(\xi^{*}_{i})\cap\partial\Omega\big{)}\times\big{(}t_{i}^{*}-d,\,t_{i}^{*}+d\big{)},\,0\big{)}$
$\displaystyle=$
$\displaystyle\deg\big{(}\big{(}\partial_{T}a(s_{i}),\,\partial_{t_{i}}A\big{)},\,\big{(}B_{d}(\xi^{*}_{i})\cap\partial\Omega\big{)}\times\big{(}t_{i}^{*}-d,\,t_{i}^{*}+d\big{)},\,0\big{)}$
$\displaystyle=$ $\displaystyle\text{sign}\,\det\left[\left(\begin{aligned}
\partial^{2}_{TT}a(\xi^{*}_{i})\quad&\quad\frac{1}{p}\left(\frac{4}{t^{*}_{i}}\partial_{T}a(\xi^{*}_{i})-\partial_{T}\partial_{\nu}a(\xi^{*}_{i})\right)\\\
0\quad\quad\quad&\quad\quad\quad-\frac{1}{p}\frac{4}{(t^{*}_{i})^{2}}a(\xi^{*}_{i})\\\
\end{aligned}\right)\right]=\pm 1\neq 0,$

and

$\displaystyle\deg\left(\big{(}\partial_{T(s_{l+1})}a,\ldots,\partial_{T(s_{m})}a\big{)},\,\prod_{k=l+1}^{m}\big{(}B_{d}(\xi^{*}_{k})\cap\partial\Omega\big{)},\,0\right)=\text{sign}\left(\prod_{k=l+1}^{m}\partial^{2}_{TT}a(\xi^{*}_{k})\right)=\pm
1\neq 0.$

Furthermore, using the properties of Brouwer degree, by (7.6) we conclude

$\displaystyle\,\deg\left(\nabla_{\left(T(\mathbf{s}),\mathbf{t}\right)}\widehat{F}_{p}\big{(}\mathbf{s},\mathbf{t}\big{)},\,\,\prod_{i=1}^{l}\big{(}B_{d}(\xi^{*}_{i})\cap\partial\Omega\big{)}\times\prod_{k=l+1}^{m}\big{(}B_{d}(\xi^{*}_{k})\cap\partial\Omega\big{)}\times\prod_{i=1}^{l}\big{(}t_{i}^{*}-d,\,t_{i}^{*}+d\big{)},\,\,0\right)$
$\displaystyle=$
$\displaystyle\deg\left(\big{(}\partial_{T(s_{1})}A,\,\ldots,\,\partial_{T(s_{l})}A,\,\,\partial_{T(s_{l+1})}a,\,\ldots,\,\partial_{T(s_{m})}a,\,\,\partial_{t_{1}}A,\,\ldots,\,\partial_{t_{l}}A\big{)},\right.$
$\displaystyle\left.\qquad\quad\,\,\prod_{i=1}^{l}\big{(}B_{d}(\xi^{*}_{i})\cap\partial\Omega\big{)}\times\prod_{k=l+1}^{m}\big{(}B_{d}(\xi^{*}_{k})\cap\partial\Omega\big{)}\times\prod_{i=1}^{l}\big{(}t_{i}^{*}-d,\,t_{i}^{*}+d\big{)},\,\,0\right)$
$\displaystyle=$