1. Introduction
This paper concerns the existence of multibump solutions to the following Hartree equation
\begin{align}
-\Delta u+\lambda V(x) u+\mu u=\left(\int_{\mathbb{R}^N}\frac{|u|^p}{|x-y|^{N-\alpha}}\right)|u|^{p-2}u,\ \text{in}\ \mathbb{R}^N
\end{align}under the mass constraint
\begin{align}
\int_{\mathbb{R}^N}|u|^2=\omega,
\end{align} where 
$N\geq 3$, 
$\omega,\lambda \gt 0$, 
$p\in \left(\frac{N+\alpha}{N}, \frac{N+\alpha}{N-2}\right)\setminus\left\{\frac{N+\alpha+2}{N}\right\}$, 
$0\leq V(x)\in L^{\infty}_{loc}(\mathbb{R}^N)$ and the zero set 
$V^{-1}(0)=\{x\in \mathbb{R}^N|\ V(x)=0\}$ of V have a nonempty interior 
$\Omega:=int\ V^{-1}(0)$, 
$\mu\in \mathbb{R}$ appears as a Lagrange multiplier and is part of the unknowns. It is easy to know that 
$\frac{N+\alpha+2}{N}$ is the L 2-critical exponent of (1.1) and 
$\frac{N+\alpha}{N-2}$ is the usual critical Sobolev exponent of (1.1). We call u a normalized solution to (1.1), since (1.2) imposes a normalization on its L 2-mass.
When 
$V(x)=0$, for example, the following problem
\begin{align}
-\Delta u+ \mu u=\left(\int_{\mathbb{R}^N}\frac{|u|^p}{|x-y|^{N-\alpha}}\right)|u|^{p-2}u,\ \text{in}\ \mathbb{R}^N
\end{align} where 
$N\geq 1$, 
$\mu\in \mathbb{R}$, 
$0 \lt \alpha \lt N$, 
$\frac{N+\alpha}{N}\leq p\leq \frac{N+\alpha}{N-2}$, arises in various fields of mathematical physics. The nonlinear term 
$\left(\int_{\mathbb{R}^N}\frac{|u|^p}{|x-y|^{N-\alpha}}\right)|u|^{p-2}u$ describes the interactions between the particles. Physically, 
$|x|^{-(N-\alpha)}$ is the so-called nonlocal response or Kernel, which represents the nonlocal character of the nonlinearity of the medium. Its width determines the degree of non-locality, see [Reference Lieb28]. 
$|x|^{-(N-\alpha)}$ is a short range potential, if 
$N-\alpha \gt 1$ (long range potential, if 
$N-\alpha\leq1$, see [Reference Ginibre and Velo18]), and sometimes even considers Dirac. This leads to the so-called Grose–Pitaevsky or cubic Schrödinger equation, which describes the standard Kerr-type local medium. In the physical case 
$N=3,p=2$, µ = 1 and α = 2, (1.3) appeared at least as early as in 1954, in a work by S.I. Pekar describing the quantum theory of a polariton at rest [Reference Pekar38]. In [Reference Lieb28], Choquard applied it as an approximation to the Hartree–Fock theory of one-component plasma. It also arises in multiple particle systems [Reference Gross19] and quantum mechanics [Reference Lieb28, Reference Penrose40]. In this context, (1.3) is usually called the nonlinear Schrödinger–Newton equation.
The existence and qualitative properties of solutions to the Hartree equation (1.1) and (1.3) have been studied for a few decades by variational methods. If 
$N=3,p=2$, µ = 1 and α = 2, the existence of ground states to equation (1.3) was obtained in [Reference Lieb28, Reference Lions30] by using variational methods. Concerning the qualitative properties of the ground state, uniqueness is proved in [Reference Lieb28]. In [Reference Moroz and Van Schaftingen36], Moroz and Schaftingen established some results about qualitative properties and decay estimates of ground state solutions to (1.3) with p > 1. Ye [Reference Ye51] studied the qualitative properties of the functional related to (1.1), including the existence and nonexistence of minimizers. Cingolani and Tanaka [Reference Cingolani and Tanaka10] developed the deformation arguments and proved the existence of a ground state solution for the nonlinear Hartree equation with prescribed mass, when the nonlinear term satisfies Berestycki–Lions type conditions. In [Reference Cingolani and Tanaka10, Reference Ye51], the existence of L 2-normalized solutions to equation (1.3) was investigated when 
$p\in\left(\frac{N+\alpha}{N},\frac{N+\alpha+2}{N}\right)$, i.e., in the L 2-subcritical case. Using constrained minimization method, Luo [Reference Luo33] proved the existence of a normalized solution for problem (1.3) with least energy when α = 2 and 
$2 \lt p \lt \{4,N\}$. When the nonlinear term satisfies different assumptions, we can see [Reference Bartsch, Liu and Liu2, Reference Cingolani, Gallo and Tanaka8, Reference Deng, Lu and Shuai12, Reference Li and Ye27, Reference Xia and Zhang46, Reference Yang and Zhu47] and the references therein.
When 
$V(x)\neq0$, there have been some works concerning the study of problem (1.1) under various assumptions on the potential. We refer the readers to [Reference Alves, Gao, Squassina and Yang1, Reference Cao, Wang and Zou5, Reference Chen and Guo7, Reference Cingolani and Tanaka9, Reference Ding, Gao and Yang13–Reference Gao and Wang17, Reference Guo, Hu, Peng and Shuai20–Reference Hironobu23, Reference Li and Zhang26, Reference Peng and Rizzi39, Reference Su and Shi42, Reference Wang44, Reference Yao, Chen, DRădulescu and Sun49, Reference Yao, Sun and Wu50, Reference Yu and Chen52] and the references therein. Among them, the authors are mainly interested in the problem of how the properties of the existence, multiplicity, and concentration of solutions. When 
$V(x)=k|x|^2$, Wang [Reference Wang44] proved the existence of ground state for problem (1.1) by the minimization method on the Nehari manifold. Luo [Reference Luo33] considered the existence and stability/instability of standing waves with a prescribed L 2-norm for problem (1.1) with 
$\alpha=2,2 \lt p \lt \min\{4,N\}$. Cao [Reference Cao, Wang and Zou5] et al. considered the problem with 
$a|x|^2$ replaced by a more general harmonic potential. They proved the existence and stability of standing waves with a prescribed L 2-norm for problem (1.1). In the semi-linear case, Wei and Winter [Reference Wei and Winter45] constructed families of solutions concentrating to the nondegenerate critical points of potential V by using Lyapunov–Schmidt reduction when N = 3 and 
$p=\alpha=2$. Using the nonlocal penalization method, Moroz et al. [Reference Moroz and Van Schaftingen37] proved existence of a single-peak solutions of (1.1) concentrating at a local minima of V(x) for 
$p\in\left[2,\frac{N+\alpha}{N-2}\right)$. When 
$p\in\left(\frac{N+\alpha}{N},\frac{N+\alpha}{N-2}\right)$, Cingolani and Tanaka [Reference Cingolani and Tanaka9] obtained the multiplicity of concentrating solutions by the cup-length of a critical set 
$Crit_{V_0}=\{x\in\Omega:\ V(x)=V_0\}$, where 
$\Omega\subset\mathbb{R}^N$ is a bounded set such that 
$V_0\equiv\inf\limits_{x\in\Omega}V(x) \lt \inf\limits_{x\in\partial\Omega}V(x)$. Subsequently, they [Reference Cingolani and Tanaka11] studied the existence of a family of solutions concentrating to a local maximum or a saddle point of 
$V(x)\in C^{N}(\mathbb{R}^N, \mathbb{R})$. The results for the critical problem, we can see [Reference Alves, Gao, Squassina and Yang1, Reference Gao and Yang16, Reference Guo, Hu, Peng and Shuai20–Reference He, Radulescu and Zou22, Reference Liu and Chen32, Reference Su and Liu41, Reference Yang and Ding48, Reference Zhang and Zhang54]. For more information on global well-posedness, blow-up solutions, and scattering for the Hartree equation (1.1), we refer the reader to [Reference Feng14, Reference Gao and Wang17, Reference Hironobu23, Reference Li and Zhang26, Reference Miao, Xu and Zhao34] and the references therein.
Motivated by the above papers, the aim of this paper is to consider the existence of normalized multibump solutions for equation (1.1) in the case of a steep potential well. To illustrate the results of the problem (1.1) in this direction, we assume that:
1. (V 1)
$\Omega=\cup_{1\leq j\leq\ell_0}\Omega_j$
$(\ell_0\in \mathbb{N}\setminus\{0\})$, where
$\Omega_{1},\cdots,\Omega_{\ell_0}$ are nonempty bounded domains with C 1 smooth boundary such that
$\overline{\Omega_j}\cap\overline{\Omega\setminus\Omega_j}= \emptyset$ for each
$j=1,\cdots,\ell_0$.2. (V 2)
$V\in L^{\infty}_{loc}(\mathbb{R}^N)$ and
$V(x)\geq0$ for a.e.
$x\in\mathbb{R}^N$;
$\Omega:=int V^{-1}(0)$ satisfies
$(V_1)$, and V > 0 a.e., in
$\mathbb{R}^N\setminus\Omega$.3. (V 3) There exists
$M_0 \gt 0$ such that
\begin{align*}
meas\left(\{x\in\mathbb{R}^N|\ V(x)\leq M_0\}\right) \lt \infty.
\end{align*}
Let 
$\ell\in[1, \ell_0]$ be an integer. For every fixed subset
\begin{equation*}
L:=\{j_1,j_2, \cdots, j_{\ell}\}\subset\{1, \cdots, \ell_0\},
\end{equation*} we set 
$\Omega_{L}:=\cup_{i\in L}\Omega_i.$ By scaling, we have
\begin{equation*}
v(\cdot)=\omega^{\frac{-(2+\alpha)}{2(-Np+N+\alpha+2)}}u(\omega^{\frac{-(p-1)}{-Np+N+\alpha+2}}\cdot)
\end{equation*} for every 
$u\in H^1(\mathbb{R}^N)$ with 
$\int_{\mathbb{R}^N}|u|^2=\omega$. Then (1.1) is equivalent to
\begin{align*}
\begin{cases}
-\Delta v+\lambda \omega^{\frac{-2(p-1)}{-Np+N+\alpha+2}}V\left(\omega^{\frac{-(p-1)}{-Np+N+\alpha+2}} x\right)v+\bar{\mu} v=\left(\int_{\mathbb{R}^N}\frac{|v|^p}{|x-y|^{N-\alpha}}dy\right)|v|^{p-2}v,\\ \qquad \text{in}\ \mathbb{R}^N,\\
\int_{\mathbb{R}^N}|v(x)|^{2}=1.
\end{cases}
\end{align*} where 
$\bar{\mu}=\omega^{\frac{-2(p-1)}{-Np+N+\alpha+2}}\mu$. For convenience, let
\begin{align*}
\varepsilon=\omega^{\frac{-(p-1)}{-Np+N+\alpha+2}}.
\end{align*}Then (1.1) is equivalent to
\begin{align}
\begin{cases}
-\Delta v+\lambda \varepsilon^2V\left(\varepsilon x\right)v+\bar{\mu} v=\left(\int_{\mathbb{R}^N}\frac{|v|^p}{|x-y|^{N-\alpha}}dy\right)|v|^{p-2}v, \quad \text{in}\ \mathbb{R}^N,\\
\int_{\mathbb{R}^N}|v(x)|^{2}=1.
\end{cases}
\end{align}Clearly, there holds
\begin{align*}
\varepsilon\rightarrow0\ \text{as}\ \omega\rightarrow
\begin{cases}
+\infty,\quad &\text{if}\ p\in \left(\frac{N+\alpha}{N}, \frac{N+\alpha+2}{N}\right),\\
0^+,\quad &\text{if}\ p\in \left(\frac{N+\alpha+2}{N}, \frac{N+\alpha}{N-2}\right).
\end{cases}
\end{align*} Set 
$V_{\varepsilon}(\cdot)=V(\varepsilon\cdot)$ and consider the functional corresponding to (1.4) by
\begin{equation*}
J_{\lambda,\varepsilon}(u)=\frac{1}{2}\int_{\mathbb{R}^N}\left(|\nabla u|^2+\lambda\varepsilon^2V_{\varepsilon}u^2\right)-\frac{1}{2p}\int_{\mathbb{R}^N}
\int_{\mathbb{R}^N}\frac{|u(x)|^p|u(y)|^p}{|x-y|^{N-\alpha}}dxdy,
\end{equation*}which is defined on
\begin{equation*}
E_{\varepsilon}:=\left\{u\in H^1(\mathbb{R}^N)|\ \int_{\mathbb{R}^N}V_{\varepsilon}u^2 \lt +\infty\right\}
\end{equation*}endowed with the norm
\begin{equation*}
||u||_{\varepsilon}:=\left(||u||^2+\int_{\mathbb{R}^N}\varepsilon^2 V_{\varepsilon}u^2\right)^{\frac{1}{2}},
\end{equation*} where 
$||\cdot||$ is the standard 
$H^1(\mathbb{R}^N)$ norm. Normalized solutions to (1.4) can be obtained by searching critical points of 
$J_{\lambda, \varepsilon}$ on the constraint
\begin{equation*}
M_{\varepsilon}=\left\{u\in E_{\varepsilon}|\int_{\mathbb{R}^N} |u|^2=1\right\}.
\end{equation*}From the above analysis, we have the following results.
Theorem 1.1 Suppose the conditions 
$(V_1)-(V_3)$ hold. Then for every ɛ sufficiently small, there exists 
$\lambda_{\varepsilon} \gt 0$ such that if 
$\lambda \gt \lambda_{\varepsilon}$, problem (1.1) has a solution 
$(\mu_{\lambda, \omega}, v_{\lambda, \omega})$ with 
$\mu_{\lambda, \omega} \gt 0$ and 
$v_{\lambda, \omega} \gt 0$ in 
$\mathbb{R}^N$.
According to the conditions 
$(V_1)-(V_3)$, we deduce that as 
$\lambda\rightarrow+\infty$, the limit problem related to (1.1) is
\begin{align}
\begin{cases}
-\Delta v+\mu v=\left(\int_{\mathbb{R}^N}\frac{|v|^p}{|x-y|^{N-\alpha}}dy\right)|v|^{p-2}v, \quad \text{in}\ \Omega,\\
v \gt 0\ \text{in}\ \Omega_{L},\ v=0\ \text{in}\ \Omega\setminus\Omega_{L},\ \int_{\Omega_{L}}|v(x)|^{2}=\omega.
\end{cases}
\end{align} Then it follows from the asymptotic behaviour of the solution 
$(\mu_{\lambda, \omega}, v_{\lambda, \omega})$ given by Theorem 1.1, we have the following results about problem (1.5).
Theorem 1.2 Suppose 
$(V_1)$ holds. Then for sufficiently small ɛ > 0, problem (1.5) admits a normalized solution 
$v_{\varepsilon}\in H^1_0(\Omega)$ with an associated Lagrange multiplier 
$\mu=\mu_{\omega}$. Moreover, there exists 
$x^j_{\omega}\in \Omega_j$, for each 
$j\in L$ such that when 
$\varepsilon\rightarrow0$, there holds
\begin{align*}
&\omega^{\frac{-2(p-1)}{-Np+N+\alpha+2}}\mu_{\omega}\rightarrow\mu_{\ell^{-1}},\quad \text{in}\ \mathbb{R},\\
&\omega^{-\frac{2+\alpha}{2(-Np+N+\alpha+2)}}v_{\varepsilon}(\omega^{-\frac{p-1}{-Np+N+\alpha+2}}\cdot)-
\sum\limits_{j=1}^{\ell}
W_j\left(\cdot-\omega^{-\frac{p-1}{-Np+N+\alpha+2}}x^j_{\omega}\right)\rightarrow0,\\ &\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \text{in}\ H^1(\mathbb{R}^N),
\end{align*}where
\begin{equation*}
\mu_{\ell^{-1}}:=\frac{2(-Np+N+\alpha+2p)\ell I_{\ell^{-1}}}{Np-N-2-\alpha} \gt 0,
\end{equation*}and Wj is the normalized ground state solution to
\begin{align}
\begin{cases}
-\Delta u+ \mu u=\left(\int_{\mathbb{R}^N}\frac{|u|^p}{|x-y|^{N-\alpha}}dy\right)u, \quad \text{in}\ \mathbb{R}^N,\\
\int_{\mathbb{R}^N}|u|^2=\ell^{-1}.
\end{cases}
\end{align}Remark 1.3. The results in theorems 1.1 and 1.2 hold for 
$\ell=1$.
For the nonlinear Schrödinger equation
\begin{align}
\begin{cases}
-\Delta u+\lambda V(x) u+\mu u=|u|^{p-2}u\ \text{in}\ \mathbb{R}^N\\
\int_{\mathbb{R}^N}|u|^2=\omega,
\end{cases}
\end{align} the existence and orbital stability of the normalized solutions have been studied in [Reference Bartsch, Molle, Rizzi and Verzini3, Reference Cazenave and Lions6, Reference Ikoma and Miyamoto24, Reference Jeanjean and Lu25, Reference Tang, Zhang, Zhang and Zhou43, Reference Zhang and Zhang53]. Among them, Tang et al. [Reference Tang, Zhang, Zhang and Zhou43] studied the existence of normalized multibump solutions to problem (1.7) when V satisfies 
$(V_1)-(V_3)$ by using a penalized argument [Reference Zhang and Zhang53] and combining with the unique of solutions to problem (1.7) satisfying 
$V(x)=0$. However, we have to point out that when 
$V(x)=0$, problem (1.1) has a unique solution only if 
$p=\alpha=2$, and in other cases the solution is not unique due to the existence of a non-local term, that is, a Hartree term. More importantly, the penalization technique depends on the decay property of linear equations outside small balls. To this aim, it is crucial for the exponential decay of solutions. From [Reference Moroz and Van Schaftingen36], we see that the solutions decay exponentially only if 
$p\in \left[2,\frac{N+\alpha}{N-2}\right)$. Hence, the penalization technique in [Reference Zhang and Zhang53] is not valid when p < 2.
According to 
$(V_1)$, we known that 
$\Omega=\cup_{i=1}^{\ell} \Omega_{i}(\ell\geq2)$ is composed of several disjoint bounded domains 
$\{\Omega_i\}^{\ell}_{i=1}$ and each point in 
$\Omega_{i}$ is a degenerate minimizer of V. Hence, it is quite nontrivial to prove the existence of normalized multibump solutions to the limit problem (1.5) that do not vanish over multiple components of 
$\Omega_i$. If we consider the L 2-constrained problem on each component 
$\Omega_i$ separately, then the associated Lagrange multipliers may be different. Moreover, due to the interaction in the L 2-constraint, it is difficult to guarantee the nonvanishing property of gauge resolution for multiple given components.
Motivated by [Reference Byeon and Jeanjean4], the positive ground state solutions to the following limiting Hartree system
\begin{align}
\begin{cases}
-\Delta v_i+\mu v_i=\left(\int_{\mathbb{R}^N}\frac{|v_i(y)|^p}{|x-y|^{N-\alpha}}\right)|v_i(x)|^{p-2}v_i\quad \text{in}\ \mathbb{R}^N,\\
v_i(x) \gt 0,\ \lim\limits_{x\rightarrow \infty}v_i(x)=0,\ i=1,2,\dots,\ell,\\
|v_i|^2_2=s_i,\ \sum_{i=1}^{\ell} s_i =1,
\end{cases}
\end{align}can be used as the building blocks to construct the multi-bump solutions of (1.1). However, the existence of positive ground state solutions is unknown in the literature. Unlike [Reference Byeon and Jeanjean4], µ is unknown in (1.8) under the L 2-constraint. Therefore, we cannot use their method to prove the existence of positive ground state solutions to problem (1.8). Inspired by [Reference Luo33, Reference Ye51], combining some analysis tricks, we show the existence of positive ground state solutions to system (1.8), see Lemma 2.3.
In this paper, we adopt a new approach inspired by [Reference Cingolani and Tanaka9], where we obtain an operator τɛ defined by Proposition 4.1. It has similar properties to the tail minimization operator via deformation parameters in [Reference Cingolani and Tanaka9]. Using the operator τɛ in combination with the standard deformation argument, we can directly deal with the problem of 
$J_{\lambda, \varepsilon}$ in 
$X_{\varepsilon}^d$ (see Section 3) without introducing the penalization technique or requiring uniform exponential decay of the solutions, so that we can deal with the case of 
$p\in\left(\frac{N+\alpha}{N}, \frac{N+\alpha}{N-2}\right)\setminus\left\{\frac{N+\alpha+2}{N}\right\}$.
The paper is organized as follows. In Section 2, the existence of the positive ground state solutions to (1.8) will be given. Following, an 
$\varepsilon,\lambda$-dependent concentration-compactness (see Proposition 3.3) will be established in Section 3. In Section 4, we develop a deformation parameter for the tail of 
$u\in X_{\varepsilon}^d$, that is, the outside of the small balls, to construct a map τɛ. In Section 5, we use the deformation lemma to search for critical points of 
$J_{\lambda,\varepsilon}$ restricted to Mɛ and give the proof of Theorems 1.1 and 1.2.
2. Preliminary results
In this paper, for convenience, we will use the following notations.
•
$L^q(\mathbb{R}^N) (1\leq q\leq\infty)$ is the Lebesgue space with norm
$|u|_q=\left(\int_{\mathbb{R}^N}|u|^qdx\right)^{\frac{1}{q}}$.•
$H^1(\mathbb{R}^N)$ is the Sobolev space with norm
$\|u\|=\left(\int_{\mathbb{R}^N}|\nabla u|^2+|u|^2dx\right)^{\frac{1}{2}}$.•
$B_R(y)=\{x\in \mathbb{R}^N:\ |x-y| \lt R\}$.•
$o_n(1), o_{\varepsilon}(1)$ indicate a generic infinitesimal as
$n\rightarrow\infty$,
$\varepsilon\rightarrow0$.
2.1. Estimates for non-local term
Since we will be dealing with a non-local type problem (1.1), we introduce the two inequalities, which will be frequently used in this article.
(1) (Hardy-Littlewood-Sobolev inequality (see [Reference Lieb and Loss29])) Let
$q,r \gt 1$ and
$\alpha\in (0,N)$ with
$1/q+1/r=\frac{N+\alpha}{N}$. Then there exists a constant
$C=(N,\alpha,q,r) \gt 0$ such that
\begin{align*}
\left|\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{f(x)g(y)}{|x-y|^{N-\alpha}}dx\right|\leq C|f|_q|g|_r,
\end{align*}for all
$f\in L^q(\mathbb{R}^N)$ and
$g\in L^r(\mathbb{R}^N)$.(2) ([Reference Cingolani and Tanaka9]) Assume
$q,r \gt 1$ and
$\alpha\in (0,N)$ with
$1/q+1/r \lt \frac{N+\alpha}{N}$. Then for R > 0, there exists a constant
$D_R=D_R(N,\alpha,p,r) \gt 0$ such that
$\lim\limits_{R\rightarrow\infty}D_R=0$ and
\begin{align*}
\left|\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{f(x)g(y)}{|x-y|^{N-\alpha}}dx\right|\leq D_R|f|_q|g|_r
\end{align*}for all
$f\in L^q(\mathbb{R}^N)$ and
$g\in L^r(\mathbb{R}^N)$ with
$dist(supp f, suppg)\geq R$.
2.2. Limits problems
From [Reference Byeon and Jeanjean4], we know that the normalized ground state solutions to the following Kirchhoff system
\begin{align}
\begin{cases}
-\Delta v_i+\mu v_i=\left(\int_{\mathbb{R}^N}\frac{|v_i(y)|^p}{|x-y|^{N-\alpha}}\right)|v_i(x)|^{p-2}v_i\quad \text{in}\ \mathbb{R}^N,\\
v_i(x) \gt 0,\ \lim\limits_{x\rightarrow \infty}v_i(x)=0,\ i=1,2,\dots,\ell,\\
|v_i|^2_2=s_i,\ \sum_{i=1}^{\ell} s_i =1
\end{cases}
\end{align}can be used as the building blocks to construct the multi-peak solutions of (1.1). However, the existence of positive normalized ground state solutions for system (2.1) is unknown in the literature. To study system (2.1), we first give the results about the equation
\begin{align}
\begin{cases}
-\Delta v+\mu v=\left(\int_{\mathbb{R}^N}\frac{|v(y)|^p}{|x-y|^{N-\alpha}}\right)|v(x)|^{p-2}v\quad \text{in}\ \mathbb{R}^N,\\
|v|^2_2=\omega.
\end{cases}
\end{align}Define the energy functional
\begin{align*}
I(u)=\frac{1}{2}\int_{\mathbb{R}^N}|\nabla u|^2-\frac{1}{2p}\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}
\frac{|u(x)|^p|u(y)|^p}{|x-y|^{N-\alpha}}dxdy.
\end{align*} It follows from [Reference Ye51] that for 
$\frac{N+\alpha}{N} \lt p \lt \frac{N+\alpha+2}{N}$, the minimizing problem
\begin{align*}
I_{\omega}:=\inf\left\{I(u)|\ |u|_2^2=\omega,\ |\nabla u|_2^2=\frac{Np-N-\alpha}{2p}\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u(x)|^p|u(y)|^p}{|x-y|^{N-\alpha}}dxdy\right\}
\end{align*}is achieved and
\begin{align}
I_{\omega}=\inf\left\{I(u)|\ |u|_2^2=\omega\right\} \lt 0.
\end{align} As for the L 2-supercritical case, when 
$2 \lt p \lt \min\{4,N\}$ and α = 2, Luo [Reference Luo33] proved that Iω is achieved and
\begin{align}
I_{\omega}=\inf\left\{\max\limits_{t \gt 0}I(t^{N/2}u(t\cdot))|\ |u|_2^2=\omega\right\} \gt 0.
\end{align} We have to point out that the same conclusions can be obtained when 
$p\in \left(\frac{N+\alpha+2}{N},\frac{N+\alpha}{N-2}\right)$ by using the minimization method in [Reference Luo33], that is, Iω is achieved and (2.4) still holds for each 
$p\in \left(\frac{N+\alpha+2}{N},\frac{N+\alpha}{N-2}\right)$. We will not prove it here. Since the normalized ground state solution of (2.2) is not necessarily unique, for convenience, we define
\begin{align*}
S_{\omega}&=\left\{u\in H^1(\mathbb{R}^N):\ I(u)=I_{\omega}, |u|_2^2=\omega,\ |\nabla u|_2^2=\frac{Np-N-\alpha}{2p}\right. \\ &\quad\left.\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}
\frac{|u(x)|^p|u(y)|^p}{|x-y|^{N-\alpha}}dxdy\right\},
\end{align*} where 
$p\in \left(\frac{N+\alpha}{N}, \frac{N+\alpha}{N-2}\right)\setminus\left\{\frac{N+\alpha+2}{N}\right\}$. Let 
$\omega=\ell^{-1}$, we have
\begin{align*}
\mu_{\ell^{-1}}=\frac{2(-Np+N+\alpha+2p)\ell I_{\ell^{-1}}}{Np-N-2-\alpha}.
\end{align*} Clearly, for each 
$u\in S_{\ell^{-1}}$, 
$(\mu_{\ell^{-1}},u)$ is the normalized ground state solution of (2.2) when 
$\omega=\ell^{-1}$.
Set
\begin{align*}
\mathcal{S}^{\ell}_{+}:=\left\{s=(s_1,\cdots,s_{\ell})\in (0,1)^{\ell}|\sum_{i=1}^{\ell}s_i=1\right\}.
\end{align*} We define 
$s^0=\left(\ell^{-1},\cdots,\ell^{-1} \right)$ and 
$E(s):=\sum^{\ell}_{i=1}I_{s_i}$ for 
$s=(s_1,\cdots,s_{\ell})\in \mathcal{S}^{\ell}_{+}$. Clearly, 
$(\mu_{\ell^{-1}}, (u_{\ell^{-1}},$ 
$\cdots, u_{\ell^{-1}}))$ is the solution of system (2.1), where 
$u_{\ell^{-1}}\in S_{\ell^{-1}}$. Next, we claim that 
$(\mu_{\ell^{-1}},(u_{\ell^{-1}},$ 
$\cdots,u_{\ell^{-1}}))$ is the normalized ground state solution of (2.1).
Lemma 2.2. For each 
$s\in \mathcal{S}^{\ell}_{+}\setminus\ {s^0}$ (
$\ell\geq2$), there holds
\begin{align}
E(s) \lt E(s^0),\ \text{if} \ p\in\left(\frac{N+\alpha}{N}, \frac{N+\alpha+2}{N}\right),
\end{align}
\begin{align}
E(s) \gt E(s^0),\ \text{if} \ p\in\left(\frac{N+\alpha+2}{N},\frac{N+\alpha}{N-2}\right).
\end{align}Proof. Let 
$s=(s_1,\cdots,s_{\ell})\in \mathcal{S}^{\ell}_{+}\setminus\ \{0\}$. For each u with 
$|u|_2^2=\ell^{-1}$, we define
\begin{align*}
v_i(x):=\left(s_i\ell\right)^{\frac{2+\alpha}{2(-Np+N+\alpha+2)}}u\left(\left(s_i\ell\right)^{\frac{p-1}{-Np+N+\alpha+2}}\right).
\end{align*} Then by calculating, we have 
$|v_i|_2^2=s_i\ell|u|_2^2$, 
$|\nabla v_i|_2^2=\left(s_i\ell\right)^{\frac{-Np+N+\alpha+2p}{-Np+N+\alpha+2}}|\nabla u|^2_2$ and
\begin{align*}
\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|v_i(x)|^p|v_i(y)|^p}{|x-y|^{N-\alpha}}dxdy=\left(s_i\ell\right)^{\frac{-Np+N+\alpha+2p}{-Np+N+\alpha+2}}
\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u(x)|^p|u(y)|^p}{|x-y|^{N-\alpha}}dxdy.
\end{align*} Thus, according to the definition of 
$I_{s_i}$, we deduce 
$I_{s_i}=\left(s_i\ell\right)^{\frac{-Np+N+\alpha+2p}{-Np+N+\alpha+2}}I_{\ell^{-1}}$. Then combining with (2.3) and (2.4), we calculate that 
$E(s^0)=\ell I_{\ell^{-1}}$ satisfying 
$E(s^0) \lt 0$ if 
$p\in\left(\frac{N+\alpha}{N}, \frac{N+\alpha+2}{N}\right)$ and 
$E(s^0) \gt 0$ if 
$p\in\left(\frac{N+\alpha+2}{N}, \frac{N+\alpha}{N-2}\right)$. Moreover, it holds
\begin{align}
E(s)=&\sum_{i=1}^{\ell}I_{s_i}
=E(s^0)
\sum_{i=1}^{\ell}s_i^{\frac{-Np+N+\alpha+2p}{-Np+N+2+\alpha}}(\ell)^{\frac{2p-2}{-Np+N+2+\alpha}} \nonumber\\
=&E(s^0)
\sum_{i=1}^{\ell}s_i^{1+\frac{2p-2}{-Np+N+2+\alpha}}(\ell)^{\frac{2p-2}{-Np+N+2+\alpha}}.
\end{align}Denote
\begin{align*}
d_p:=
\begin{cases}
\frac{2p-2}{-Np+N+2+\alpha},\quad &\text{if}\ p\in\left(\frac{N+\alpha}{N}, \frac{N+\alpha+2}{N}\right),\\
1+\frac{2p-2}{-Np+N+2+\alpha},\quad &\text{if}\ p\in\left(\frac{N+\alpha+2}{N}, \frac{N+\alpha}{N-2}\right).
\end{cases}
\end{align*} It follows from the Höder inequality that for 
$p\in\left(\frac{N+\alpha}{N}, \frac{N+\alpha+2}{N}\right)$,
\begin{align}
\left[\sum_{i=1}^{\ell}s_i^{1+d_p}(\ell)^{d_p}\right]^
{\frac{1}{1+d_p}}=&\left[\sum_{i=1}^{\ell}s_i^{1+d_p}
(\ell)^{d_p}\right]^{\frac{1}{1+d_p}}
\nonumber\\
\geq&\sum_{i=1}^{\ell}s_i (\ell)^{d_p/(1+d_p)}(\ell)^{-d_p/(1+d_p)}=1
\end{align} and equality holds if and only if 
$s_i=\ell^{-1}$, 
$i=1,2,\cdots,\ell$. In a similar way, for 
$p\in\left(\frac{N+\alpha+2}{N}, \frac{N+\alpha}{N-2}\right)$, we also obtain
\begin{align}
\left[\sum_{i=1}^{\ell}(\ell)^{1+d_p}(s_i)^{d_p}\right]^{\frac{1}{1+d_p}}\geq 1
\end{align} and equality holds if and only if 
$s_i=\ell^{-1}$, 
$i=1,2,\cdots,\ell$. According to (2.7)–(2.9), we complete the proof.
Lemma 2.3. If 
$(\mu, (u_1,\cdots,u_{\ell}))$ is the solution of system (2.1), then we have 
$\sum_{i=1}^{\ell}I(u_i)\geq E(s^0)$. Moreover, equalities hold if and only if 
$u_i\in S_{\ell^{-1}}$ for each 
$i=1,\cdots,\ell$.
Proof. Let 
$|u_i|^2_2=s_i$, 
$\sum\limits_{i=1}^{\ell}s_i=1$ for each 
$i=1,\cdots,\ell$. For convince, we denote
\begin{align*}
A=\frac{2(-Np+N+\alpha+2p)}{-Np+N+2+\alpha}.
\end{align*} If 
$p\in\left(\frac{N+\alpha}{N},\frac{N+\alpha+2}{N}\right)$, then A > 0 and a direct calculation that
\begin{align}
\mu=-A\frac{I(u_i)}{s_i}\leq -A \frac{I_{s_i}}{s_i}=-A\left(s_i\ell\right)^{\frac{-Np+N+\alpha+2p}{-Np+N+\alpha+2}}\frac{I_{\ell^{-1}}}{s_i}
\end{align}and
\begin{align}
\mu=-A\sum_{i=1}^{\ell}I(u_i) \leq-A\sum_{i=1}^{\ell}I_{s_i}=-AE(s).
\end{align} where 
$i=1,\cdots,\ell$. If 
$s_i=\ell^{-1}$ for each 
$i=1,\cdots,\ell$, by 
$I(u_i)\geq I_{\ell^{-1}}$, we easily obtain 
$E(s)\geq E(s^0)$ and equalities holds if and only if 
$I(u_i)= I_{\ell^{-1}}$, i.e., 
$u_i\in S_{\ell^{-1}}$. Otherwise, there exists i 0 such that 
$s_{i_0} \lt \ell^{-1}$. Combining with (2.10) and (2.11), we obtain that 
$\mu=-A\sum_{i=1}^{\ell}I(u_i)\leq -A s_{i_0}^{\frac{2p-2}{-Np+N+\alpha+2}}\ell^{\frac{-Np+N+\alpha+2p}{-Np+N+\alpha+2}}I_{\ell^{-1}} \lt -A \ell I_{\ell^{-1}}$. This implies 
$\sum\limits_{i=1}^{\ell}I(u_i) \gt \ell I_{\ell^{-1}}=E(s^0)$. If 
$p\in \left(\frac{N+\alpha+2}{N},\frac{N+\alpha}{N-2}\right)$, then by Lemma 2.2, we have 
$\sum\limits_{i=1}^{\ell}I(u_i)
\geq\sum\limits_{i=1}^{\ell}I_{s_i}=E(s)\geq E(s^0)$.
3. Concentration-compactness type results
Without loss of generality, we may suppose that 
$L=\{1,2,\cdots,\ell\}$ and give the following notations: For 
$i=1,\cdots, \ell$, 
$\varepsilon\in(0,1)$, and t > 0,
\begin{align*}
&\Omega_{i,\varepsilon}:=\{x\in\mathbb{R}^N|\varepsilon x\in \Omega_i\}, \quad \!\Omega_{L,\varepsilon}:=\{x\in\mathbb{R}^N|\varepsilon x\in \Omega_L\}, \!\!\quad\Omega_{\varepsilon}:=\{x\in\mathbb{R}^N|\varepsilon x\in \Omega\},\\
&\Omega^{t}_{i,\varepsilon}:=\{x\in\mathbb{R}^N|dist(\varepsilon x,\Omega_i) \lt t\},\quad \Omega^{t}_{L,\varepsilon}:=\{x\in\mathbb{R}^N|dist(\varepsilon x,\Omega_L) \lt t\}.
\end{align*} If follows from 
$(V_1)$ that we may assume that
\begin{align*}
r:=\frac{1}{20}\min\limits_{1\leq\ell\leq\ell_0}dist(\overline{\Omega}_i,\overline{\Omega\setminus\Omega_i}) \gt 0.
\end{align*}Define
\begin{equation*}X_\varepsilon=\left\{\sum_{i=1}^\ell(\xi_\varepsilon w_i)(\cdot-x_i/\varepsilon)\vert\;x_i\in\Omega_i,\;w_{i}\in S_{\ell^{-1}}\right\},\end{equation*} where 
$\xi_{\varepsilon}(\cdot)=\xi(\varepsilon\cdot)$ with 
$\xi\in C^{\infty}_0(B_{2r}(0), [0,1])$ satisfying ξ = 1 in 
$B_r(0)$ and 
$|\nabla \xi|\leq 2/r$. Note that 
$supp\xi_{\varepsilon}(\cdot-x_i/\varepsilon)\cap\cup_{j\neq i}\Omega_{j,\varepsilon}=\emptyset$ for 
$x_i\in\overline{\Omega}_{i}$. Set
\begin{align*}
d_{\ell}:=\min\left\{1,\left\{\left(\frac{1}{2A_p}\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u(x)|^p|u(y)|^p}
{|x-y|^{N-\alpha}}\right)^{\frac{1}{2p}}:\ u\in S_{\ell^{-1}}\right\}\right\},
\end{align*} where Ap is the Sobolev embedding constant from 
$H^1(\mathbb{R}^N)$ to 
$L^p(\mathbb{R}^N)$. For sufficiently small ɛ > 0, we will find a solution to (1.1) in the following set
\begin{align*}
X_{\varepsilon}^d=\left\{u\in M_{\varepsilon}|\ \inf\limits_{v\in X_{\varepsilon}}\|u-v\|_{\varepsilon}\leq d\right\}.
\end{align*} where 
$0 \lt d \lt d_{\ell}$.
Remark 3.1. For 
$0 \lt d \lt d_{\ell}$, we have 
$\Arrowvert u\Arrowvert\leq\limsup_{\varepsilon\rightarrow0^+}\Arrowvert u\Arrowvert_\varepsilon\leq\ell\max_{i=1,\cdots,\ell}\Arrowvert w_{i}\Arrowvert+1$ for 
$u\in X_{\varepsilon}^d$. Here, 
$w_{i}\in S_{\ell^{-1}}$ is defined as in Xɛ. Moreover,
\begin{align*}
\sup\{\|u\|_{\varepsilon}| \varepsilon\in(0,1), u\in X_{\varepsilon}^d\} \lt +\infty.
\end{align*} In fact, since 
$d_{\ell}\leq1$ and 
$\|V_{\varepsilon}\|_{L^{\infty}(\Omega_{L,\varepsilon}^{2r}\setminus\Omega_{L,\varepsilon})}
=\|V\|_{L^{\infty}
(\Omega_{L}^{2r}\setminus\Omega_{L})} \lt +\infty$, we have
\begin{equation*}\Arrowvert\xi_\varepsilon w_{i}\Arrowvert\rightarrow\Arrowvert w_{i}\Arrowvert\;\text{and}\;\sup_{x_\varepsilon\in{\overline\Omega}_{L,\varepsilon}}\int_{\mathbb{R}^N}\varepsilon^2V_\varepsilon(x)\vert(\xi_\varepsilon u(x-x_\varepsilon))\vert^2\rightarrow0\;\text{as}\;\varepsilon\rightarrow0.\end{equation*}So the conclusion holds.
Define
\begin{align*}
G(u):=\frac{1}{2}\int_{\mathbb{R}^N}|u|^2,\ \text{for}\ u\in E_{\varepsilon}.
\end{align*} Then we give an 
$\varepsilon,\lambda$-dependent concentration-compactness type results, which will be useful to develop deformation theory in Section 5.
Proposition 3.2. For any 
$\varepsilon\in(0,1)$ fixed, the Palais–Smale condition holds for 
$J_{\lambda,\varepsilon}$ in the set 
$X_{\varepsilon}^d$. That is, 
$\lambda_{\varepsilon}^{\ast} \gt 0$ exists such that for any 
$\lambda \gt \lambda_{\varepsilon}^{\ast}$, if 
$\{u_n\}\subset X_{\varepsilon}^d$ satisfies
\begin{align}
&\lim\sup_{n\rightarrow \infty} J_{\lambda, \varepsilon}(u_n)\leq E(s^0),
\end{align}
\begin{align}
&\lim_{n\rightarrow\infty}\inf_{\mu_n\in \mathbb{R}}\|J'_{\lambda, \varepsilon}(u_n)+\mu_n G'(u_n)\|_{E_{\varepsilon}^*}=0,
\end{align}then un has a strong convergent subsequence in Eɛ.
Proof. From Remark 3.1, we easily obtain that 
$\{u_n\}$ is bounded in Eɛ. Since
\begin{align*}
J_{\lambda,\varepsilon}(u_n)+\frac{1}{2p}\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_n|^p|u_n|^p}{|x-y|^{N-\alpha}}
=\frac{1}{2}|\nabla u_n|^2_2+\frac{1}{2}\int_{\mathbb{R}^N}\lambda \varepsilon^2 V_{\varepsilon}(x) u_n^2+o_n(1),
\end{align*}combining with the Hardy-Little-wood-Sobolev inequality and the assumption (3.1), we obtain that for large n,
\begin{align*}
\frac{1}{2}\int_{\mathbb{R}^N}\lambda \varepsilon^2 V_{\varepsilon}(x) u_n^2\leq J_{\lambda,\varepsilon}(u_n)+C\|u_n\|^{2p}\leq C_1
\end{align*}and
\begin{align*}
|\mu_n|=&|\mu_n G'(u_n)u_n|=|J'_{\lambda,\varepsilon}(u_n)u_n+o_n(1)|\\
\leq& C\left(\|u_n\|_{\varepsilon}^2+\|u_n\|_{\varepsilon}^{2p}+\int_{\mathbb{R}^N}\lambda \varepsilon^2 V_{\varepsilon}(x) u^2\right)\leq C_1,
\end{align*} where the constant 
$C_1 \gt 0$ independent of ɛ and λ. Since 
$\{u_n\}$ is bounded in Eɛ, we assume, up to a subsequence, 
$u_n\rightharpoonup u_0$ in Eɛ.
Next, we claim 
$u_n\rightarrow u_0$ in Eɛ. We first claim that 
$\lambda_{\varepsilon}^{\ast} \gt 0$ exists such that for any 
$\lambda \gt \lambda_{\varepsilon}^{\ast}$ and δ > 0, R > 1 exists such that
\begin{align}
\lim\sup_{n\rightarrow \infty}\int_{\mathbb{R}^N\backslash B_{2R}(0)}\left( |\nabla u_n|^2+\left(\frac{\lambda}{2} \varepsilon^2 V_{\varepsilon}(x) +1\right)u^2_n\right) \lt \delta.
\end{align} we fix a cut-off function 
$\overline{\psi}\in C^{\infty}(\mathbb{R}^N, [0,1])$ satisfying 
$\overline{\psi}(x)=1$ in 
$|x|\geq2$, 
$\overline{\psi}(x)=0$ in 
$|x|\leq1$ and 
$|\nabla \overline{\psi}(x)|\leq2$ for all 
$x\in \mathbb{R}^N$. For convenience, set 
$\overline{\psi}_{R}(\cdot)=\overline{\psi}(\cdot/R)$. Then we text (3.2) by 
$\overline{\psi}_{R}^2u_n$ to get
\begin{align}
&\int_{\mathbb{R}^N}\left[ |\nabla (\overline{\psi}_R u_n)|^2+(\lambda \varepsilon^2 V_{\varepsilon}(x) +\mu_n)(\overline{\psi} _Ru_n)^2\right] \nonumber\\
=&\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{\overline{\psi}_{R}^2|u_n(x)|^p|u_n(y)|^p}{|x-y|^{N-\alpha}}+\int_{\mathbb{R}^N}|\nabla \overline{\psi}_{R}|^2u_n^2+o_n(1).
\end{align}Using the two inequalities of Proposition 2.1, we derive that
\begin{align*}
\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{\overline{\psi}^2_{R}|u_n(x)|^p|u_n(y)|^p}{|x-y|^{N-\alpha}}
\leq &\int_{\mathbb{R}^N\setminus B_{2R}(0)}\int_{\mathbb{R}^N\setminus B_{2R}(0)}\frac{|u_n(x)|^p|u_n(y)|^p}{|x-y|^{N-\alpha}}\\ &+
2\int_{R \lt |x| \lt 2R}\int_{\mathbb{R}^N}\frac{\overline{\psi}^2_R(x)|u_n(x)|^p|u_n(y)|^p}{|x-y|^{N-\alpha}}\\
&+\int_{B_R(0)}\int_{\mathbb{R}^N\setminus B_{2R}(0)}\frac{|u_n(x)|^p|u_n(y)|^p}{|x-y|^{N-\alpha}}\\
\leq& C|u_n|^{2p}_{L^{2Np/(N+\alpha)}(\mathbb{R}^N\setminus B_{2R}(0))}+D_R\\
\leq& C\|\overline{\psi}_R u_n\|^{2}| u_n|^{2(p-1)}_{L^{2Np/(N+\alpha)}(\mathbb{R}^N\setminus B_{2R}(0))}+D_R
\end{align*}and
\begin{align*}
\int_{\mathbb{R}^N}|\nabla \overline{\psi}_{R}|^2u_n^2=\int_{R \lt |x| \lt 2R}|\nabla \overline{\psi}_{R}|^2u_n^2\leq\frac{C}{R}.
\end{align*} where C is a positive constant independent of ɛ, R, and large n, 
$\lim\limits_{R\rightarrow\infty}D_R=0$. Then combining with (3.4), we obtain
\begin{align}
&\int_{\mathbb{R}^N}\left[ |\nabla (\overline{\psi}_R u_n)|^2+(\lambda \varepsilon^2 V_{\varepsilon}(x) +\mu_n)(\overline{\psi}_Ru_n)^2\right] \nonumber\\
\leq& C_2\|\overline{\psi}_R u_n\|^{2}| u_n|^{2(p-1)}_{L^{2Np/(N+\alpha)}(\mathbb{R}^N\setminus B_{2R}(0))}+D_R+\frac{C_2}{R}+o_n(1),
\end{align} where C 2 is a constant independent of 
$\varepsilon, \lambda, R$ and large n. According to the condition 
$(V_3)$, we denote
\begin{align*}
K_{0}(R):=\left(\mathbb{R}^N\setminus B_{R}(0)\right)\cap\{x\in \mathbb{R}^N| V_{\varepsilon}(x)\leq M_0\}
\end{align*}and
\begin{align*}
L_{0}(R):=\left(\mathbb{R}^N\setminus B_{R}(0)\right)\cap\{x\in \mathbb{R}^N| V_{\varepsilon}(x) \gt M_0\}.
\end{align*} Then we deduce that for every 
$\varepsilon, \lambda \gt 0$,
\begin{align}
\lim_{R\rightarrow \infty}meas\ K_{0}(R)=0,\quad
\lim_{R\rightarrow \infty}\limsup_{n\rightarrow\infty}\int_{K_{0}(R)}|u_n|^{\frac{2Np}{N+\alpha}}=0
\end{align}and
\begin{align}
&\lim_{R\rightarrow \infty}\limsup_{n\rightarrow\infty}\int_{K_{0}(R)}\left(\lambda \varepsilon^2V_{\varepsilon}(x)+|\mu_n|\right)(\overline{\psi}_Ru_n)^2 \nonumber\\
\leq&\lim_{R\rightarrow \infty}\limsup_{n\rightarrow\infty}\int_{K_{0}(R)}\left(\lambda\varepsilon^2M_0+C_1\right)(\overline{\psi}_Ru_n)^2=0.
\end{align}It follows from the Gagliardo–Nirenberg inequality that
\begin{align}
\int_{L_{0}(R)}|u_n|^{\frac{2Np}{N+\alpha}}\leq&C_3|\nabla u_n|_{L^2(L_{0}(R))}^{\frac{(Np-N-\alpha)N}{N+\alpha}}
|u_n|_{L^2(L_{0}(R))}^{\frac{(-Np+N+\alpha+2p)N}{N+\alpha}} \nonumber\\
\leq&C_3\left(M_0^{-1}\lambda^{-1}\varepsilon^{-2}\right)^{\frac{(-Np+N+\alpha+2p)N}{N+\alpha}},
\end{align} where C 3 is a positive constant independent of ɛ, λ and large n. Using (3.6) and (3.8), and fixing 
$\lambda_{\varepsilon}^{\ast}$ large enough such that for 
$\lambda \gt \lambda_{\varepsilon}^{\ast}$,
\begin{align}
\lim_{R\rightarrow\infty}\limsup_{n\rightarrow\infty}|u_n|^{2(p-1)}_{L^{\frac{2Np}{N+\alpha}}(\mathbb{R}^N\setminus B_{2R}(0))} \lt (2C_2)^{-1},\ \frac{\lambda}{2}\varepsilon^2 M_0-C_1 \gt 1.
\end{align}
\begin{align}
&\limsup\limits_{R\rightarrow\infty}\limsup\limits_{n\rightarrow\infty}\int_{\mathbb{R}^N}(\lambda \varepsilon^2 V_{\varepsilon}(x) +\mu_n)(\overline{\psi}_Ru_n)^2\nonumber\\
&\quad=\limsup\limits_{R\rightarrow\infty}\limsup\limits_{n\rightarrow\infty}\int_{L_{0}(R)}(\lambda \varepsilon^2 V_{\varepsilon}(x) +\mu_n)(\overline{\psi}_Ru_n)^2 \nonumber\\
&\quad\geq \limsup\limits_{R\rightarrow\infty}\limsup\limits_{n\rightarrow\infty}\int_{L_{0}(R)}\left(\frac{1}{2}\lambda \varepsilon^2 V_{\varepsilon}(x) +1\right)(\overline{\psi}_Ru_n)^2 \nonumber\\
&\quad=\limsup\limits_{R\rightarrow\infty}\limsup\limits_{n\rightarrow\infty}\int_{\mathbb{R}^N}\left(\frac{1}{2}\lambda \varepsilon^2 V_{\varepsilon}(x) +1\right)(\overline{\psi}_Ru_n)^2.
\end{align}Thus, it follows from (3.5), (3.9), and (3.10) that
\begin{align*}
\limsup_{R\rightarrow\infty}\limsup_{n\rightarrow\infty}\int_{\mathbb{R}^N}\left[ |\nabla (\overline{\psi}_R u_n)|^2+\left(\frac{\lambda }{2} \varepsilon^2 V_{\varepsilon}(x) +1\right)(\overline{\psi}_Ru_n)^2\right]=0,
\end{align*}which implies (3.3). Then by the Rellich theorem and (3.3), we deduce that
\begin{align}
\int_{\mathbb{R}^N}V_{\varepsilon}(x)|u_n-u_0|^2\rightarrow0\ \text{and}\ u_n\rightarrow u_0\ \text{in}\ L^q(\mathbb{R}^N), q\in \left[2,\frac{N+\alpha}{N-2}\right).
\end{align} According to this, we text (3.2) by 
$u_n-u_0$ and get
\begin{align*}
\int_{\mathbb{R}^N}\nabla u_n\nabla(u_n-u_0)=&-\int_{\mathbb{R}^N}\lambda \varepsilon^2 V_{\varepsilon}(x)u_n(u_n-u_0)-\mu_n(u_n-u_0)u_n
\\ &+\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_n|^pu_n^{p-1}(u_n-u_0)}{|x-y|^{N-\alpha}}\\
=&o_n(1),
\end{align*}which implies
\begin{align}
\int_{\mathbb{R}^N}|\nabla u_n-\nabla u_0|^2\rightarrow 0.
\end{align} By (3.11)–(3.12), we easily get 
$u_n\rightarrow u_0$ in Eɛ. Thus, we complete the proof.
Since 
$meas(\Omega_{L,\varepsilon}^{8r}\setminus\Omega_{L,\varepsilon}) \lt +\infty$ and 
$V_{\varepsilon} \gt 0$ a.e. in 
$\Omega_{L,\varepsilon}^{8r}\setminus\Omega_{L,\varepsilon}$, we choose 
$\alpha_{\varepsilon}\in(0,1)$ such that
\begin{align}
meas(A^{\varepsilon}:=\{x\in\Omega_{L,\varepsilon}^{8r}\setminus\Omega_{L,\varepsilon}:\ V_{\varepsilon}(x) \lt \alpha_{\varepsilon}\})\rightarrow0\ \text{as}\ \varepsilon\rightarrow0.
\end{align}Then we have the following result.
Proposition 3.3. For sufficiently small d > 0, 
$J_{\lambda,\varepsilon}$ has the following property: let 
$\{\varepsilon_n\}$, 
$\{\lambda_n\}$, 
$\{\mu_n\}\subset\mathbb{R}$, and 
$\{u_n\}\subset X_{\varepsilon}^d$ satisfying that
\begin{align}
&\varepsilon_n\rightarrow 0^+,\quad \lambda_n \varepsilon_n^2\alpha_{\varepsilon_n}\rightarrow+\infty,\nonumber\\
&\lim_{n\rightarrow\infty}J_{\lambda_n,\varepsilon_n}(u_n)\leq E(s^0),
\end{align}
\begin{align}
&\|J'_{\lambda_n,\varepsilon_n}(u_n)+\mu_nG'(u_n)\|_{E^*_{\varepsilon_n}}\rightarrow0,
\end{align} as 
$n\rightarrow \infty$. Then 
$\mu_n\rightarrow \mu_{\ell^{-1}}$ and for each 
$i=1,\cdots,\ell$, there exist 
$q_{n,i}\in\Omega_{i,\varepsilon_n}$, 
$u_i\in S_{\ell^{-1}}$ such that
\begin{align*}
dist (q_{n,i},\partial\Omega_{i,\varepsilon_n})\rightarrow \infty\quad \text{and}\quad \|u_n-\sum^{\ell}_{i=1}(\xi_{\varepsilon_n}u_i(\cdot-q_{n,i}))\|_{\varepsilon_n}\rightarrow0.
\end{align*}Proof. Since 
$u_n\in X_{\varepsilon}^d$, there exist 
$q_{n,i}\in\overline{\Omega}_{i,\varepsilon_n}$, 
$u_{i}\in S_{\ell^{-1}}$ and 
$h_n\in E_{\varepsilon_n}$ satisfying 
$\|h_n\|_{\varepsilon_n}\leq d$ such that 
$u_n=\sum_{i=1}^\ell(\xi_{\varepsilon_n}u_i)(\cdot-q_{n,i})+h_n$. By Remark 3.1, we obtain that 
$\|u_n\|$ and 
$\|u_n\|_{\varepsilon_n}$ are bounded for large n. Similar to Proposition 3.2, by (3.14), we also deduce that for large n,
\begin{align*}
|\mu_n|,\ \int_{\mathbb{R}^N}\lambda_n \varepsilon_n^2V_{\varepsilon_n}(x)u_n^2\leq C
\end{align*}where the constant C > 0 is independent of n. Then by (3.13),
\begin{align*}
\int_{\mathbb{R}^N\setminus\Omega_{L,\varepsilon}}|u_n|^2\leq& \int_{\mathbb{R}^N\setminus\Omega^{8r}_{L,\varepsilon}}|u_n|^2+\int_{\Omega_{L,\varepsilon}^{8r}\setminus(A^{\varepsilon_n}\cup
\Omega_{L,\varepsilon})}u_n^2+\int_{A^{\varepsilon_n}}|u_n|^2\\ \leq&
\int_{\mathbb{R}^N\setminus\Omega^{8r}_{L,\varepsilon_n}}|u_n|^2+(\lambda\varepsilon^2\alpha_{\varepsilon_n})^{-1}+|meas A^{\varepsilon_n}|^{\frac{\alpha+2}{N+\alpha}}|u_n|^2_{\frac{N+\alpha}{N-2}},
\end{align*} where 
$meas A^{\varepsilon_n}\rightarrow0$ as 
$n\rightarrow\infty$. Now we claim that 
$\lim\limits_{n\rightarrow\infty}\int_{\mathbb{R}^N\setminus\Omega^{8r}_{L,\varepsilon_n}}u_n^2=0$. Otherwise, by Lions’ Lemma, there is 
$x_{\varepsilon_n}\in \mathbb{R}^N\setminus\Omega^{8r}_{L,\varepsilon_n}$ such that
\begin{align*}
\limsup\limits_{n\rightarrow\infty}\int_{B_1(x_{\varepsilon_n})}|u_n|^2 \gt 0.
\end{align*} Since 
$|q_{n,i}-x_{\varepsilon_n}|\rightarrow\infty$, we have
\begin{equation*}\vert\sum_{i=1}^\ell(\xi_{\varepsilon_n}u_{i})(\cdot-q_{n,i})\vert_{L^2(B_1(x_{\varepsilon_n}))}^2=o_n(1).\end{equation*} It follows from 
$\Arrowvert u_n-\sum_{i=1}^\ell(\xi_{\varepsilon_n}u_{i})(\cdot-q_{n,i})\Arrowvert_{B_1(x_{\varepsilon_n})}\leq2d$ that we obtain a contradiction for sufficiently small d. Hence, we conclude that 
$\int_{\mathbb{R}^N\setminus\Omega_{L,\varepsilon_n}^{8r}}|u_n|^2=o_n(1)$.
Since 
$|\mu_n|$ and 
$\|u_n\|$ are bounded, up to a subsequence, we may assume that 
$\mu_n\rightarrow\mu_0$ in 
$\mathbb{R}$ and
\begin{align*}
u_n(\cdot+q_{n,i})\rightharpoonup v_i\quad& \text{in}\ H^1(\mathbb{R}^N)\ \text{and}\ L^q(\mathbb{R}^N), q\in\left[2,\frac{N+\alpha}{N-2}\right)\\
u_n(\cdot+q_{n,i})\rightarrow v_i\quad& \text{in}\ L^q_{loc}(\mathbb{R}^N), q\in \left[2,\frac{N+\alpha}{N-2}\right)\\
u_n(\cdot+q_{n,i})\rightarrow v_i\quad& a.e.\ \text{in}\ \mathbb{R}^N.
\end{align*}Obviously, there holds
\begin{align}
&\lim_{R\rightarrow\infty}\lim_{n\rightarrow\infty}\int_{B_R(0)}\int_{B_R(0)}\frac{|u_n(\cdot+q_{n,i})|^p
|u_n(\cdot+q_{n,i})|^p}{|x-y|^{N-\alpha}} \nonumber\\
=&
\lim_{R\rightarrow\infty}\int_{B_R(0)}\int_{B_R(0)}\frac{|v_i|^p|v_i|^p}{|x-y|^{N-\alpha}}
=\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|v_i|^p|v_i|^p}{|x-y|^{N-\alpha}}.
\end{align} On the other hand, since 
$|q_{n,i}-q_{n,j}|\rightarrow\infty$ if i ≠ j, we get that for each fixed R > 0,
\begin{equation*}\int_{B_R(0)}\int_{B_R(0)}\frac{|u_n(\cdot+q_{n,i})|^p|u_n(\cdot+q_{n,i})|^p}{|x-y|^{N-\alpha}}
\geq\int_{B_R(0)}\int_{B_R(0)}\frac{|u_i|^p|u_i|^p}{|x-y|^{N-\alpha}}-A_pd^{2p}+o_n(1).\end{equation*}For d sufficiently small, according to the above inequality and combining with (3.16), we deduce that
\begin{align}\int_{R^N}\int_{R^N}\frac{\left|v_i(x)\right|^p\left|v_i(y)\right|^p}{\left|x-y\right|^{N-\alpha}}& \geq \int_{R^N}\int_{R^N}\frac{\left|u_i\right|^p\left|u_i\right|^p}{\left|x-y\right|^{N-\alpha}}-A_pd^{2p} \nonumber \\ & \geq \frac12\int_{R^N}\int_{R^N}\frac{\left|u_i\right|^p\left|u_i\right|^p}{\left|x-y\right|^{N-\alpha}}\geq0.\end{align}By Fatou’s lemma, we have
\begin{align}
|v_i|_2\leq\lim\limits_{n\rightarrow \infty}\inf|u_n|^2=1 \lt 2.
\end{align}The following proof proceeds in four claims.
Claim 1 
$\lim\limits_{n\rightarrow\infty}dist (q_{n,i},\partial \Omega_{i,\varepsilon_n})\rightarrow \infty$.
Suppose, on the contrary, that for some i, up to a subsequence, 
$\lim\limits_{n\rightarrow \infty}dist (q_{n,i},\partial \Omega_{i,\varepsilon_n})=M\in[0,\infty)$, and for some 
$q_i\in\partial \Omega_i$, 
$\varepsilon_n q_{n,i}\rightarrow q_i$ as 
$n\rightarrow \infty$. By a transformation, we may assume that 
$e_N:=(0,\cdots,0,1)$ is inward normal to 
$\partial \Omega_i$ at qi. Denote 
$A=\{x\in \mathbb{R}^N|x_N \gt -M\}$ and let 
$B+q_{n,i}=\{x+q_{n,i}|x\in B\}$ where 
$B\subset\rm{int}(\mathbb{R}^N\setminus A)$ is a closed ball. Then we claim
\begin{align}
\lim\limits_{n\rightarrow\infty}\inf dist (B+q_{n,i},\Omega_{i,\varepsilon_n}) \gt 0.
\end{align} Through the compactness and indirect argument, we assume 
$|x+q_{n,i}-x_n|\rightarrow0$ for some 
$x\in B$ and 
$x_n\in\Omega_{i,\varepsilon_n}$. Let 
$q_n\in \partial \Omega_{i,\varepsilon_n}$ be such that 
$|q_{n,i}-q_n|=dist (q_{n,i}, \partial \Omega_{i,\varepsilon_n})\rightarrow M$. Clearly, it holds 
$\{q_{n,i}, x+q_{n,i},x_n\}\subset B_{|x|+M+1}(q_n)$ for large n. Since 
$\varepsilon_nq_n\rightarrow q_i$, we derive that 
$(q_{n,i}-q_n)\cdot e_N\rightarrow M$ and 
$\sup\{(q_n-y)\cdot e_N|y\in \Omega_{i, \varepsilon_n}\cap B_{|x|+M+1}(q_n)\}\rightarrow 0$. By calculating, we conclude that
\begin{align*}
0=\lim\limits_{n\rightarrow \infty}\sup(x+q_{n,i}-x_n)\cdot e_N=x\cdot e_N+M+\lim\limits_{n\rightarrow\infty}\sup(q_n-x_n)\cdot e_N \lt 0,
\end{align*}which contradicts with (3.19). Since
\begin{align*}
0=&\lim\limits_{n\rightarrow\infty}|u_n|_{L^2(\mathbb{R}^N\setminus\Omega_{\varepsilon_n})}\leq \lim\limits_{n\rightarrow\infty}|u_n|_{L^2(B+q_{n,i})}\\
=&\lim\limits_{n\rightarrow\infty}|u_n(\cdot+q_{n,i})|_{L^2(B)}=|v_i|_{L^2(B)},
\end{align*} we obtain that 
$v_i\in H^1_0(A)$ by the arbitrary choice of B. Next, we text (3.15) by 
$\psi(\cdot-q_{n,i})$ for any 
$\psi\in C^{\infty}_0(A)$. Similar to (3.19), we show that 
$supp \psi(\cdot-q_{n,i})\subset\Omega_{L,\varepsilon_n}$ for large n,
\begin{align*}
o_n(1)=&\int_{\mathbb{R}^N}\nabla u_n\nabla\psi(x-q_{n,i})+\int_{\mathbb{R}^N}\lambda_n \varepsilon^2_nV_{\varepsilon_n}(x)u_n\psi(x-q_{n,i})+\mu_n\\ &\times\int_{\mathbb{R}^N}u_n\psi(x-q_{n,i})
-\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_n(y)|^p|u_n(x)|^{p-2}u_n(x)\psi(x-q_{n,i})}{|x-y|^{N-\alpha}}\\
=&\int_{\mathbb{R}^N}\nabla u_n(x+q_{n,i})\nabla\psi(x)+\mu_n\int_{\mathbb{R}^N}u_n(x+q_{n,i})\psi(x)\\
&-\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_n(y+q_{n,i})|^p|u_n(x+q_{n,i})|^{p-2}u_n(x+q_{n,i})\psi(x)}{|x-y|^{N-\alpha}}+o_n(1)
\\=&\int_{A}\nabla v_i\nabla\psi(x)+\mu_0\int_{A}v_i\psi(x)
-\int_{A}\int_{A}\frac{|v_i|^p|v_i|^{p-2}\psi(x)v_i(x)}{|x-y|^{N-\alpha}}+o_n(1).
\end{align*} According to this, we deduce that vi is a nontrivial solution to 
$-\Delta v_i+\mu_0v_i=\left(\int_{A}\frac{|v_i|^{p}}{|x-y|^{N-\alpha}}\right)|v_i|^{p-2}v_i$, which is impossible since this equation has no nontrivial solution on a half space by [Reference Moroz and Van Schaftingen35].
Claim 2. 
$\mu_0 \gt 0$ and 
$v_i \gt 0$ for 
$i=1,\cdots,\ell$.
Since 
$\lim\limits_{n\rightarrow\infty}dist (q_{n,i},\partial\Omega_{i,\varepsilon_n})=\infty$, it holds
\begin{align}
-\Delta v_i+\mu_0v_i=\left(\int_{\mathbb{R}^N}\frac{|v_i|^{p}}{|x-y|^{N-\alpha}}\right)|v_i|^{p-2}v_i,\quad \text{in}\ H^1(\mathbb{R}^N).
\end{align} By the Pohozaev identity, we deduce 
$\mu_0 \gt 0$. Now we claim 
$v_i \gt 0$. Since u > 0, we deduce that for any d < 1,
\begin{equation*}\int_{\mathbb{R}^N}\vert u_n^-\vert^2\leq\int_{\mathbb{R}^N}\left|u_n-\sum_{j=1}^\ell(\xi_{\varepsilon_n}u_i)(\cdot-q_{n,i})\right|^2=\int_{\mathbb{R}^N}\vert h_n\vert^2\leq d^2\end{equation*}From Fatou’s Lemma, one has
\begin{align}
|v_i^-|_2^2\leq\liminf\limits_{n\rightarrow\infty}|u_n^-(\cdot+q_{n,j})|^2_2 \lt d.
\end{align} By Remark 3.1, 
$\Arrowvert v_i^-\Arrowvert\leq\ell\max_{i=1,\cdots,\ell}\Arrowvert u_i\Arrowvert+d$. According to equation (3.20), we get
\begin{equation*}\int_{\mathbb{R}^N}\vert\nabla v_i^-\vert^2+\mu(v_i^-)^2=\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{\vert v_i^-\vert^p\vert v_i^-\vert^p}{\vert x-y\vert^{N-\alpha}}dxdy\leq C_\beta\vert(v_i^-)^2\vert_2^{2p}+D\beta\vert\nabla v_i^-\vert_2^2,\end{equation*} where β is chosen so that 
$D\beta=\beta(\ell\Arrowvert v_i^-\Arrowvert+1)^{2(p-1)}=1$. This implies
\begin{align}
\mu\leq C_{\beta}\left(\int_{\mathbb{R}^N}(v_i^-)^2\right)^{p-1}.
\end{align} Combining with (3.21)–(3.22), we deduce 
$v_i^-=0$ by controlling d small. Using the maximum principle, we have 
$v_i \gt 0$.
Claim 3. 
$u_n-\sum^{\ell}_{i=1}v_i\left(\cdot-q_{n,i}\right)\rightarrow 0$ in 
$L^{q}(\mathbb{R}^N)$ for any 
$q\in \left[2, \frac{N+\alpha}{N-2}\right)$.
Suppose, on the contrary, that according to the Lions’ vanishing lemma [Reference Lions31], there is a point sequence 
$\{x_n\}\subset\mathbb{R}^N$ that satisfies
\begin{align}
\limsup\limits_{n\rightarrow\infty}\int_{B_1(x_n)}\big|u_n-\sum^{\ell}_{i=1}v_i(\cdot-q_{n,i})\big|^2 \gt 0.
\end{align}For any fixed R > 0, there holds
\begin{align*}
\int_{B_R(q_{n,i})}|u_n-\sum^{\ell}_{j=1}v_j(\cdot-q_{n,j})|^2=
\int_{B_R(0)}|u_n(\cdot+q_{n,i})-v_i|^2+o_n(1)\rightarrow0.
\end{align*} Combining with (3.23), we have 
$|x_n-q_{n,i}|\rightarrow \infty$ for 
$i=1,2,\cdots,\ell$. Hence, we deduce 
$\|v_i(\cdot+x_n-q_{n,i})\|_{L^2(B_1(0))}\rightarrow0$ and
\begin{align}
\limsup_{n\rightarrow\infty}|u_n(\cdot+x_n)|_{{L^2}(B_1(0))} \gt 0.
\end{align} We assume, up to a subsequence, 
$u_n(\cdot+x_n)\rightharpoonup v_0$ in 
$H^1(\mathbb{R}^N)$. It follows from (3.24) that 
$v_0\neq0$.
For R > 0, we set 
$\psi_R=1-\overline{\psi}(\cdot/R)\in C^{\infty}_0(B_{2R}(0))$ where 
$\overline{\psi}$ is defined as Proposition 3.2. Then we test (3.15) by 
$\psi_{R,n}u_n:=\psi_{R}(\cdot-x_n)u_n$ and obtain that
\begin{align*}
o_n(1)=&\int_{\mathbb{R}^N}\left(|\nabla u_n|^2+\mu_n u_n^2+\lambda_n\varepsilon_n^2V_{\varepsilon_n}u_n^2\right)\psi_{R,n}+\int_{\mathbb{R}^N} u_n\nabla u_n\nabla \psi_{R,n}\\ &-
\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_n(y)|^p|u_n(x)|^p\psi_{R,n}}{|x-y|^{N-\alpha}}\\ \geq&
\int_{\mathbb{R}^N}\left(|\nabla u_n|^2+\mu_n u_n^2\right)\psi_{R,n}-
\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_n(y)|^p|u_n(x)|^p\psi_{R,n}}{|x-y|^{N-\alpha}}-\frac{C}{R}.
\end{align*} By the Rellich theorem and 
$\mu_n\rightarrow \mu_0$,
\begin{align*}
&\int_{\mathbb{R}^N}\mu_n u_n^2\psi_{R,n}-\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_n(y)|^p|u_n(x)|^p\psi_{R,n}}{|x-y|^{N-\alpha}}\\=&
\int_{\mathbb{R}^N}\mu_n u_n(\cdot+x_n)^2\psi_{R}-\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}
\frac{|u_n(y+x_n)|^p|u_n(x+x_n)|^p\psi_{R}}{|x-y|^{N-\alpha}}\\
\rightarrow&\int_{\mathbb{R}^N}\mu_0 v_0^2\psi_{R}-\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|v_0|^p|v_0|^p\psi_{R}}{|x-y|^{N-\alpha}}.
\end{align*}Then we get that
\begin{align*}
&\int_{\mathbb{R}^N}\left(|\nabla v_0|^2+\mu_0 v_0^2\right)\psi_{R}-\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|v_0(y)|^p|v_0(x)|^p\psi_{R}(x)}{|x-y|^{N-\alpha}}\\
\leq&\limsup\limits_{n\rightarrow\infty}\int_{\mathbb{R}^N}\left(|\nabla u_n|^2+\mu_n u_n^2\right)\psi_{R}-\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_n(y)|^p|u_n(x)|^p\psi_{R,n}}{|x-y|^{N-\alpha}}\\
\leq&\frac{C}{R}.
\end{align*} Passing to the limit 
$R\rightarrow \infty$ in the above equation yields that it holds
\begin{align*}
\int_{\mathbb{R}^N}|\nabla v_0|^2+\mu_0 v_0^2\leq\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|v_0(x)|^p|v_0(y)|^p}{|x-y|^{N-\alpha}}.
\end{align*}Arguing in a similar way to (3.17) and (3.22) yields
\begin{equation*}\vert v_i\vert_2\geq\vert u_i\vert_2-d=\ell^\frac{-1}2-d\end{equation*}and
\begin{align*}
\mu_0\leq C_{\beta}\left(\int_{\mathbb{R}^N}| v_0|^2\right)^{p-1}.
\end{align*}Then when d sufficiently small,
\begin{align*}
\sum\limits_{i=1}^{\ell}|v_i|^2_2+|v_0|^2_2\geq 1-2d\ell^{1/2}+\ell d^2+\ell|v_0|^2_2 \gt 1.
\end{align*}This is impossible since
\begin{align*}
1=\lim\limits_{n\rightarrow \infty}|u_n|^2_2\geq&\lim\limits_{R\rightarrow\infty}\limsup\limits_{n\rightarrow\infty}\left(\sum\limits_{i=1}^{\ell}
\|u_n(\cdot+q_{n,i})\|^2_{L^2(B_R(0))}
\,{+}\,\|u_n(\cdot+x_{n})\|^2_{L^2(B_R(0))}\!\right)\\=&\sum\limits_{i=1}^{\ell}|v_i|^2_2+|v_0|^2_2.
\end{align*}Thus, claim 3 is proved.
Claim 4. 
$\|u_n-\sum\limits_{i=1}^{\ell}\left(\xi_{\varepsilon_n}v_i\right)\left
(\cdot-q_{n,i}\right)\|_{\varepsilon_n}\rightarrow0$.
Let 
$R_n:=\frac{1}{2}\lim\limits_{n\rightarrow\infty}\min\limits_{1\leq i\leq \ell}dist \left(q_{n,i},\partial\Omega_{i,\varepsilon_n}\right)$. According to Claim 1, we have 
$R_n\rightarrow \infty$. Set 
$\widetilde{\psi}_n=1-\overline{\psi}(\cdot/R_n)\in C^{\infty}_0(B_{2R_n}(0))$. It is easy to check that 
$\|\left(\xi_{\varepsilon_n}v_i\right)\left(\cdot-q_{n,i}\right)-(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})\|_{\varepsilon_n}\rightarrow0$. Thus, we only need to claim
\begin{align}
\|u_n-\sum\limits_{i=1}^{\ell}(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})\|_{\varepsilon_n}\rightarrow 0.
\end{align} Noting that 
$supp\ \widetilde{\psi}_n(\cdot-q_{n,i})\subset B_{2R_n}(q_{n,i})\subset\subset\Omega_{i,\varepsilon_n}$ and (3.15) is tested by 
$u_n-\sum\limits_{i=1}^{\ell}(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})$, as well as claim 3, we obtain
\begin{align*}
&\int_{\mathbb{R}^N}\nabla u_n\nabla \left(u_n-\sum_{i=1}^{\ell}(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})\right)\\ &+\int_{\mathbb{R}^N}(\lambda_n\varepsilon_n^2V_{\varepsilon_n}+\mu_n)u_n\left(u_n-\sum\limits_{i=1}^{\ell}
(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})\right)\leq o_n(1).
\end{align*} Due to 
$\widetilde{\psi}_nv_j\rightarrow v_j$ and 
$u_n(\cdot+q_{n,j})\rightharpoonup v_j$ in 
$H^1(\mathbb{R}^N)$ for each 
$j=1,\cdots,\ell$, and 
$|q_{n,j}-q_{n,i}|\rightarrow\infty$ for i ≠ j, it follows that
\begin{align*}
&\int_{\mathbb{R}^N}\nabla (\widetilde{\psi}_nv_j)(\cdot-q_{n,j})\nabla\left(u_n-\sum\limits_{i=1}^{\ell}(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})\right)
\\=&\int_{\mathbb{R}^N}\nabla v_j\nabla\left(u_n(\cdot+q_{n,i})-v_j-\sum^{\ell}_{i=1,i\neq j}v_i(\cdot+q_{n,j}-q_{n,i})\right)+o_n(1)\\=&\int_{\mathbb{R}^N}\nabla v_j\nabla\left(u_n(\cdot+q_{n,j})-v_j\right)+o_n(1)=o_n(1).
\end{align*} Moreover, it follows from 
$V_{\varepsilon_n}=0$ in 
$supp\ \varphi_{n}(\cdot-q_{n,j})\subset\Omega_{\varepsilon_n}$ that
\begin{align*}
&\int_{\mathbb{R}^N}(\lambda_n\varepsilon_n^2V_{\varepsilon_n}+\mu_n)\left(u_n-\sum\limits_{i=1}^{\ell}
(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})\right)(\widetilde{\psi}_nv_j)(\cdot-q_{n,j})\\=&\int_{\mathbb{R}^N}\mu_n\left(u_n-\sum\limits_{i=1}^{\ell}
(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})\right)(\widetilde{\psi}_nv_j)(\cdot-q_{n,j})=o_n(1).
\end{align*}Thus, from the above results, we conclude
\begin{align*}
&\int_{\mathbb{R}^N}\left|\nabla\left(u_n-\sum_{i=1}^{\ell}(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})\right)\right|^2+(\lambda_n\varepsilon_n^2V_{\varepsilon_n}+\mu_n)
\left|u_n-\sum\limits_{i=1}^{\ell}
(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})\right|^2\\ &\quad\leq o_n(1).
\end{align*}This proves claim 4.
Now we are ready to complete the proof. By claim 4, we have
\begin{align*}
\sum\limits_{i=1}^{\ell}|v_i|^2_2=\lim_{n\rightarrow\infty}\left|\sum\limits_{i=1}^{\ell}
(\xi_{\varepsilon_n}v_i)(\cdot-q_{n,i})\right|^2_2=1
\end{align*}and together with Lemma 2.3,
\begin{align*}
\lim\limits_{n\rightarrow\infty}J_{\lambda,\varepsilon}(u_n)=&\lim\limits_{n\rightarrow\infty}\left(\frac{1}{2}\left| \nabla \sum_{i=1}^{\ell}(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})\right|^2_2-\frac{1}{2p}\right.\\
&\times\left.\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}
\frac{\left|\left(\sum_{i=1}^{\ell}(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})\right)\right|^p\left|
\left(\sum_{i=1}^{\ell}(\widetilde{\psi}_nv_i)(\cdot-q_{n,i})\right)\right|^p}{|x-y|^{N-\alpha}}\right)\\
=&\sum_{i=1}^{\ell}\left(\frac{1}{2}| \nabla v_i|^2_2-\frac{1}{2p}\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}
\frac{|v_i|^p|v_i|^p}{|x-y|^{N-\alpha}}\right)\geq E(s^0).
\end{align*} This indicates 
$v_i\in S_{\ell^{-1}}$ and the proof is completed.
4. Deformation outside small balls
To find the critical point of 
$J_{\lambda,\varepsilon}$ in 
$X_{\varepsilon}^d$, it is crucial to control the size of u of outside of the set 
$\Omega^{8r}_{L,\varepsilon}$. We notice that the penalization technique described in [Reference Tang, Zhang, Zhang and Zhou43] can only handle the case of 
$p\geq $2, so we take a different approach to construct a special deformation flow outside a small ball, i.e., we construct a map τɛ as shown below.
Define
\begin{align*}
Y_{\varepsilon}^d=\left\{u\in E_{\varepsilon}| \inf\limits_{v\in X_{\varepsilon}}\|u-v\|_{\varepsilon}\leq d\right\}.
\end{align*}Proposition 4.1. For sufficiently small ɛ > 0, there exists a map 
$\tau_{\varepsilon}:X_{\varepsilon}^d\rightarrow Y_{\varepsilon}^d$ such that
(i) τɛ is continuous.
(ii)
$\tau_{\varepsilon}(u)=u$ if
$u(x)=0$ for all
$(\Omega_{L,\varepsilon}^{8r})^c$.(iii) For each
$u\in X_{\varepsilon}^d$,
(4.1)\begin{align}
&\tau_{\varepsilon}(u)=u(x)\ \text{for}\ x\in\Omega^{2r}_{L,\varepsilon},
\end{align}(4.2)\begin{align}
&|\nabla \tau_{\varepsilon}(u)|^2_{L^2\left((\Omega_{L,\varepsilon}^{8r})^c\right)}\leq \sigma_{\varepsilon},
\end{align}(4.3)\begin{align}
&J_{\lambda, \varepsilon}(\tau_{\varepsilon}(u))\leq J_{\lambda,\varepsilon}(u).
\end{align}
Here 
$\sigma_{\varepsilon} \gt 0$ is independent of u and satisfies 
$\lim\limits_{\varepsilon\rightarrow0}\sigma_{\varepsilon}=0$.
In what follows, we construct 
$\Gamma_{\varepsilon}(u): X_{\varepsilon}^{d}\rightarrow \mathbb{R}$ by
\begin{align*}
&\Gamma_{\varepsilon}(u)=|\nabla u|^2_{L^2((\Omega_{L,\varepsilon}^{8r})^c)},
\end{align*}which will play an important role to build τɛ.
Lemma 4.2. For any 
$u\in X_{\varepsilon}^{d}$, there exists a 
$\zeta_{u}\in {H^1_0\left((\Omega^{2r}_{L,\varepsilon})^c\right) }$ such that
(i)
(4.4)\begin{align}
\|\zeta_{u}\|\leq 2d.
\end{align}(ii)
(4.5)\begin{align}
J'_{\lambda,\varepsilon}(u)\zeta_{u}\geq a\left(\Gamma_{\varepsilon}(u)-\sigma_{\varepsilon}\right),
\end{align}where σɛ is independent of
$p, u$ and satisfies
$\lim\limits_{\varepsilon\rightarrow0}\sigma_{\varepsilon}=0$.(iii)
(4.6)\begin{align}
\Gamma'_{\varepsilon}(u)\zeta_{u}\geq a\left(\Gamma_{\varepsilon}(u)-\sigma_{\varepsilon}\right).
\end{align}where
$a, \sigma_{\varepsilon} \gt 0$ are as in (ii).
Proof. (i) Denote
\begin{align*}
\zeta_{u}(x)=\sum_{i=1}^{\ell}\varphi_{\varepsilon,i}(x)u(x)\in {H^1_0\left((\Omega^{2r}_{L,\varepsilon})^c\right) },
\end{align*} where 
$\varphi_{\varepsilon,i}(x)\in C_0^{\infty}\left(\mathbb{R}^N,[0,1]\right)$ is a cut-off function satisfying 
$\varphi_{\varepsilon,i}=1$ on 
$(\Omega_{\varepsilon,i}^{8r})^c$ and 
$\varphi_{\varepsilon,i}=0$ on 
$\Omega^{2r}_{\varepsilon,i}$ for each 
$i=1,\cdots,\ell$. Clearly, for ɛ > 0 small, it holds 
$\|\zeta_{u}\|_{{H^1_0\left((\Omega^{2r}_{L,\varepsilon})^c\right) }}\leq 2d$.
(ii) From Propositions 2.1, we calculate that
\begin{align*}
\int_{\mathbb{R}^N}\int_{(\Omega^{2r}_{L,\varepsilon})^c}
\frac{|u(x)|^p|u(y)|^{p}}{|x-y|^{N-\alpha}}=&\int_{\mathbb{R}^N}\int_{{(\Omega_{L,\varepsilon}^{8r})^c}}
\frac{|u(x)|^p|u(y)|^{p}}{|x-y|^{N-\alpha}}\\ &+
\int_{\mathbb{R}^N}\int_{\Omega_{L,\varepsilon}^{8r}\setminus\Omega^{2r}_{L,\varepsilon}}
\frac{|u(x)|^p|u(y)|^{p}}{|x-y|^{N-\alpha}},\\
\leq&\int_{{(\Omega_{L,\varepsilon}^{8r})^c}}\int_{{(\Omega_{L,\varepsilon}^{8r})^c}}
\frac{|u(x)|^p|u(y)|^{p}}{|x-y|^{N-\alpha}}\\ &+
2\int_{\mathbb{R}^N}\int_{\Omega_{L,\varepsilon}^{8r}\setminus\Omega^{2r}_{L,\varepsilon}}
\frac{|u(x)|^p|u(y)|^{p}}{|x-y|^{N-\alpha}}\\ &+\int_{\Omega^{2r}_{L,\varepsilon}}\int_{{(\Omega_{L,\varepsilon}^{8r})^c}}
\frac{|u(x)|^p|u(y)|^{p}}{|x-y|^{N-\alpha}}\\
\leq & C_{\varepsilon}+C\|u\|^{2p}_{H^1_0((\Omega_{L,\varepsilon}^{8r})^c)}.
\end{align*} where 
$\lim\limits_{\varepsilon\rightarrow 0}C_{\varepsilon}=0$. Combining with the above inequality, a direct calculation gives
\begin{align}
J'_{\lambda,\varepsilon}(u)\zeta_{u}=&\int_{\mathbb{R}^N}\nabla u\nabla \sum_{i=1}^{\ell}\varphi_{\varepsilon,i}(x)u(x)dx+\int_{\mathbb{R}^N}\lambda V_{\varepsilon}(x) u \sum_{i=1}^{\ell}\varphi_{\varepsilon,i}(x)u(x)dx \nonumber\\
&-\int_{\mathbb{R}^N}\int_{(\Omega^{2r}_{L,\varepsilon})^c}
\frac{|u(x)|^p|u(y)|^{p}\sum_{i=1}^{\ell}\varphi_{\varepsilon,i}(y)}{|x-y|^{N-\alpha}} \nonumber\\
\geq &\int_{{(\Omega_{L,\varepsilon}^{8r})^c}}|\nabla u|^2dx+\int_{\Omega_{L,\varepsilon}^{8r}\setminus\Omega_{L,\varepsilon}}u\nabla u\nabla \sum_{i=1}^{\ell}\varphi_{\varepsilon,i}(x)dx \nonumber\\
&-\int_{\mathbb{R}^N}\int_{(\Omega^{2r}_{L,\varepsilon})^c}
\frac{|u(x)|^p|u(y)|^{p}}{|x-y|^{N-\alpha}} \nonumber\\
\geq &\int_{{(\Omega_{L,\varepsilon}^{8r})^c}}|\nabla u|^2dx-C_{\varepsilon}
-D_{\varepsilon}-C\|u\|^{2p}_{H^1_0((\Omega_{L,\varepsilon}^{8r})^c)},
\end{align} where 
$\lim\limits_{\varepsilon\rightarrow 0}D_{\varepsilon}=0$. Controlling d small, there is a > 0 such that
\begin{align}
\int_{{(\Omega_{L,\varepsilon}^{8r})^c}}|\nabla u|^2dx-C\|u\|^{2p}_{H^1_0((\Omega_{L,\varepsilon}^{8r})^c)}\geq a \int_{{(\Omega_{L,\varepsilon}^{8r})^c}}|\nabla u|^2dx.
\end{align} Let 
$\sigma_{\varepsilon}=\left(C_{\varepsilon}+D_{\varepsilon}\right)/a$ and combining with (4.7)–(4.8), (4.5) is proved.
(iii) (4.6) can be shown as in (ii).
Let 
$u\in X _{\varepsilon}^{d}$ and 
$X(u)=\frac{\zeta_{u}}{\|\zeta_{u}\|}$. By Lemma 4.2, we have the following corollary.
Corollary 4.3. There is a locally Lipschitz vector field
\begin{align*}
X(u):X_{\varepsilon}^{d}\rightarrow {H^1_0\left((\Omega^{2r}_{L,\varepsilon})^c\right) }
\end{align*} such that, for all 
$u\in X_{\varepsilon}^{d}$,
(i)
$\|X(u)\|_{H^1}\leq 1$.(ii)
$J'_{\lambda,\varepsilon}(u)X(u),
\Gamma'_{\varepsilon}(u)X(u) \gt 0$ for
$\Gamma_{\varepsilon}(u) \gt \sigma_{\varepsilon}$. Here
$\sigma_{\varepsilon} \gt 0$ is as in Lemma 4.2.
We consider the following equation in H: for 
$u\in X_{\varepsilon}^{d}$,
\begin{align}
\begin{cases}
\frac{dW}{dT}=-\rho(\Gamma_{\varepsilon}(W))X(W),\\
W(0,u)=u,
\end{cases}
\end{align} where the function 
$\rho(r)\in C^{\infty}([0,\infty),[0,1])$ satisfies
\begin{align*}
\rho(r)=
\begin{cases}1,\ \text{for}\ r\geq2\sigma_{\varepsilon},\\
0,\ \text{for}\ r\leq\sigma_{\varepsilon}.
\end{cases}
\end{align*} Obviously, problem (4.9) has a unique solution, which we may assume to be 
$W(T,u)$. We next investigate some properties of 
$W(T,u)$.
Lemma 4.4. Let 
$W(T,u)$ be the unique solution of (4.9). Then we have the following conclusions:
(i)
$W(T,u)\in Y_{\varepsilon}^{d}$ for all T > 0 and
$u\in X_{\varepsilon}^{d}$.(ii) There holds
\begin{align*}
&\frac{d}{dT}\Gamma_{\varepsilon}(W(T,u))\leq 0,\ \text{for all}\ u\in X_{\varepsilon}^{d},\\
&\frac{d}{dT}\Gamma_{\varepsilon}(W(T,u)) \lt 0,\ \text{if}\ \Gamma_{\varepsilon}(W) \gt \sigma_{\varepsilon}.
\end{align*}(iii) There holds
\begin{align*}
\frac{d}{dT}J_{\lambda,\varepsilon}(\widetilde{W}(T,u)(x))\leq0,
\end{align*}where
\begin{align*}
\widetilde{W}(T,u)(x)=
\begin{cases}
u(x)\ &\text{for}\ x\in \Omega_{L,\varepsilon}^{2r},\\
W(T,u)\ &\text{for}\ x\in \left(\Omega_{L,\varepsilon}^{2r}\right)^c.
\end{cases}
\end{align*}(iv) There exists
$T_{\varepsilon} \gt 0$ such that for all
$u\in X_{\varepsilon}^{d}$,
\begin{align*}
\Gamma_{\varepsilon}(W(T_{\varepsilon},u))\leq \sigma_{\varepsilon}.
\end{align*}
Proof. First, we show (ii). Obviously, there holds 
$\frac{d}{dT}\Gamma_{\varepsilon}(W(T,u))\leq 0$ for all 
$u\in X_{\varepsilon}^{d}$. It follows from 
$\Gamma_{\varepsilon}(W) \gt \sigma_{\varepsilon}$ that 
$\rho(\Gamma_{\varepsilon}(W)) \gt 0$. Hence, combining with Lemma 4.2,
\begin{align*}
\frac{d}{dT}\Gamma_{\varepsilon}(W(T,u))=-\rho(\Gamma_{\varepsilon}(W))\Gamma'_{\varepsilon}(W)X(W) \lt 0.
\end{align*} By (ii), we get that 
$\Gamma_{\varepsilon}(W(T,u))$ is non-increasing. This suggests that 
$W(T,u)\in Y_{\varepsilon}^{d}$ for 
$T\in [0,\infty)$. Hence, (i) follows. (iv) is easily inferred from (ii). Since 
$W(T,u)|_{x\in \cup_{i=1}^{\ell}\partial\Omega_{\varepsilon,i}^{2r}}$ is independent of T, one has 
$\widetilde{W}(T,u)(x)\in H^1_0(\mathbb{R}^N)$. Hence, from Proposition 4.3 (ii), we can see that (iii) holds.
The proof of Proposition 4.1: For any 
$u\in X_{\varepsilon}^d$, we define 
$\tau_{\varepsilon}=\widetilde{W}(T_{\varepsilon},u)$. It follows from Lemma 4.4 that τɛ satisfies (4.1)–(4.3) mentioned in Proposition 4.1.
5. Existence of
$\ell$-bump critical points
In this section, we complete the proof of Theorems 1.1 and 1.2.
Proposition 5.1. There is 
$\varepsilon_0 \gt 0$ such that for 
$\varepsilon\in (0, \varepsilon_0)$ and 
$\lambda \gt \overline{\lambda}_{\varepsilon}$, 
$J_{\lambda,\varepsilon}|_{M_{\varepsilon}}$ has a critical point 
$u_{\lambda, \varepsilon}\in X_{\varepsilon}^d$ with an associated Lagrange multiplier 
$\mu_{\lambda,\varepsilon}$ satisfying
\begin{align*}
&\lim\limits_{\varepsilon\rightarrow0}J_{\lambda,\varepsilon}(u_{\lambda,\varepsilon})= E(s^0),\\
&\lim\limits_{\varepsilon\rightarrow0}\sup\limits_{\lambda \gt \overline{\lambda}_{\varepsilon}}
|\mu_{\lambda,\varepsilon}-\mu_{\ell^{-1}}|=0.
\end{align*}Proof of proposition 5.1
Choosing d > 0 smaller if necessary, as in Proposition 3.3 we can see for some 
$\nu_{d}, k_{\ell} \gt 0$ independent of ɛ that for ɛ > 0 small,
\begin{align}
\|J'_{\lambda,\varepsilon}|_{M_{\varepsilon}}(u)\|_{T^*_{u}( M_{\varepsilon})}\geq \nu_{d}\quad \text{for}\ u\in \left(X_{\varepsilon}^d\setminus X_{\varepsilon}^{d/4}\right)\cap\left[ J_{\lambda,\varepsilon}\leq E(s^0)+k_{\ell}\right].
\end{align} If 
$J_{\lambda,\varepsilon}(u)$ does not have any critical point with 
$J_{\lambda,\varepsilon}\leq E(s^0)+k_{\ell}$ in 
$X_{\varepsilon}^d$, then from Proposition 3.2, there exists 
$\kappa_{\lambda,\varepsilon} \gt 0$ such that
\begin{align}
\|J'_{\lambda,\varepsilon}\|_{M_{\varepsilon}}(u)\|_{T^*_u( M_{\varepsilon})}\geq \kappa_{\lambda,\varepsilon}\quad \text{for}\ u\in X_{\varepsilon}^d\cap\left[J_{\lambda,\varepsilon}\leq E(s^0)+k_{\ell}\right].
\end{align}For convenience, we set
\begin{align*}
&\left[J_{\lambda,\varepsilon}\leq b_1\right]:=\left\{u\in E_{\varepsilon}| J_{\lambda,\varepsilon}(u)\leq b_1\right\},\\
&\left[b_2\leq J_{\lambda,\varepsilon}\right]:=\left\{u\in E_{\varepsilon}| J_{\lambda,\varepsilon}(u)\geq b_2\right\},\\
&\left[b_2\leq J_{\lambda,\varepsilon}\leq b_1\right]:=\left\{u\in E_{\varepsilon}| b_2\leq J_{\lambda,\varepsilon}(u)\leq b_1\right\}
\end{align*} for any 
$b_1, b_2\in\mathbb{R}$.
If 
$J_{\lambda,\varepsilon}(u)$ does not have any critical point with 
$J_{\lambda,\varepsilon}\leq E(s^0)+k_{\ell}$ in 
$X_{\varepsilon}^d$, by (5.1)–(5.2), we can find a deformation flow 
$\eta: M_{\varepsilon}\rightarrow M_{\varepsilon}$ be defined as the following lemma.
Lemma 5.2. Under the assumption (5.1) and (5.2), for any 
$k\in \left(0, \min\left\{k_{\ell},\frac{d\nu_{d}}{16}\right\}\right)$, there exists an 
$\varepsilon_{k} \gt 0$ such that for any 
$\varepsilon\in(0,\varepsilon_k)$ and 
$\lambda\in (\overline{\lambda}_{\varepsilon},+\infty)$, there is a deformation 
$\eta: M_{\varepsilon}\rightarrow M_{\varepsilon}$ satisfying
(1)
$\eta(u)=u$ if
$u\in M_{\varepsilon}\setminus\left(X_{\varepsilon}^d\cap\left[E(s^0)-2k_{\ell}\leq J_{\lambda,\varepsilon}\leq E(s^0)+k_{\ell}\right]\right)$.(2)
$J_{\lambda,\varepsilon}(\eta(u))\leq J_{\lambda,\varepsilon}(u)$ if
$u\in M_{\varepsilon}$.(3)
$\eta(u)\in X_{\varepsilon}^d\cap\left[ J_{\lambda,\varepsilon}\leq E(s^0)-k\right]$ if
$u\in X_{\varepsilon}^{d/4}\cap\left[J_{\lambda,\varepsilon}\leq E(s^0)+k\right]$.
Proof. Set 
$\nu_{\varepsilon}: M_{\varepsilon}\rightarrow E_{\varepsilon}$ is locally Lipschitz continuous tangent vector field satisfying 
$\nu_{\varepsilon}(u)\in T_{u} M_{\varepsilon}$ for every 
$u\in M_{\varepsilon}$ and
\begin{align}
&\|\nu_{\varepsilon}\|\leq1\ \text{and}\ \langle J'_{\lambda,\varepsilon}|_{M_{\varepsilon}}(u),\nu_{\varepsilon}(u)\rangle\geq 0\ \text{for}\ u\in M_{\varepsilon}, \nonumber\\
&\langle J'_{\lambda,\varepsilon}|_{M_{\varepsilon}}(u),\nu_{\varepsilon}(u)\rangle\geq
\frac{\kappa_{\lambda,\varepsilon}}{2}\ \text{for}\ u\in X_{\varepsilon}^d\cap \left[J_{\lambda,\varepsilon}
\leq E(s^0)+k_{\ell}\right], \nonumber\\
&\langle J'_{\lambda,\varepsilon}|_{M_{\varepsilon}}(u),\nu_{\varepsilon}(u)\rangle\geq \frac{\nu_{d}}{2}\ \text{for}
\ u\in (X_{\varepsilon}^d\setminus X_{\varepsilon}^{d/4})\cap\left[J_{\lambda,\varepsilon}\leq E(s^0)+k_{\ell}\right].
\end{align} Then for every 
$u\in M_{\varepsilon}$, the flow 
$\overline{\eta}(t,u)$ is defined as a solution of
\begin{align*}
\begin{cases}
\frac{\rm{d}\overline{\eta}}{\rm{d}t}(t,u)=-\zeta_1(\overline{\eta}(t,u))\zeta_2(\overline{\eta}(t,u))
\nu_{\varepsilon}(\overline{\eta}(t,u)),\\
\overline{\eta}(0,u)=u,
\end{cases}
\end{align*} where locally Lipschitz functions 
$\zeta_1, \zeta_2: M_{\varepsilon}\rightarrow [0,1]$ are defined as
\begin{align*}
\zeta_1(u)=\begin{cases} 1,&\quad\ u\in X_{\varepsilon}^{d/2},\\ 0,&\quad\ u\notin X_{\varepsilon}^d,\end{cases}\quad
\zeta_2(u)=\begin{cases} 1,&\quad\ u\in \left[E(s^0)-k_{\ell}\leq J_{\lambda,\varepsilon}\right],\\ 0,&\quad\ u\notin \left[E(s^0)-2k_{\ell}\leq J_{\lambda,\varepsilon}\right].\end{cases}
\end{align*} It is easy to obtain that 
$\overline{\eta}\in C([0,+\infty)\times M_{\varepsilon}, M_{\varepsilon})$ satisfies
\begin{align*}
\begin{cases}
\overline{\eta}(\cdot,u)\in X_{\varepsilon}^d\quad \text{if}\ u\in X_{\varepsilon}^d,\\
\overline{\eta}(\cdot,u)\equiv u\quad \text{if}\ u\notin \left[E(s^0)-2k_{\ell}\leq J_{\lambda,\varepsilon}\right]\cap X_{\varepsilon}^d,
\end{cases}
\end{align*} For 
$u\in M_{\varepsilon}$ and 
$t\geq0$, we have
\begin{align}
\frac{\rm{d}}{\rm{d}t} J_{\lambda,\varepsilon}(\overline{\eta}(t,u))=-\zeta_1(\overline{\eta}(t,u))\zeta_2(\overline{\eta}(t,u))\langle J'_{\lambda,\varepsilon}|_{M_{\varepsilon}}(\overline{\eta}(t,u)),\nu_{\varepsilon}(t,u))\rangle\leq0.
\end{align} We set 
$\eta:=\overline{\eta}\left(\frac{4k}{\kappa_{\lambda,\varepsilon}},u\right)$ and note that it satisfies the desired first two conclusions. To show (3), we argue by contradiction and assume 
$u_{\varepsilon}\in X_{\varepsilon}^{d/4}\cap\left[J_{\lambda,\varepsilon}\leq E(s^0)+k\right]$ exists such that
\begin{align}
J_{\lambda,\varepsilon}(\overline{\eta}(t,u))\in\left( E(s^0)-k,E(s^0)+k\right]\quad \text{for all}\ t\in\left[0,\frac{4k}{\kappa_{\lambda,\varepsilon}}\right].
\end{align} If 
$\overline{\eta}\in X_{\varepsilon}^{d/2}$ for every 
$t\in \left[0,\frac{4k}{\kappa_{\lambda,\varepsilon}}\right]$, it follows from (5.3) and (5.4) that
\begin{align*}
J_{\lambda,\varepsilon}(\eta(u_{\varepsilon}))=&J_{\lambda,\varepsilon}
\left(\overline{\eta}\left(\frac{4k}{\kappa_{\lambda,\varepsilon}},u_{\varepsilon}\right)\right)\\
=&\int^{\frac{4k}{\kappa_{\lambda,\varepsilon}}}_0\frac{\rm{d}}{\rm{d}t}J_{\lambda,\varepsilon}
(\overline{\eta}(t,u_{\varepsilon}))dt+J_{\lambda,\varepsilon}(u_{\varepsilon})
\\ \leq&-\int_{0}^{\frac{4k}{\kappa_{\lambda,\varepsilon}}}
\frac{\kappa_{\lambda,\varepsilon}}{2}dt+E(s^0)+k\\
=&E(s^0)-k.
\end{align*} This contradicts (5.5). If 
$\overline{\eta}(t,u_{\varepsilon})\notin X_{\varepsilon}^{d/2}$ for some 
$t\in\left[0,\frac{4k}{\kappa_{\lambda,\varepsilon}}\right]$, then we can find an interval 
$[t_1, t_2]\subset\left[0,\frac{4k}{\kappa_{\lambda,\varepsilon}}\right]$ such that 
$\overline{\eta}(t_1,\mu_{\varepsilon})\in\partial X_{\varepsilon}^{d/4}$, 
$\overline{\eta}(t_2,u_{\varepsilon})\in\partial X_{\varepsilon}^{d/2}$ and 
$\overline{\eta}(t,u_{\varepsilon})\subset X_{\varepsilon}^{d/2}\setminus X_{\varepsilon}^{d/4}$ for all 
$t\in (t_1, t_2)$. Since
\begin{align*}
t_2-t_1\geq\int_{t_1}^{t_2}\|\nu_{\varepsilon}(\overline{\eta}(t,u_{\varepsilon}))\|dt\geq \left\|\int_{t_1}^{t_2}\frac{\rm{d}\overline{\eta}}{\rm{d}t}(t,u_{\varepsilon})dt\right\|=
\overline{\eta}(t_2,u_{\varepsilon})-\overline{\eta}(t_1,u_{\varepsilon})\geq d/4,
\end{align*}according to (5.3) and (5.4), there holds
\begin{align*}
J_{\lambda,\varepsilon}(\eta(u_{\varepsilon}))=J_{\lambda,\varepsilon}\left(\overline{\eta}
\left(\frac{4k}{\kappa_{\lambda,\varepsilon}},u_{\varepsilon}\right)\right)
\leq& J_{\lambda,\varepsilon}(\overline{\eta}(t_2,u_{\varepsilon}))\leq\int_{t_1}^{t_2}\frac{\rm{d}}{\rm{d}t}
J_{\lambda,\varepsilon}(\overline{\eta}(t,u_{\varepsilon}))\\ &+
J_{\lambda,\varepsilon}(\overline{\eta}(t_1,u_{\varepsilon}))\\
\leq& -\frac{\nu_{d}}{2}(t_2-t_1)+E(s^0)+k\\
\leq& -2k+E(s^0)+k\\ \leq&E(s^0)-k,
\end{align*}which also contradicts (5.5). Therefore (c) holds.
As shown by (2.3) and (2.4), the geometric structure of the energy functional I(u) depends on p. Moreover, from Lemma 2.2, we also see that the magnitude relation between E(s) and 
$E(s^0)$ depends on p. Thus, we study it in two parts.
Step 1: L 2-subcritical case 
$p\in \left(\frac{N+\alpha}{N}, \frac{N+\alpha+2}{N}\right)$
We define the 
$(\ell-1)$-dimensional initial path by
\begin{align*}
\widetilde{\chi}_{\varepsilon}(s)=C_{\varepsilon}\sum_{i=1}^{\ell}(\widetilde{\varphi}_{\varepsilon} w_{s_i})(\cdot-q_{i,\varepsilon})\quad \text{for}\ s=(s_1, \cdots,s_{\ell})\in \mathcal{S}_{\delta},
\end{align*} where 
$\mathcal{S}_{\delta}:=\{s\in \mathcal{S}^{\ell}_+|\ |s-s^0|\leq\delta, \delta \gt 0\}$, 
$\theta:=\frac{1}{4}\min\limits_{1\leq i\leq\ell}dist (q_i, \partial \Omega_i)$, 
$w_{s_i}\in S_{s_i}$, 
$\widetilde{\varphi}_{\varepsilon}(\cdot)=\widetilde{\varphi}(\varepsilon\cdot)$ is defined for some 
$\widetilde{\varphi}\in C^1_0(B_{2\theta}(0),[0,1])$ satisfying 
$\widetilde{\varphi}=1$ in 
$B_{\theta}(0)$ and 
$|\nabla \widetilde{\varphi}|\leq 2/\theta$, and
\begin{align*}
C_{\varepsilon}:=\left|\sum^{\ell}_{i=1}\left(\widetilde{\varphi}_{\varepsilon}w_{s_i}(\cdot-q_{i,\varepsilon})\right)\right|^{-1}_{2}.
\end{align*}Fix δ > 0 small enough (independent of ɛ small) such that
\begin{align*}
\widetilde{\chi}_{\varepsilon}(s)\in X_{\varepsilon}^{d/4}\quad \text{for}\ s\in \mathcal{S}_{\delta}
\end{align*}and
\begin{align*}
C_{\varepsilon}\rightarrow 1\quad \text{as}\ \varepsilon\rightarrow 0\ \text{uniformly in}\ \mathcal{S}_{\delta}.
\end{align*} According to 
$V_{\varepsilon}=0$ in 
$supp\ \widetilde{\varphi}_{\varepsilon}(\cdot-q_{i,\varepsilon})\subset\Omega_{i,\varepsilon}$ for each 
$i=1,\cdots,\ell$, we obtain that 
$J_{\lambda,\varepsilon}(\widetilde{\chi}_{\varepsilon}(s))$ is independent of λ and
\begin{align*}
J_{\lambda,\varepsilon}(\widetilde{\chi}_{\varepsilon}(s))\rightarrow E(s)\ \text{as}\ \varepsilon\rightarrow 0\ \text{uniformly for}\ s\in \mathcal{S}_{\delta}
\end{align*} Moreover, using 
$E(s) \lt E(s^0)$ introduced in Lemma 2.2, we have
\begin{align}
\lim\limits_{\varepsilon\rightarrow 0}\max_{s\in \mathcal{S}_{\delta}} J_{\lambda,\varepsilon}(\widetilde{\chi}_{\varepsilon}(s))=E(s^0).
\end{align} Fix 
$k\in(0,\min\{k_{\ell}, ,\frac{d\nu_{d}}{16}\})$ such that
\begin{align}
\lim\limits_{\varepsilon\rightarrow 0}\max\limits_{s\in S_{\partial \mathcal{S}_{\delta}}} J_{\lambda,\varepsilon}(\widetilde{\chi}_{\varepsilon}(s))=\max\limits_{s\in S_{\partial \mathcal{S}_{\delta}}} E(s) \lt E(s^0)-2k.
\end{align}Next, we will have the contradiction in the L 2-subcritical case. By (5.6) and Lemma 5.2 (3), for small ɛ > 0, we have
\begin{align}
J_{\lambda,\varepsilon}(\eta(\widetilde{\chi}_{\varepsilon}(s)))\leq E(s^0)-k \quad \text{and}\ \eta(\widetilde{\chi}_{\varepsilon}(s))\in X_{\varepsilon}^d\ \text{for}
\ s\in \mathcal{S}_{\delta}.
\end{align} When 
$\ell\geq2$, it follows from (5.8) and lemma 5.2 (1) that for smaller ɛ > 0, there holds
\begin{align*}
\eta(\widetilde{\chi}_{\varepsilon}(s))=\widetilde{\chi}_{\varepsilon}(s)\quad \text{if}\ s\in\partial \mathcal{S}_{\delta}.
\end{align*}Set
\begin{align*}
\Psi_{\varepsilon,i}(u)=\left(\int_{\Omega_{i,\varepsilon}^{8r}}|u|^2\right)\left(\int_{\Omega_{L,\varepsilon}^{8r}}|u|^2\right)^{-1}\quad \forall\ u\in M_{\varepsilon},\ i=1,\cdots,\ell.
\end{align*} To get a contradiction, we now claim that there exists 
$s^1\in \mathcal{S}_{\delta}$ such that
\begin{align}
\Psi_{\varepsilon,i}(\tau_{\varepsilon}(\eta(\widetilde{\chi}_{\varepsilon}(s^1))))=\ell^{-1}\quad \text{for each}\ i=1,\cdots,\ell.
\end{align}By Proposition 4.1, we have
\begin{align}
\begin{cases}
J_{\lambda,\varepsilon}\circ\tau_{\varepsilon}\circ\eta\circ\widetilde{\chi}_{\varepsilon}(s)\leq E(s^0)-k,\ \text{if}\ s\in \mathcal{S}_{\delta}, \\
\|\nabla \tau_{\varepsilon}(\eta(\widetilde{\chi}_{\varepsilon}(s)))\|_{L^2((\Omega^{8r}_{L,\varepsilon})^c)}=o_{\varepsilon}(1),\\
\tau_{\varepsilon}\circ\eta\circ\widetilde{\chi}_{\varepsilon}(s)(x)=\eta\circ\widetilde{\chi}_{\varepsilon}(s)(x),\ \text{if}\ x\in \Omega_{L,\varepsilon}^{2r}, \\
\tau_{\varepsilon}\circ\eta\circ\widetilde{\chi}_{\varepsilon}(s)(x)=\widetilde{\chi}_{\varepsilon}(s)(x),\ \text{if}\ s\in\partial \mathcal{S}_{\delta}\ \text{and}\ x\in \Omega_{L,\varepsilon}^{2r}.
\end{cases}
\end{align} We use the Brouwer degree to prove (5.9) holds. We consider the map 
$f_{\varepsilon}:\mathcal{S}_{\delta}\rightarrow \mathbb{R}^{\ell}$ defined by
\begin{align*}
f_{\varepsilon}(s)=(\Psi_{\varepsilon,1}\circ\widetilde{\chi}_{\varepsilon}(s),\cdots,\Psi_{\varepsilon,\ell}\circ\widetilde{\chi}_{\varepsilon}(s)).
\end{align*}Clearly,
\begin{align*}
\Psi_{\varepsilon,i}\circ\widetilde{\chi}_{\varepsilon}(s)\rightarrow |w_{s_i}|^2_2=s_i\quad \text{as}\ \varepsilon\rightarrow 0\ \text{uniformly for every}\ s\in \mathcal{S}_{\delta}.
\end{align*} Hence 
$s^0\not\in f_{\varepsilon}(\partial \mathcal{S}_{\delta})$ for ɛ small and the Brouwer degree
\begin{align*}
\deg(f_{\varepsilon}, \mathcal{S}_{\delta}, s^0)=\deg(id, \mathcal{S}_{\delta}, s^0)=1\quad \text{for}\ \varepsilon\ \text{small}.
\end{align*}Then we set
\begin{align*}
h_{\varepsilon}(s):=(\Psi_{\varepsilon,1}\circ\tau_{\varepsilon}
\circ\eta\circ\widetilde{\chi}_{\varepsilon}(s),\cdots,\Psi_{\varepsilon,\ell}\circ\tau_{\varepsilon}\circ\eta\circ\widetilde{\chi}_{\varepsilon}(s)).
\end{align*} Combining with (5.10), we deduce 
$s^1\not\in g_{\varepsilon}(\partial \mathcal{S}_{\delta})$ for ɛ small and the Brower degree
\begin{align*}
\deg(h_{\varepsilon},\mathcal{S}_{\delta},s^0)=\deg(f_{\varepsilon}, \mathcal{S}_{\delta}, s^0)=1\quad \text{for}\ \varepsilon\ \text{small},
\end{align*}and the claim is verified. Now we denote
\begin{align*}
U_{\varepsilon}:=
\begin{cases}
\tau_{\varepsilon}(\eta(\widetilde{\chi}_{\varepsilon}(1)))\quad \text{if}\ \ell=1,\\
\tau_{\varepsilon}(\eta(\widetilde{\chi}_{\varepsilon}(s^1)))\quad \text{if}\ \ell\geq 2,
\end{cases}
U_{i,\varepsilon}:=\xi_{i,\varepsilon}U_{\varepsilon},
\end{align*} where 
$\xi_{i,\varepsilon}\in C^{\infty}_0(\Omega^{10r}_{i,\varepsilon},[0,1])$ is a cut-off function such that 
$\xi_{i,\varepsilon}=1$ on 
$\Omega_{i,\varepsilon}^{8r}$ and 
$|\nabla \xi_{i,\varepsilon}|\leq\varepsilon/r$, for each 
$i=1,2,\cdots,\ell$. It follows from 
$\|\nabla U_{\varepsilon}\|_{L^2(\Omega_{L,\varepsilon}^{8r})^c}=o_{\varepsilon}(1)$ and the fact that 
$\Omega_{i,\varepsilon}^{8r}\cap\Omega_{j,\varepsilon}^{10r}=\emptyset$ for each 
$i\not\equiv j$, we conclude that
\begin{align}
|U_{i,\varepsilon}|^2_2=\Psi_{\varepsilon,i}(U_{\varepsilon})+o_{\varepsilon}(1)=\ell^{-1}
+o_{\varepsilon}(1).
\end{align}Thus
\begin{align}
I_{\ell^{-1}}\leq I\left(\ell^{-1/2}U_{i,\varepsilon}|U_{i,\varepsilon}|^{-1}_2\right)=I(U_{i,\varepsilon})+o_{\varepsilon}(1).
\end{align}Using (5.8), (5.11) and (5.12), we calculate that
\begin{align*}
E(s^0)-k\geq& J_{\lambda,\varepsilon}(U_{\varepsilon})\geq I(U_{\varepsilon})\\
\geq&\sum_{i=1}^{\ell}\frac{1}{2}\int_{\Omega_{i,\varepsilon}^{8r}}|\nabla U_{\varepsilon}|^2-\frac{1}{2p}
\int_{\Omega_{i,\varepsilon}^{8r}}\int_{\Omega_{i,\varepsilon}^{8r}}\frac{|U_{\varepsilon}|^p|U_{\varepsilon}|^p}{|x-y|^{N-\alpha}}+o_{\varepsilon}(1)\\
=&\sum_{i=1}^{\ell}I(U_{i,\varepsilon})+o_{\varepsilon}(1)\\
\geq&\ell I_{\ell^{-1}}+o_{\varepsilon}(1)\\
=&E(s^0)+o_{\varepsilon}(1),
\end{align*}which is a contradiction.
Claim 2: L 2-supercritical case 
$p\in\left(\frac{N+\alpha+2}{N},\frac{N+\alpha}{N-2}\right)$
Similar to Claim 1, according to the geometry of the energy functional I(u), we set the initial path by
\begin{align*}
\vartheta_{\varepsilon}(t):=\widetilde{C_{\varepsilon}}\sum\limits_{i=1}^{\ell}
\left(t_i^{N/2}\widetilde{\varphi}_{\varepsilon}(t_i(\cdot - q_{i,\varepsilon}))u_i\left(t_i(\cdot -q_{i,\varepsilon})\right)\right)\quad \text{for}\ t=(t_1,\cdots,t_{\ell})\in\mathbb{R}^{\ell}_{+},
\end{align*} where 
$u_i\in S_{\ell^{-1}}$, 
$q_i\in \Omega_i$ and 
$q_{i,\varepsilon}=\frac{1}{\varepsilon}q_i$ for 
$i=1,2,\cdots, \ell$, 
$\widetilde{\varphi}_{\varepsilon}(\cdot)=\widetilde{\varphi}(\varepsilon\cdot)$ with 
$\widetilde{\varphi}\in C^1_0(B_{2\theta}(0),[0,1])$ satisfying 
$\widetilde{\varphi}=1$ in 
$B_{\theta}(0)$ and 
$|\nabla\widetilde{\varphi}|\leq 2/\theta$,
\begin{align*}
\widetilde{C_{\varepsilon}}:=\left|\sum\limits_{i=1}^{\ell}
\left(t_i^{N/2}\widetilde{\varphi}_{\varepsilon}(t_i(\cdot - q_{i,\varepsilon}))u_i\left(t_i(\cdot -q_{i,\varepsilon})\right)\right))\right|^{-1}_2.
\end{align*} Then we fix 
$\sigma^{\ast} \gt 0$ small enough (independent of ɛ small) such that
\begin{align*}
\vartheta_{\varepsilon}(t)\in X_{\varepsilon}^{d/4}\quad \text{for}\ t\in[1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell}.
\end{align*} Obviously, there holds 
$\lim\limits_{\varepsilon\rightarrow0}\widetilde{C_{\varepsilon}}=1$ and
\begin{align*}
\lim\limits_{\varepsilon\rightarrow0}J_{\lambda,\varepsilon}(\vartheta_{\varepsilon}(t))&=\sum\limits_{i=1}^{\ell} I\left(t_i^{N/2}u_i(t_i\cdot)\right)\\&=\sum\limits_{i=1}^{\ell}\left(\frac{1}{2}t_i^2|\nabla u_i|_2^2-\frac{1}{2p}t_i^{\alpha-1}\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_i(x)|^p|u_i(y)|^p}{|x-y|^{N-\alpha}}\right)
\end{align*} for any 
$t\in[1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell}$. Since 
$I(u_i) \gt I\left(t_i^{N/2}u_i(t_i\cdot)\right)$ for 
$t_i\in [1-\sigma^{\ast},1+\sigma^{\ast}]\setminus\{1\}$, we can fix 
$k\in (0,\min\{k_{\ell},\frac{d\nu_{d}}{16}\})$ such that
\begin{align}
\lim\limits_{\varepsilon\rightarrow0}\max\limits_{\partial[1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell}}J_{\lambda,\varepsilon}\circ\vartheta_{\varepsilon} \lt E(s^0)-2k \lt \lim\limits_{\varepsilon\rightarrow0}\max\limits_{[1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell}}
J_{\lambda,\varepsilon}\circ\vartheta_{\varepsilon}=E(s^0).
\end{align}It follows from (5.13) and Lemma 5.2 that
\begin{align}
J_{\lambda,\varepsilon}(\eta(\vartheta_{\varepsilon}(t)))\leq E(s^0)-k,\ \eta(\vartheta_{\varepsilon}(t))\in X_{\varepsilon}^d\ \text{for}\ t\in
[1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell}
\end{align}and
\begin{align*}
\eta(\vartheta_{\varepsilon}(t))=\vartheta_{\varepsilon}(t)\quad \text{if}\ t\in \partial [1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell}.
\end{align*} For each 
$i=1,\cdots,\ell$, define
\begin{align*}
\overline{\Psi}_{\varepsilon,i}(u)=\left(\int_{\Omega^{8r}_{i,\varepsilon}}|\nabla u|^2\right)^{\frac{1}{-Np+2+N+\alpha}}\left(\frac{1}{2p}
\int_{\Omega^{8r}_{i,\varepsilon}}\int_{\Omega^{8r}_{i,\varepsilon}}\frac{|\nabla u|^p|\nabla u|^p}{|x-y|^{N-\alpha}}\right)^{\frac{-1}{-Np+2+N+\alpha}}
\end{align*} Now we claim that there exists 
$t^1\in [1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell}$ such that
\begin{align}
\overline{\Psi}_{\varepsilon,i}(\chi(\eta(\vartheta_{\varepsilon}(t^1))))=1\quad \text{for each}\ i=1,\cdots,\ell.
\end{align}From Proposition 4.1, we obtain that
\begin{align}
\begin{cases}
J_{\lambda,\varepsilon}\circ\tau_{\varepsilon}\circ\eta\circ\vartheta_{\varepsilon}(t)\leq E(s^0)-k,\ \text{if}\ t\in [1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell},\\
\|\nabla \tau_{\varepsilon}(\eta(\vartheta_{\varepsilon}(t)))\|_{L^2((\Omega^{8r}_{L,\varepsilon})^c)}=o_{\varepsilon}(1),\\
\tau_{\varepsilon}\circ\eta\circ\vartheta_{\varepsilon}(t)(x)=
\eta\circ\vartheta_{\varepsilon}(t)(x),\ \text{if}\ x\in \Omega_{L,\varepsilon}^{2r},\\
\tau_{\varepsilon}\circ\eta\circ\vartheta_{\varepsilon}(t)(x)
=\vartheta_{\varepsilon}(t)(x),\ \text{if}\ t\in \partial [1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell}\ \text{and}\ x\in \Omega_{L,\varepsilon}^{2r}.
\end{cases}
\end{align} Define the map 
$\widetilde{f_{\varepsilon}}:[1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell}\rightarrow \mathbb{R}^{\ell}$:
\begin{align*}
\widetilde{f_{\varepsilon}}(t)=\left(\overline{\Psi}_{\varepsilon,1}\circ\vartheta_{\varepsilon}(t)
,\cdots,\overline{\Psi}_{\varepsilon,\ell}\circ\vartheta_{\varepsilon}(t)\right).
\end{align*} Since 
$supp \widetilde{\varphi}_{\varepsilon}(\cdot-q_{i,\varepsilon})\subset\Omega_{i,\varepsilon}\subset\Omega_{i,\varepsilon}^{8r}$ and 
$\Omega_{i,\varepsilon}^{8r}\cap\Omega_{j,\varepsilon}^{8r}=\emptyset$ for each i ≠ j, we deduce
\begin{align*}
\lim\limits_{\varepsilon\rightarrow0}\overline{\Psi}_{\varepsilon,i}\circ\vartheta_{\varepsilon}(t)=t_i
\end{align*} uniformly for 
$t\in [1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell}$. Then 
$t^0:=(1,\cdots,1)\in \mathbb{R}^{\ell}\setminus \widetilde{f_{\varepsilon}}(\partial [1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell})$ for ɛ small and the Brower degree
\begin{align*}
\deg(\widetilde{f_{\varepsilon}}, [1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell},t^0)=\deg(id, [1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell},t^0)=1
\quad \text{for}\ \varepsilon\ \text{small}.
\end{align*}Then we denote
\begin{align*}
\widetilde{h}_{\varepsilon}:=\left(\overline{\Psi}_{\varepsilon,1}\circ\tau_{\varepsilon}
\circ\eta\circ\vartheta_{\varepsilon}(t),\cdots,\overline{\Psi}_{\varepsilon,\ell}
\circ\tau_{\varepsilon}\circ\eta\circ\vartheta_{\varepsilon}(t)\right).
\end{align*} It follows from (5.16) that 
$t^0:=(1,\cdots,1)\in \mathbb{R}^{\ell}\setminus \widetilde{f_{\varepsilon}}(\partial [1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell})$ for ɛ small and the Brower degree
\begin{align*}
\deg(\widetilde{h}_{\varepsilon}, [1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell},t^0)=\deg(\widetilde{f_{\varepsilon}}, [1-\sigma^{\ast},1+\sigma^{\ast}]^{\ell},t^0)=1\quad \text{for}\ \varepsilon\ \text{small}.
\end{align*}Hence, the claim is proved. Now denote
\begin{align*}
\widetilde{U}_{\varepsilon}:=\tau_{\varepsilon}\left(\eta(\vartheta_{\varepsilon}(t^1))\right)
,\quad\widetilde{U}_{i,\varepsilon}
:=\xi_{i,\varepsilon}\widetilde{U}_{\varepsilon}\left(\sum\limits_{i=1}^{\ell}
|\xi_{i,\varepsilon}\widetilde{U}_{\varepsilon}|^2_2\right)^{-\frac{1}{2}}
\end{align*} where 
$\xi_{i,\varepsilon}\in C^{\infty}_{0}(\Omega_{i,\varepsilon}^{10r}[0,1])$ is a cut-off function such that 
$\xi_{i,\varepsilon}=1$ on 
$\Omega_{i,\varepsilon}^{8r}$ and 
$|\nabla \xi_{i,\varepsilon}|\leq\varepsilon/r$, for each 
$i=1,\cdots,\ell$. Since 
$\|\nabla \widetilde{U}_{\varepsilon}\|_{L^2((\Omega^{8r}_{L,\varepsilon})^c)}=o_{\varepsilon}(1)$, we deduce that
\begin{align}
\sum\limits_{i=1}^{\ell}|\xi_{i,\varepsilon}\widetilde{U}_{\varepsilon}|^2_2=1+o_{\varepsilon}(1).
\end{align}Combining with (5.15) and (5.17), we obtain that
\begin{align}
t_{i,\varepsilon}:=\overline{\Psi}_{\varepsilon,i}(\widetilde{U}_{\varepsilon})+o_{\varepsilon}(1)=1+o_{\varepsilon}(1).
\end{align}A direct calculation shows that
\begin{align*}
s^*:=\left(\left|t_{1,\varepsilon}^{-\frac{N}{2}}\widetilde{W}_{1,\lambda,\varepsilon}(t_{1,\varepsilon}^{-1}\cdot)\right|^2_2,\cdots,
\left|t_{\ell,\varepsilon}^{-\frac{N}{2}}\widetilde{W}_{\ell,\lambda,\varepsilon}
(t_{\ell,\varepsilon}^{-1}\cdot)\right|^2_2\right)\in \mathcal{S}_+^{\ell}
\end{align*}Hence, using Lemma 2.2, equation (5.17) and (5.18), we deduce that
\begin{align}
E(s^0)\leq E(s^*)\leq \sum\limits_{i=1}^{\ell} I\left(t_{i,\varepsilon}^{-\frac{N}{2}}\widetilde{U}_{i,\varepsilon}(t_{i,\varepsilon}^{-1}\cdot)\right)
=\sum\limits_{i=1}^{\ell}I(\widetilde{U}_{i,\varepsilon}
)+o_{\varepsilon}(1).
\end{align}However, according to (5.13), (5.17) and (5.19) we conclude that
\begin{align*}
E(s^0)-k\geq J_{\lambda,\varepsilon}(\widetilde{U}_{\varepsilon})\geq I(\widetilde{U}_{\varepsilon})
\geq&\sum_{i=1}^{\ell}\frac{1}{2}\int_{\Omega_{i,\varepsilon}^{8r}}|\nabla \widetilde{U}_{\varepsilon}|^2-\frac{1}{2p}
\int_{\Omega_{i,\varepsilon}^{8r}}\int_{\Omega_{i,\varepsilon}^{8r}}
\frac{|\widetilde{U}_{\varepsilon}|^p|\widetilde{U}_{\varepsilon}|^p}{|x-y|^{N-\alpha}}\\ &+o_{\varepsilon}(1)\\
=&\sum\limits_{i=1}^{\ell}I(\widetilde{U}_{i,\varepsilon}
)+o_{\varepsilon}(1)\\ \geq& E(s^0)+o_{\varepsilon}(1),
\end{align*}which is a contradiction.
This completes the proof of Proposition 5.1.
Proof of Theorems 1.1 and 1.2
Following the proof of Proposition 5.1, for ɛ sufficiently small and λ large enough, 
$J_{\lambda,\varepsilon}$ has a critical point 
$u_{\lambda,\varepsilon}$ with the associated Lagrange multipler 
$\mu_{\lambda,\varepsilon}$ in 
$X_{\varepsilon}^d$ satisfying 
$\lim\limits_{\varepsilon\rightarrow0}J_{\lambda,\varepsilon}(u_{\lambda,\varepsilon})=E(s^0)$. According to the strong maximum principle, we obtain that 
$u_{\lambda,\varepsilon} \gt 0$. Then we complete the proof of Theorem 1.1 by setting
\begin{align*}
\overline{u}_{\lambda,\varepsilon}=\omega^{\frac{2+\alpha}{2(-Np+N+\alpha+2)}}u_{\lambda,\varepsilon}\left(\omega^{\frac{2(p-1)}{-Np+N+\alpha+2}}\right),\quad \overline{\mu}_{\lambda,\varepsilon}=\omega^{\frac{2(p-1)}{-Np+N+\alpha+2}}\mu_{\lambda,\varepsilon}\\
\end{align*} Next, we give the proof of Theorem 1.2. For fixed ɛ small, we consider a sequence 
$\lambda_n\rightarrow\infty$. Let 
$(\mu_n,u_n):=(\mu_{\lambda_n,\varepsilon},u_{\lambda_n,\varepsilon})$ be the positive solution of problem (1.4). Similar to Proposition 3.2, we get 
$\|u_n\|_{\varepsilon}, \int_{\mathbb{R}^N}\lambda_n\varepsilon^2V_{\varepsilon}(x)u_{\varepsilon}^2$ and µn are bounded. This shows that, up to a subsequence, we may assume that 
$u_n\rightharpoonup u_{\varepsilon}$ in Eɛ and 
$\mu_n\rightarrow\mu_{\varepsilon}$ in 
$\mathbb{R}$. Using Fatou’s lemma, we have
\begin{align*}
\int_{\mathbb{R}^N}\varepsilon^2V_{\varepsilon}(x)u_{\varepsilon}^2\leq \liminf\limits_{n\rightarrow \infty}\int_{\mathbb{R}^N}\varepsilon^2V_{\varepsilon}(x)u_{n}^2\leq \frac{\liminf\limits_{n\rightarrow \infty}\int_{\mathbb{R}^N}\lambda_n\varepsilon^2V_{\varepsilon}(x)u_{n}^2}{\lambda_n}=0,
\end{align*} which implies 
$u_{\varepsilon}=0$ a.e. in 
$\mathbb{R}^N\setminus\Omega_{\varepsilon}$. Next, we claim 
$u_n\rightarrow u_{\varepsilon}$ in Eɛ. Since un solves
\begin{align*}
-\Delta u_n+\lambda_n \varepsilon^2 V_{\varepsilon}(x)u_n+\mu_nu_n=\left(\int_{\mathbb{R}^N}\frac{|u_n|^p}{|x-y|^{N-\alpha}}\right)|u_n|^{p-2}u_n,\ \text{in}\ \mathbb{R}^N,
\end{align*} we deduce that 
$u_{\varepsilon}\in H^1_0(\Omega_{\varepsilon})$ is the solution of
\begin{align}
-\Delta u_{\varepsilon}+\mu_{\varepsilon}u_{\varepsilon}=\left(\int_{\mathbb{R}^N}\frac{|u_{\varepsilon}|^p}
{|x-y|^{N-\alpha}}\right)|u_{\varepsilon}|^{p-2}u_{\varepsilon},\ \text{in}\ \Omega_{\varepsilon}
\end{align}and
\begin{align}
&\int_{\mathbb{R}^N}|\nabla (u_n-u_{\varepsilon})|^2+\int_{\mathbb{R}^N}\left(\lambda_n\varepsilon^2 V_{\varepsilon}(x)+\mu_n\right)(u_n-u_{\varepsilon})^2\nonumber\\ &\quad=\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}
\frac{|u_n-u_{\varepsilon}|^p|u_n-u_{\varepsilon}|^p}
{|x-y|^{N-\alpha}}+o_n(1).
\end{align}By the Hardy-Littlewood-Sobolev inequality,
\begin{align*}
\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}
\frac{|u_n-u_{\varepsilon}|^p|u_n-u_{\varepsilon}|^p}
{|x-y|^{N-\alpha}}\leq C\|u_n-u_{\varepsilon}\|^{2p}\leq C\|u_n-u_{\varepsilon}\|_{H^1(\mathbb{R}^N\setminus\Omega_{\varepsilon})}^{2p-2}\|u_n-u_{\varepsilon}\|^2.
\end{align*}Combining with (5.21), we have
\begin{align*}
\|u_n-u_{\varepsilon}\|^2 \leq C\|u_n-u_{\varepsilon}\|_{H^1(\mathbb{R}^N\setminus\Omega_{\varepsilon})}^{2p-2}\|u_n-u_{\varepsilon}\|^2.
\end{align*} Controlling ɛ small enough such that 
$C\|u_n-u_{\varepsilon}\|_{H^1(\mathbb{R}^N\setminus\Omega_{\varepsilon})}^{2p-2} \lt 1$, then we have 
$\|u_n-u_{\varepsilon}\|^2=o_n(1)$. Hence, we deduce 
$u_n\rightarrow u_{\varepsilon}$ in Eɛ and 
$|u_{\varepsilon}|_2^2=1$. Clearly, it holds
\begin{align*}
\int_{\mathbb{R}^N}\lambda_n\varepsilon^2V_{\varepsilon}(x)u_n^2= \liminf\limits_{n\rightarrow \infty}\int_{\mathbb{R}^N\setminus\Omega_{\varepsilon}}\lambda_n\varepsilon^2V_{\varepsilon}(u_{n}-u_{\varepsilon})^2\rightarrow0.
\end{align*} Similar to the proof of claim 2 in Proposition 3.3 to 
$u_{\varepsilon}|_{\Omega_{\varepsilon}\setminus\Omega_{L,\varepsilon}}$, we obtain 
$u_{\varepsilon}=0$ in 
$\mathbb{R}^N\setminus\Omega_{\varepsilon}$.
Now, we study the properties of µɛ and uɛ. Note that up to a subsequence 
$\mu_{\varepsilon}\rightarrow\mu_{\ell^{-1}}$. Set
\begin{align*}
Q_{R,\varepsilon}:=\{x\in \Omega_{L,\varepsilon}:\ dist(x, \partial\Omega_{\varepsilon})\leq R\}
\end{align*}for every R > 0. Similar to the proof of Claim 1 in Proposition 3.3, we obtain
\begin{align}
\lim\limits_{\varepsilon\rightarrow0}\int_{Q_{R,\varepsilon}}|u_{\varepsilon}|^q=0
\end{align} for 
$1 \lt q \lt \frac{N+\alpha}{N-2}$. Let 
$\varphi_{\varepsilon,R}\in C^{\infty}(\mathbb{R}^N,[0,1])$ satisfying
\begin{align*}
\begin{cases}
\varphi_{\varepsilon,R}=1,\ \text{in}\ \mathbb{R}^N\setminus\Omega_{L,\varepsilon},\\
\varphi_{\varepsilon,R}=0,\ \text{in}\ \Omega_{L,\varepsilon}\setminus Q_{R,\varepsilon},\\
|\nabla \varphi_{\varepsilon,R}|\leq \frac{2}{R}.
\end{cases}
\end{align*} Clearly, it holds 
$\mu_{\varepsilon} \gt \frac{1}{2}\mu_{\ell^{-1}}$ for ɛ small. Multiplying (5.20) by 
$\varphi_{\varepsilon,R}^2u_{\varepsilon}$ and integrating, we deduce for ɛ small enough,
\begin{align*}
&\int_{\mathbb{R}^N}\left[|\nabla (\varphi_{\varepsilon,R}u_{\varepsilon})|^2+\frac{1}{2}\mu_{\ell^{-1}}(\varphi_{\varepsilon,R}u_{\varepsilon})^2\right]\\
\leq&\int_{\mathbb{R}^N}\left[\nabla u_{\varepsilon}\nabla (\varphi_{\varepsilon,R}^2u_{\varepsilon})+\mu_{\varepsilon}(\varphi_{\varepsilon,R}u_{\varepsilon})^2\right]
+\int_{\mathbb{R}^N}|\nabla \varphi_{\varepsilon,R}|^2u_{\varepsilon}^2\\
=&\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_{\varepsilon}(y)|^p
|u_{\varepsilon}(x)|^p\varphi^2_{\varepsilon,R}(x)}{|x-y|^{N-\alpha}}+\int_{\mathbb{R}^N}|\nabla \varphi_{\varepsilon,R}|^2u_{\varepsilon}^2\\
\leq & \left(\int_{Q_{R,\varepsilon}}|u_{\varepsilon}|^{\frac{2Np}{N-\alpha}}\right)^{\frac{N-\alpha}{2Np}}
\left(\int_{\mathbb{R}^N}|u_{\varepsilon}|^{\frac{2Np}{N-\alpha}}\right)^{\frac{N-\alpha}{2Np}}+\int_{Q_{R,\varepsilon}}|\nabla \varphi_{\varepsilon,R}|^2u_{\varepsilon}^2\\ \leq & \left(\int_{Q_{R,\varepsilon}}|u_{\varepsilon}|^{\frac{2Np}{N-\alpha}}\right)^{\frac{N-\alpha}{2Np}}
\left(\int_{\mathbb{R}^N}|u_{\varepsilon}|^{\frac{2Np}{N-\alpha}}\right)
^{\frac{N-\alpha}{2Np}}+\int_{Q_{R,\varepsilon}}\frac{4}{R^2}u_{\varepsilon}^2.
\end{align*}Combining with the above inequality and (5.22), one has
\begin{equation}
\limsup_{\varepsilon\rightarrow0}\|\varphi_{\varepsilon,R}u_{\varepsilon}\|\leq\frac{C}{R}.
\end{equation} For any sequence 
$\varepsilon_n\rightarrow 0$, we define 
$\widetilde{\lambda}_n$ such that 
$\widetilde{\lambda}_n\varepsilon_n^2\alpha_{\varepsilon_n}\rightarrow+\infty$. Clearly, 
$(1-\varphi_{\varepsilon_n,R})v\in H^1_0(\Omega_{L,\varepsilon_n})$ for any 
$v\in E_{\varepsilon_n}$. By calculating,
\begin{align*}
&\left|\left\langle J'_{\widetilde{\lambda}_n,\varepsilon_n}(u_{\varepsilon_n}),v\right\rangle+\mu_{\varepsilon_n}\langle G'(u_{\varepsilon_n}),v\rangle\right|\\
=&\left|\int_{\mathbb{R}^N}\left(\nabla u_{\varepsilon_n}\nabla v+\mu_{\varepsilon_n} u_{\varepsilon_n}v\right)-\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_{\varepsilon_n}(y)|^p|u_{\varepsilon_n}|^{p-2}u_{\varepsilon_n}v}{|x-y|^{N-\alpha}}\right|\\
=&\left|\int_{\mathbb{R}^N}\left(\nabla u_{\varepsilon_n}\nabla(\varphi_{\varepsilon_n,R} v)+\mu_{\varepsilon_n} u_{\varepsilon_n}\varphi_{\varepsilon_n,R}v\right)\right.\\ &-
\left.\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_{\varepsilon_n}(y)|^p
|u_{\varepsilon_n}(x)|^{p-2}u_{\varepsilon_n}(x)\varphi_{\varepsilon_n,R}v}{|x-y|^{N-\alpha}}\right|\\
=&\left|\int_{\mathbb{R}^N}\left[\nabla (\varphi_{\varepsilon_n,R}u_{\varepsilon_n})\nabla v+\nabla \varphi_{\varepsilon_n,R}\left(\nabla u_{\varepsilon_n}v-u_{\varepsilon_n}\nabla v\right)+\mu_{\varepsilon_n}u_{\varepsilon_n}\varphi_{\varepsilon_n,R}v\right]\right.\\
&\left.-\int_{\mathbb{R}^N}\int_{\mathbb{R}^N}\frac{|u_{\varepsilon_n}(y)|^p
|u_{\varepsilon_n}(x)|^{p-2}u_{\varepsilon_n}(x)\varphi_{\varepsilon_n,R}v}{|x-y|^{N-\alpha}}\right|\\
\leq&C\left(\|\varphi_{\varepsilon_n,R}u_{\varepsilon_n}\|+\frac{1}{R}\|u_{\varepsilon_n}\|+\|u_{\varepsilon_n}\|^2\|
\varphi_{\varepsilon_n,R}u_{\varepsilon_n}\|^2\right)\|v\|_{\varepsilon_n}.
\end{align*} Combining with (5.23) and passing to the limit as 
$R\rightarrow0$, it holds
\begin{equation*}
\lim_{R\rightarrow\infty}\|J'_{\widetilde{\lambda}_n,\varepsilon_n}(u_{\varepsilon_n})+\mu_{\varepsilon_n} G'(u_{\varepsilon_n})\|_{E^*_{\varepsilon_n}}=0.
\end{equation*}Since
\begin{equation*}
\lim_{n\rightarrow\infty}J_{\widetilde{\lambda}_n,\varepsilon_n}(u_{\varepsilon_n})=E(s^0),
\end{equation*} by proposition 3.3 to 
$\mu_{\varepsilon_n}$ and 
$u_{\varepsilon_n}$, we obtain that there exist 
$q_{n,i}\in\Omega_{i,\varepsilon_n}$, 
$u_i\in S_{\ell^{-1}}$ such that
\begin{equation*}
dist (q_{n,i},\partial\Omega_{i,\varepsilon_n})\rightarrow\infty \quad \text{and} \ \|u_{\varepsilon_n}-\sum\limits_{i=1}^{\ell}(\xi_{\varepsilon_n}u_i)(\cdot-q_{n,i})
\|_{{\varepsilon}_n}\rightarrow0\ \text{as}\ n\rightarrow\infty.
\end{equation*}By the following scaling,
\begin{align*}
\overline{u}_{\varepsilon}=\omega^{\frac{2+\alpha}{2(-Np+N+\alpha+2)}}u_{\varepsilon}\left(\omega^{\frac{2(p-1)}{-Np+N+\alpha+2}}\right),\quad \overline{\mu}_{\varepsilon}=\omega^{\frac{2(p-1)}{-Np+N+\alpha+2}}\mu_{\varepsilon}.
\end{align*}we complete the proof of Theorem 1.2.
Acknowledgements
H. Zhang is supported by Hunan Provincial Innovation Foundation for Postgraduate (No. CX20230219) and the Fundamental Research Funds for the Central Universities of Central South University (No. 1053320220488). H. Chen is supported by the National Natural Science Foundation of China (Grant No. 12071486).
