A fountain of positive Bubbles on a Coron's Problem for a Competitive Weakly Coupled Gradient System

TL;DR

This paper constructs solutions with multiple positive bubbles concentrating at a point in a critical elliptic system with a small hole, revealing complex bubble interactions and blow-up behavior.

Contribution

It introduces a new type of non-synchronized bubble solutions in a coupled elliptic system with a shrinking domain hole, using Lyapunov-Schmidt reduction.

Findings

Existence of fountain-like positive bubble solutions

Detailed analysis of bubble interaction balancing boundary effects

Demonstration of complex blow-up phenomena in coupled systems

Abstract

We consider the following critical elliptic system: \begin{equation*} \begin{cases} -\Delta u_i=\mu_i u_i^{3}+\beta u_i^{ } \sum\limits_{j\neq i} u_j^{2} \quad \hbox{in}\ \Omega_\varepsilon \\ u_i=0 \hbox{ on } \partial\Omega_\varepsilon , \qquad u_i>0 \hbox{ in } \Omega_\varepsilon \end{cases}\qquad i=1,\ldots, m, \end{equation*} in a domain $\Omega_\varepsilon \subset \mathbb{R}^4$ with a small shrinking hole $B_\varepsilon(\xi_0)$. For $\mu_i>0$, $\beta<0$, and $\varepsilon>0$ small, we prove the existence of a non-synchronized solution which looks like a fountain of positive bubbles, i.e. each component $u_i$ exhibits a towering blow-up around $\xi_0$ as $\varepsilon \to 0$. The proof is based on the Ljapunov-Schmidt reduction method, and the velocity of concentration of each layer within a given tower is chosen in such a way that the interaction between bubbles of different components balance the interaction of the first bubble of each component with the boundary of the domain, and in addition is dominant when compared with the interaction of two consecutive bubbles of the same component.

Equations

$$\{\begin{array}{l}-\Delta u_i={\mu }_i{u}_i^3+\beta u_i\sum _{j\ne i}{u}_j^{2}in{\Omega }_{\epsilon }\\ u_i=0\; on\partial {\Omega }_{\epsilon },u_i>0\; in{\Omega }_{\epsilon }\end{array}i=1,\dots ,m,$$

$$\{\begin{array}{l}-\Delta u_i={\mu }_i{u}_i^p+\beta u_i^{\frac{p-1}{2}}\sum _{j\ne i}{u}_j^{\frac{p+1}{2}}in\Omega \\ u_i=0\; on\partial \Omega ,u_i>0\; in\; \Omega \end{array}i=1,\dots ,m,$$

$$-\Delta u=u^{p}in\Omega ,u=0on\partial \Omega ,u>0in\Omega .$$

$$-\sqrt{{\mu }_1{\mu }_2}<\beta <\min \{{\mu }_1,{\mu }_2\}\text{or}\beta >\max \{{\mu }_1,{\mu }_2\},$$

$$\Omega \subset {\mathbb{R}}^4\text{ bounded domain, symmetric with respect to one of its points ξ0∈Ω,}$$

$$\{\begin{array}{l}-\Delta u_i={\mu }_i{u}_i^3+\beta u_i\sum _{j\ne i}{u}_j^{2}in{\Omega }_{\epsilon }\\ u_i=0\; on\partial {\Omega }_{\epsilon },u_i>0\; in{\Omega }_{\epsilon }\end{array}i=1,\dots ,m,$$

$$U_{\delta ,\xi }={\alpha }_4\frac{\delta }{{\delta }^2+|x-\xi {|}^2},\delta >0,x,\xi \in {\mathbb{R}}^N$$

$$-\Delta U=U^3,U\in {\mathcal{D}}^{1,2}({\mathbb{R}}^4).$$

$$-\Delta W=-\Delta U_{\delta ,\xi }=U_{\delta ,\xi }^3\text{ in }{\Omega }_{\epsilon },W=0\text{ on }\partial {\Omega }_{\epsilon }.$$

$$u_{i,\epsilon }={\mu }_i^{-\frac{1}{p-1}}\sum _{j\in I_i}{P}_{\epsilon }{U}_{{\delta }_j^{\epsilon },{\xi }_0}+{\varphi }_i^{\epsilon },i=1,\dots ,m$$

$${\delta }_j^{\epsilon }=d_j^{\epsilon }{\epsilon }^{\frac{j}{k+1}}{\left(\log \frac{1}{\epsilon }\right)}^{\frac{1}{2}-\frac{j}{k+1}}\text{ for some }{d}_j^{\epsilon }\to d_j^{\ast },j=1,\dots ,k$$

$$d_j^{\ast }={\Gamma }^{\frac{i}{2(k+1)}}{\left(A^2\tau (0)\right)}^{\frac{i}{2(k+1)}-\frac{1}{2}}{\left(\frac{|\beta |{\alpha }_4^4|{\mathbb{S}}^3|}{k+1}\right)}^{\frac{1}{2}-\frac{i}{k+1}}$$

$$\Vert {\varphi }_i^{\epsilon }{\Vert }_{H_0^1({\Omega }_{\epsilon })}\to 0\text{ as }\epsilon \to 0,i=1,\dots ,m.$$

$$H_{\epsilon ,{\xi }_0}=\{u\in H_0^1({\Omega }_{\epsilon }):u(x)=u(2{\xi }_0-x)\forall x\in {\Omega }_{\epsilon }\}.$$

$$m=2,k⩾2,\text{ and }N=4\text{(so that p=3)}.$$

$$I_1=\{\text{odd numbers between 1 and k}\},I_2=\{\text{even numbers between 1 and k}\}.$$

$$\{\begin{array}{l}-\Delta u_1={\mu }_1{u}_1^3+\beta u_1{u}_2^{2}in{\Omega }_{\epsilon },\\ -\Delta u_2={\mu }_2{u}_2^3+\beta u_1^2{u}_{2}in{\Omega }_{\epsilon }\\ u_1=u_2=0on\partial {\Omega }_{\epsilon },u_1,u_2>0\; in{\Omega }_{\epsilon }.\end{array}$$

$$u_{1,\epsilon }={\mu }_1^{-\frac{1}{2}}\sum _{j\in I_1}{P}_{\epsilon }{U}_{{\delta }_j^{\epsilon },{\xi }_0}+{\varphi }_1^{\epsilon }\text{and}{u}_{2,\epsilon }={\mu }_2^{-\frac{1}{2}}\sum _{j\in I_2}{P}_{\epsilon }{U}_{{\delta }_j^{\epsilon },{\xi }_0}+{\varphi }_2^{\epsilon }$$

$${\delta }_j^{\epsilon }=d_j^{\epsilon }{\epsilon }^{\frac{j}{k+1}}{\left(\log \frac{1}{\epsilon }\right)}^{\frac{1}{2}-\frac{j}{k+1}}\text{ for some }{d}_j^{\epsilon }\to d_j^{\ast },j=1,\dots ,k$$

$$d_j^{\ast }={\Gamma }^{\frac{j}{2(k+1)}}{\left(A^2\tau (0)\right)}^{\frac{j}{2(k+1)}-\frac{1}{2}}{\left(\frac{|\beta |{\alpha }_4^4|{\mathbb{S}}^3|}{k+1}\right)}^{\frac{1}{2}-\frac{j}{k+1}}$$

$$\Vert {\varphi }_1^{\epsilon }{\Vert }_{H_0^1({\Omega }_{\epsilon })}\to 0,\Vert {\varphi }_2^{\epsilon }{\Vert }_{H_0^1({\Omega }_{\epsilon })}\to 0\text{ as }\epsilon \to 0.$$

$${\mu }_1={\mu }_2=1,\text{ and take }{\xi }_0=0\in \Omega ,$$

$$J_{\epsilon }(u_1,\dots ,u_m)=\sum _{i=1}^m{\int }_{{\Omega }_{\epsilon }}\left(\frac{|\nabla u_i{|}^2}{2}-\frac{{\mu }_i(u_i^{+}{)}^{p+1}}{p+1}\right)-\frac{2\beta }{p+1}{\sum _{i,j=1}}_{i<j}^m{\int }_{{\Omega }_{\epsilon }}|u_i{|}^{\frac{p+1}{2}}|u_j{|}^{\frac{p+1}{2}}.$$

$$-\Delta u_i={\mu }_i(u_i^{+}{)}^p+\beta \sum _{j\ne i}{u}_i|u_i{|}^{\frac{p-3}{2}}|u_j{|}^{\frac{p+1}{2}},i=1,\dots ,m.$$

$$0⩾-{\int }_{{\Omega }_{\epsilon }}|\nabla u_i^{-}{|}^2=-\beta \sum _{j\ne i}{\int }_{{\Omega }_{\epsilon }}|u_i^{-}{|}^{\frac{p+1}{2}}|u_j{|}^{\frac{p+1}{2}}⩾0.$$

$$\{\begin{array}{l}-\Delta u_i={\mu }_i{u}_i^p+{\beta }_i{u}_i^{q_i}\sum _{j\ne i}{u}_j^{q_j}in{\Omega }_{\epsilon }\\ u_i=0\; on\partial {\Omega }_{\epsilon },u_i>0\; in{\Omega }_{\epsilon }\end{array}i=1,\dots ,m,$$

$$U_{\delta ,0}(x)={\alpha }_4\frac{\delta }{{\delta }^2+|x{|}^2}$$

$$-\Delta V=3U_{\delta }^{2}V$$

$$\frac{\partial U_{\delta }}{\partial \delta }(x)={\alpha }_4\frac{|x{|}^2-{\delta }^2}{({\delta }^2+|x{|}^2{)}^2},\frac{\partial U_{\delta }}{\partial {\xi }_i}(x)=2{\alpha }_4\frac{\delta x_i}{({\delta }^2+|x{|}^2{)}^2},i=1,\dots ,4.$$

$${\mathcal{D}}_s^{1,2}({\mathbb{R}}^4):=\{\psi \in {\mathcal{D}}^{1,2}({\mathbb{R}}^4):\psi (-x)=\psi (x)\forall x\in {\mathbb{R}}^4\}$$

Full text

[TABLE]

in a domain $\Omega_{\varepsilon}\subset\mathbb{R}^{4}$ with a small shrinking hole $B_{\varepsilon}(\xi_{0})$. For $\mu_{i}>0$, $\beta<0$, and $\varepsilon>0$ small, we prove the existence of a non-synchronized solution which looks like a fountain of positive bubbles, i.e. each component $u_{i}$ exhibits a towering blow-up around $\xi_{0}$ as $\varepsilon\to 0$. The proof is based on the Ljapunov-Schmidt reduction method, and the velocity of concentration of each layer within a given tower is chosen in such a way that the interaction between bubbles of different components balance the interaction of the first bubble of each component with the boundary of the domain, and in addition is dominant when compared with the interaction of two consecutive bubbles of the same component.

Acknowlegments. A. Pistoia is partially supported by Sapienza research grant “Nonlinear PDE’s in geometry and physics”. N. Soave is partially supported by the PRIN-2015KB9WPT$\_$010 Grant: “Variational methods, with applications to problems in mathematical physics and geometry”. H. Tavares is partially supported by the Portuguese government through FCT - Fundação para a Ciência e a Tecnologia, I.P., both under the project PTDC/MAT-PUR/28686/2017 and through the grant UID/MAT/04561/2013. The second and third author are also supported by the ERC Advanced Grant 2013 n. 339958 “Complex Patterns for Strongly Interacting Dynamical Systems - COMPAT”

1. Introduction

This paper deals with the existence of solutions to the elliptic critical system

[TABLE]

when $\Omega$ is a bounded smooth domain in $\mathbb{R}^{N}$, and $p=\frac{N+2}{N-2}=2^{*}-1$, with $2^{*}$ critical Sobolev exponent. Thinking at $u_{i}$ as a density function (which is natural since (1.1) is studied in connection with problems in nonlinear optics and Bose-Einstein condensation), the sign of the real parameters $\mu_{i}$ describes the self-interaction between particles of the same density $u_{i}$, and will always be positive: that is, we have attractive self-interaction. On the contrary, the coupling parameter $\beta$, which describes the interaction between particles of different densities, will always be negative: that is, we have repulsive mutual interaction.

The system (1.1) has the trivial solution, i.e. all the components $u_{i}$ vanish. It can also have a semi-trivial solution, i.e. only $\ell<m$ components vanish. It is clear that in this case (1.1) reduces to a system with $m-\ell$ nontrivial components, so we are naturally lead to find fully nontrivial solutions, namely solutions where all the components are nontrivial. In fact, we will be concerned with positive solutions, namely fully nontrivial solutions with $u_{i}>0$ for every $i$.

It is useful to point out that (1.1) can have solutions with synchronized components, i.e. all the components satisfy $u_{i}=s_{i}u$ for some $s_{i}\in\mathbb{R}$ and $u$ solves the single equation

[TABLE]

For instance, if the number of components is $m=2$, the space dimension is $N=4$ (so that $p=3$), and

[TABLE]

then a solution of (1.2) gives rise to a synchronized solution. In this way, results available for the single equation can be translated in terms of (1.1): for instance, if $\Omega$ has nontrivial $\mathbb{Z}_{2}-$homology, then the celebrated Bahri-Coron’s result [2] claims the existence of a positive solution for (1.2), and in turn this gives existence of a synchronized solution for (1.1). It is worthwile to recall also the Coron’s result [11], where the case of a domain with a small hole has been considered, namely $\Omega$ is replaced by $\Omega_{\varepsilon}:=\Omega\setminus\overline{B_{\varepsilon}(\xi_{0})}$, and problem (1.2) has a solution which blows-up at $\xi_{0}$ as $\varepsilon\to 0$ (see also [16, 22]). Again, this family of solutions can be used to construct an associated family of synchronized solutions for (1.1).

The assumptions on the domain are natural, since, exactly as in the scalar case, a Pohozaev-type identity shows that there is no solution if $\Omega$ is starshaped (see for instance [7, p. 519] or [8]).

The above discussion induced the first two authors to investigate the following problem: does (1.1) have non-synchronized solutions? An affirmative answer is given in [20], where (1.1) is posed in a domain $\Omega_{\varepsilon}\subset\mathbb{R}^{N}$, with $N=3,4$, having $\kappa$ distinct holes; that is, $\Omega_{\varepsilon}:=\Omega\setminus\cup_{i=1}^{\kappa}\overline{B_{\varepsilon}(\xi_{i})}$, with $2\leqslant\kappa\leqslant m$; for a quite general choice of interaction terms $\beta_{ij}$ (which can be both of cooperative type, and of competitive type), Pistoia and Soave proved existence and concentration results of solutions whose components are splitted in several groups $G_{1},...,G_{\kappa}$, in such a way that each component within a given group $G_{i}$ concentrates around a point $\xi_{i}$ in a somehow synchronized fashion (in the sense that the velocity of concentration of different components belonging the same group is the same), while the different groups concentrate around different points. In particular, the main results in [20] regard the case when at least two components concentrate around different points, and hence cannot be synchronized.

In view of the above discussion, it is natural to ask the following question: if the domain has only one small hole, is it still possible to find a non-synchronized solution? The main purpose of this paper is to give a positive answer for $\beta<0$ and $N=4$ - so that $p=3$ (for a discussion of the cases $N=3$ or other dimensions, see Remark 1.9 below). More precisely, we take

[TABLE]

i.e. $x\in\Omega$ if and only if $2\xi_{0}-x\in\Omega$, and consider the following elliptic problem with $m\in\mathbb{N}$ equations:

[TABLE]

where $\Omega_{\varepsilon}$ is a domain with one hole, $\Omega_{\varepsilon}:=\Omega\setminus\overline{B_{\varepsilon}(\xi_{0})}\subset\mathbb{R}^{4}$, and $B_{\varepsilon}(\xi_{0})$ denotes the open ball of $\mathbb{R}^{4}$ centered at $\xi_{0}$ with radius $\varepsilon$. Throughout this paper we take $\mu_{i}>0$, the so called focusing case, and $\beta<0$, which means that the coupling terms in (1.4) are of competitive type.

We find solutions of (1.4) which look like a fountain of bubbles, namely their components are a superposition of bubbles centered at $\xi_{0}$ with different rates of concentration. In particular, all the components have a towering blow-up point at $\xi_{0}.$ This new phenomena is quite surprising, since it is in sharp contrast with the case of the single equation for which positive solutions cannot have neither clustering or towering blow-up points, i.e. at every blow-up point there is at most one bubble concentrating there (see Schoen [23]). We also mention that it is somehow unexpected that in a competitive regime (with a possibly large $|\beta|$) we find solutions whose components concentrate at the same point; this is only possible because the concentration rates are different, and in particular such solutions are not synchronized.

In order to state our results we need to introduce some notations. We define

[TABLE]

(a bubble) with $\alpha_{4}=2\sqrt{2}$: these functions are all the positive solutions of the problem

[TABLE]

(see [1, 4, 24]). Also, we denote by $P_{\varepsilon}:\mathcal{D}^{1,2}(\mathbb{R}^{4})\to H^{1}_{0}(\Omega_{\varepsilon})$ the projection map and we define the projection of the bubble defined in (1.5) as $W:=P_{\varepsilon}U_{\delta,\xi}\in H^{1}_{0}(\Omega_{\varepsilon})$, which is the unique solution of

[TABLE]

Take $k\in\mathbb{N}$ (the total number of bubbles) larger than or equal to $m$ . Consider $I_{1},\ldots,I_{m}\subset\{1,\ldots,k\}$ satisfying the following properties:

1. (1)

$1\in I_{1}$; 2. (2)

$I_{i}\neq\emptyset$ for every $i=1,\ldots,m$; 3. (3)

$I_{i}\cap I_{j}=\emptyset$ whenever $i\neq j$; 4. (4)

$I_{1}\cup\ldots\cup I_{m}=\{1,\ldots,k\}$; 5. (5)

for every $j\in\{1,\ldots,k\}$ and $i\in\{1,\ldots,m\}$, if $j\in I_{i}$ then $j-1,j+1\not\in I_{i}$.

Observe that one considers condition (1) without loss of generality, simply to fix ideas and simplify some statements. Conditions (2)-(3)-(4) imply that $I_{1},\ldots,I_{m}$ form a partition of $\{1,\ldots,k\}$, while condition (5) means that each set $I_{i}$ does not contain two consecutive integers. Our main result is the following

Theorem 1.1**.**

Take $\Omega$ satisfying (1.3) and let $\mu_{i}>0$, $\beta<0$. For any integer $k\geqslant m$ and for every partition $I_{1},\ldots,I_{m}$ of $\{1,\ldots,k\}$ satisfying (1)–(5), there exists $\varepsilon_{0}>0$ such that for any $\varepsilon\in(0,\varepsilon_{0})$ problem (1.4) has a solution (symmetric with respect to $\xi_{0}$) of the form

[TABLE]

with

[TABLE]

for

[TABLE]

(see the upcoming (1.16) and (1.17) for the expressions of the constants $A,\Gamma$) and

[TABLE]

Remark 1.2**.**

As stated in the theorem, each component of the solution, $u_{i,\varepsilon}$, belongs to the space

[TABLE]

Since also $U_{\delta,\xi}$ is symmetry with respect to $\xi_{0}$, then $P_{\varepsilon}U_{\delta_{j}^{\varepsilon},\xi_{0}}\in H_{\varepsilon,\xi_{0}}$, as well as the remainder terms $\phi_{i}^{\varepsilon}$.

In order to better explain our result, let us take a particular case of (1.4) and Theorem 1.1:

[TABLE]

and the following partition of $\{1,\ldots,k\}$:

[TABLE]

Clearly, $I_{1},I_{2}$ satisfies conditions (1)–(5), and it is actually the only admissible partition for $m=2$. Problem (1.4) now reads as

[TABLE]

In this particular situation, Theorem 1.1 can be stated in the following way.

Theorem 1.3**.**

Take $\Omega$ satisfying (1.3) and let $\mu_{1},\mu_{2}>0$, $\beta<0.$ For any integer $k\geqslant 2$, let $I_{1},I_{2}$ be respectively the set of all odd and even numbers between 1 and k, as in (1.9). Then there exists $\varepsilon_{0}>0$ such that for any $\varepsilon\in(0,\varepsilon_{0})$ problem (1.10) has a solution (symmetric with respect to $\xi_{0}$) of the form

[TABLE]

with

[TABLE]

for

[TABLE]

and

[TABLE]

In order to avoid insignificant technicalities that would make the presentation harder to follow, we will simply prove Theorem 1.3; in order to convince the reader that the proof of Theorem 1.4 follows precisely in the same way we will make some remarks along the paper (see Remarks 2.2, 3.7, 4.8 and 5.2).

Our result is inspired by the construction performed by Musso and Pistoia in [18] and Ge, Musso and Pistoia in [12], where the authors built sign-changing solutions to Coron’s problem whose shape resembles a superposition of bubbles centered at the point $\xi_{0}$ with alternating sign and with different rate of concentration. The proof here also follows the same scheme which is based on a Ljapunov-Schmidt procedure: we find a good first order approximation term (see (3.2)), we perform a linear theory for the linearized system around the ansatz (see Proposition 3.2), we reduce the problem to a finite dimensional one (see Proposition 3.1) and finally we study the reduced problem (see Section 4). However, the main steps of our proof require rather delicate and careful estimates, see for instance the estimates involving the interacting term in the study of the linear part in Subsection 3.1, the asymptotic expansion of the interaction energy (Lemma 4.4), and the estimate of the remainder term in Lemma 4.6. Indeed, the interaction between bubbles of different components has to balance the interaction of the first bubble of each component with the boundary of the domain, and most of all it has to be dominant compared with the interaction of two consecutive bubbles of the same component. Actually, this is possible because of the presence of an $|\log\varepsilon|$-order term which turns out to be crucial in our construction (see estimate (4.30)).

Remark 1.4**.**

We prove the existence of solutions which look like fountains of positive bubble all centered at the point $\xi_{0}$ when $\Omega$ is symmetric with respect $\xi_{0}.$ It is clear that using the same arguments of Ge, Musso and Pistoia [12] we can remove the symmetry assumption, just centering all the bubbles $U_{\delta_{i},\xi_{i}}$ at suitable points $\xi_{i}=\xi_{i}(\varepsilon)$ which approach $\xi_{0}$ with a suitable rate as $\varepsilon\to 0.$

Remark 1.5**.**

For the sake of completeness, we also mention some recent results concerning the existence of solutions to system (1.1) when $\Omega$ is the whole space $\mathbb{R}^{N}$. As far as we know, all the results deal with systems with only two components. Guo, Li and Wei in [15] established the existence of infinitely many positive nonradial solutions of (1.1), only when $N=3,$ in the competitive case. Peng, Peng and Wang discussed in [19] uniqueness of the least energy solution for $\beta>0$, and the non-degeneracy of the manifold of the synchronized positive solutions. Clapp and Pistoia in [10] proved that system (1.1) in any dimension has infinitely many fully nontrivial solutions, which are not conformally equivalent. Gladiali, Grossi and Troestler in [14, 13] obtained radial and nonradial solutions to some critical systems like (1.1) using bifurcation methods.

Remark 1.6**.**

A Brezis-Nirenberg type problem has been studied for systems, see for instance [7, 8, 6] for existence results, while for concentration and blow-up type results see [5, 21].

Remark 1.7**.**

As already mentioned, appropriate assumptions on $\beta$ allows to obtain a synchronized solution to (1.1) if $\Omega$ has nontrivial $\mathbb{Z}_{2}-$homology. We conjecture that system (1.1) has at least one (actually we would say infinitely many) positive non-synchronized solution if $\Omega$ has nontrivial $\mathbb{Z}_{2}-$homology (as in Bahri-Coron’s result for the single equation (1.2)) and $\beta<0$ is arbitrary. A first attempt in this direction is due to Clapp and Faya [9], who establish the existence of a prescribed number of fully nontrivial solutions to the system with only two components under suitable symmetry assumptions on the topologically nontrival domain $\Omega.$

We would like to remark that the difficulty in finding positive solutions to system (1.1), even with only two components, is similar to the difficulty in finding sign-changing solutions for the single equation (1.2). One key point is the blow-up analysis of solutions: in the case of positive solutions the blow-up, whenever it occurs, is isolated and simple, while in the case of sign-changing solution multiple bubbling naturally appears.

Without loss of generality, we will work from now on with

[TABLE]

assuming that $\Omega$ is symmetric with respect to the origin. Observe that we are conduced to such situation by eventually replacing $u_{i}$ with $\mu_{i}^{-\frac{1}{2}}u_{i}(x+\xi_{0})$.

Remark 1.8**.**

Solutions of (1.4) correspond to critical points with nontrivial components of the $C^{1}$–energy functional $J_{\varepsilon}:H^{1}_{0}(\Omega;\mathbb{R}^{m})\to\mathbb{R}$ defined by

[TABLE]

Indeed, if $(u_{1},\ldots,u_{m})$ is a critical point of $J_{\varepsilon}$, then it satisfies

[TABLE]

Multiplying this equation by $u_{i}^{-}$ and integrating by parts yields (since $\beta<0$)

[TABLE]

If $u_{i}\not\equiv 0$, then by the maximum principle we deduce that $u_{i}>0$.

Remark 1.9**.**

The Sobolev critical exponent is defined only for $N\geqslant 3$. On the other hand, for $p$ defined as before, the right hand sides of (1.1) are $C^{1}$ nonlinearities if and only if we have $\frac{p-1}{2}\geqslant 1$, if and only if $N\leqslant 4$. Therefore, it is reasonable to work in dimension $N=3$ or $N=4$. Here we chose to deal with the case $N=4$ only since it requires less technicalities: all the exponents are positive integers, which makes some expansions explicit. Using Taylor expansions we could have takled the case $N=3$. We conjecture that in this case the main results (and in particular the rates) would be the same.

Remark 1.10**.**

A similar approach could also be used to find solutions for critical systems in pierced domains when the interaction term is more in general like (e.g. Lotka-Volterra systems)

[TABLE]

when $\mu_{i}>0,$ $\beta_{i}<0$ and $q_{i},q_{j}>1.$ In the non-variational cases, one has to replace the asymptotic estimates on the energy of Section 4 with an argument that simply uses the system like in [17, Section 2].

Notations

Working with dimension $N=4$, we deal with the following bubbles concentrated at the origin

[TABLE]

(where $\alpha_{4}=2\sqrt{2}$), which we denote also by $U_{\delta}$; in many cases we deal with different concentration parameters $\delta_{i}$, $i=1,\dots,k$, and we shall simply write $U_{\delta_{i}}=U_{i}$. These correspond to all positive solutions of $-\Delta U=U^{3}$ in $\mathbb{R}^{4}$ which are symmetric with respect to the origin. It is well known (see [3]) that the space of solutions of the linearized equation

[TABLE]

has dimension $4+1=5$ in $\mathcal{D}^{1,2}(\mathbb{R}^{5})$, being spanned by

[TABLE]

Therefore, the space of solutions to (1.14) which belong to

[TABLE]

has dimension 1, being spanned by $\frac{\partial U_{\delta}}{\partial\delta}$. For future convenience, we observe that

[TABLE]

We take the following inner product and norm in $H^{1}_{0}(\Omega_{\varepsilon})$:

[TABLE]

and the standard $L^{p}$ norm by $\|\cdot\|_{p}$ (we omit the dependence on $\varepsilon$ for simplicity).

The Green function of the Laplace operator in $\Omega$ with Dirichlet boundary conditions is denoted by $G(x,y)$, and can be decomposed as

[TABLE]

where $\gamma_{4}:=(2|\partial B_{1}|)^{-1}$, and $H$ is the regular part of $G$ which, for every $x\in\Omega$, satisfies

[TABLE]

The Robin function of $\Omega$ is defined as $\tau(x):=H(x,x)$, and satisfies $\tau(x)\to+\infty$ as $\mathrm{dist}(x,\partial\Omega)\to 0$. Throughout the paper, we will always label the following constants:

[TABLE]

[TABLE]

and use $B_{\varepsilon},\partial B_{\varepsilon}$ instead of $B_{\varepsilon}(0),\partial B_{\varepsilon}(0)$ respectively. We will denote the $L^{p}(\Omega_{\varepsilon})$ norms by $\|\cdot\|_{L^{p}}$, while $\|u\|_{H^{1}_{0}}^{2}:=\int_{\Omega_{\varepsilon}}|\nabla u|^{2}$ for every $u\in H^{1}_{0}(\Omega_{\varepsilon})$.

2. The ansatz and reduction scheme

Recall that, without loss of generality, we assume (1.12); due to the symmetry, by the principle of symmetric criticality we can work in the space

[TABLE]

We deal with solutions of

[TABLE]

where $f:\mathbb{R}\to\mathbb{R}$, $f(s):=(s^{+})^{3}$. Denote by $\mathcal{I}^{*}:L^{\frac{4}{3}}(\Omega_{\varepsilon})\to H^{1}_{0}(\Omega_{\varepsilon})$ the adjoint operator of the canonical Sobolev embedding $\mathcal{I}:H^{1}_{0}(\Omega_{\varepsilon})\to L^{4}(\Omega_{\varepsilon})$. This means that $v:=\mathcal{I}^{*}u$ can be defined as the (unique) weak solution of

[TABLE]

Observe that, if $u$ is symmetric with respect to the origin, so is $\mathcal{I}^{*}u$. The operator $\mathcal{I}^{*}$ is continuous: there exists $C>0$, independent of $\varepsilon$, such that

[TABLE]

Using this operator, we can rewrite (2.2) as

[TABLE]

Denote $U_{j}:=U_{\delta_{j}}$ for $j=1,\ldots,k$. Our ansatz is the following: for any integer $k\geqslant 2$, we look for a solution of (2.2) in $H_{\varepsilon}$ of the form

[TABLE]

where

[TABLE]

$\mathbf{d}=(d_{1},\ldots,d_{k})$ belongs to the set

[TABLE]

and $\phi_{1},\phi_{2}\in H_{\varepsilon}$.

Remark 2.1**.**

For future reference, we collect in this remark several important relations between the different rates $\delta_{j}$. Given $\eta>0$, we have

[TABLE]

as $\varepsilon\to 0$, uniformly for $\mathbf{d}\in X_{\eta}$.

For each $\varepsilon>0$ small, our aim is to find $\eta>0$, $\mathbf{d}\in X_{\eta}$ and $\phi_{1},\phi_{2}\in H_{\varepsilon}$ such that, for $i,j=1,2$, $i\neq j$,

[TABLE]

Given $\varepsilon>0$ and $d_{1},\ldots,d_{k}>0$, for $\delta_{i}$ defined as before define

[TABLE]

(recall the Notation section) and

[TABLE]

Observe that $\mathbf{K}_{\mathbf{d},\varepsilon}^{\perp}=K_{1}^{\perp}\times K_{2}^{\perp}$. Moreover, consider the projection maps

[TABLE]

We can rewrite (2.11) as a system of 4 equations: for $i,j=1,2$, $j\neq i$,

[TABLE]

[TABLE]

In the next section, given $\varepsilon,\eta>0$ sufficiently small and $\mathbf{d}\in X_{\eta}$, we find a unique $(\phi_{1},\phi_{2})=(\phi_{1}^{\mathbf{d},\varepsilon},\phi_{2}^{\mathbf{d},\varepsilon})\in K_{\mathbf{d},\varepsilon}^{\perp}$ solution to (2.10). By plugging this result in (2.9), we end up having a problem with unknown $\mathbf{d}\in\mathbb{R}^{k}$ (thus a finite dimensional problem), which can be stated in terms of a reduced energy. We analyse this reduced energy in Section 4.

Remark 2.2**.**

For the general system (1.4) and given a partition $I_{1},\ldots,I_{m}$ of $\{1,\ldots,k\}$, the ansatz is exactly the same: $u_{i}=\sum_{j\in I_{i}}U_{i}+\phi_{i}$, for $i=1,\ldots,m$, where $\phi_{i}\in K_{i}$. We denote in this case $\mathbf{K}_{\mathbf{d},\varepsilon}=K_{1}^{\perp}\times\ldots\times K_{m}^{\perp}$, and split the system of $m$ equations:

[TABLE]

($i=1,\ldots,m$) in $2m$ equations using the projection maps $\Pi_{i}$ and $\Pi_{i}^{\perp}$.

3. Reduction to a Finite Dimensional Problem

In this section we study the solvability of (2.10). We rewrite (2.10) as

[TABLE]

where $L$ stays for the linear part

[TABLE]

$N$ stays for the nonlinear part

[TABLE]

and $R$ is the remainder term

[TABLE]

where the last equality is a consequence of the definitions of $\mathcal{I}^{*}$ and of $f$ (analogue expressions hold for $L^{2}_{\mathbf{d},\varepsilon}$, $N^{2}_{\mathbf{d},\varepsilon}$ and $R^{2}_{\mathbf{d},\varepsilon}$).

We also define

[TABLE]

and $\mathbf{R}_{\mathbf{d},\varepsilon}$ and $\mathbf{N}_{\mathbf{d},\varepsilon}$ in an analogue way.

Proposition 3.1**.**

Let $\beta<0$. Then for every $\eta>0$ sufficiently small there exists $\varepsilon_{0}>0$ and $C>0$ such that, whenever $\varepsilon\in(0,\varepsilon_{0})$ and $\mathbf{d}\in X_{\eta}$, there exists a unique function $\boldsymbol{\phi}=\boldsymbol{\phi}^{\mathbf{d},\varepsilon}\in\mathbf{K}^{\perp}_{\mathbf{d},\varepsilon}$ solving the equation

[TABLE]

and satisfying

[TABLE]

Moreover, the map $X_{\eta}\to\mathbf{K}^{\perp}_{\mathbf{d},\varepsilon},$ $\mathbf{d}\mapsto\boldsymbol{\phi}^{\mathbf{d},\varepsilon}$ is of class $\mathcal{C}^{1}.$

The proof of the proposition takes the rest of this section, and is divided into several intermediate lemmas.

3.1. Study of the linear part

As a first step, it is important to understand the solvability of the linear problem associated with (3.1), i.e.

[TABLE]

Proposition 3.2**.**

For every $\eta>0$ small enough there exists $\varepsilon_{0}>0$ small, and $C>0$, such that if $\varepsilon\in(0,\varepsilon_{0})$ then

[TABLE]

for every $\mathbf{d}\in X_{\eta}$. Moreover, $\mathbf{L}_{\mathbf{d},\varepsilon}$ is invertible in $\mathbf{K}_{\mathbf{d},\varepsilon}^{\perp}$, with continuous inverse.

The long proof proceeds by contradiction. For a fixed $\eta>0$ small, let us suppose that there exist sequences

[TABLE]

such that

[TABLE]

as $n\to\infty$, where we wrote $K_{i,n}:=K_{i,\mathbf{d}_{n},\varepsilon_{n}}$ and $\mathbf{L}_{n}:=\mathbf{L}_{\mathbf{d}_{n},\varepsilon_{n}}$ for short. In the same spirit, in this proof we write $P_{n}:=P_{\varepsilon_{n}}$, $U_{i,n}:=U_{\delta_{i,n},0}$, $\psi_{i,n}:=\psi_{\delta_{i,n},0}$, and $\Omega_{n}:=\Omega_{\varepsilon_{n}}$.

Let $\mathbf{h}_{n}:=\mathbf{L}_{n}(\boldsymbol{\phi}_{n})$. Then, by definition of $\mathbf{L}_{n}$,

[TABLE]

(an analogue equation holds for $\phi_{2,n}$) for some $w_{i,n}\in K_{i,n}$.

Lemma 3.3**.**

$\|w_{i,n}\|_{H_{0}^{1}(\Omega_{n})}\to 0$* as $n\to\infty$.*

Proof.

We focus on $w_{1,n}$, the proof for $w_{2,n}$ is analogue. As $w_{1,n}\in K_{1,n}=\textrm{span}\{P_{n}\psi_{j,n}:j\in I_{1}\}$, there exist constants $c_{j,n}$ such that

[TABLE]

Now we consider the scalar product in $H_{0}^{1}(\Omega_{\varepsilon})$ of both sides in (3.7) with $\delta_{i,n}P_{n}\psi_{i,n}$, with $i\in I_{1}$: as $h_{1,n},\phi_{1,n}\in K_{1,n}^{\perp}$, we obtain

[TABLE]

The left hand side can be estimated using [12, Remark 5.2] and (1.15) (see also [21, p. 417], noting that therein $\psi_{i,n}$ corresponds to $\delta_{i,n}\psi_{i,n}$ in the present paper) and obtaining

[TABLE]

as $n\to\infty$, where

[TABLE]

The first integral on the right hand side in (3.8) can be estimated as in [12, Formula (5.7)]:

[TABLE]

as $n\to\infty$. We have now to estimate the interaction terms. To this purpose, we observe that by Hölder and Sobolev inequality, and by (1.15),

[TABLE]

as $n\to\infty$. The precise rate of the higher order terms ($h.o.t.$) does not play any role, and in any case can be derived using Lemmas A.1 and A.2. Moreover, the leading integral on the right hand side can be estimated using Lemma A.4, obtaining

[TABLE]

Coming back to (3.11), we have

[TABLE]

as $n\to\infty$, which proves that the second integral on the right hand side in (3.8) is of order $o(\delta_{i,n})$. As far as the third integral is concerned, we note that

[TABLE]

as $n\to\infty$. The last inequality follows by Lemma A.6 if $h\neq i$, and by Lemma A.4 if $h=i$. In any case

[TABLE]

as $n\to\infty$. To sum up, by expanding (3.8), we proved that for every index $i\in I_{1}$ it results that

[TABLE]

as $n\to\infty$. From this and by Cramer’s rule, we deduce that $c_{i,n}\to 0$ for every $i\in I_{1}$. From this, the conclusion $\|w_{1,n}\|\to 0$ follows. ∎

Let us set now $z_{i,n}:=\phi_{i,n}-h_{i,n}-w_{i,n}$. Notice that, since $\|h_{i,n}\|_{H^{1}_{0}(\Omega_{n})},\|w_{i,n}\|_{H^{1}_{0}(\Omega_{n})}\to 0$, we have $\|z_{1,n}\|_{H^{1}_{0}(\Omega_{n})}^{2}+\|z_{2,n}\|_{H^{1}_{0}(\Omega_{n})}^{2}\to 1$. In terms of $z_{i,n}$, equation (3.7) can be rewritten as

[TABLE]

Of course, a similar equation holds for $z_{2,n}$.

Lemma 3.4**.**

It results that at least one of the following lower estimates holds:

[TABLE]

or

[TABLE]

Proof.

Since $\|z_{1,n}\|_{H^{1}_{0}(\Omega_{n})}^{2}+\|z_{2,n}\|_{H^{1}_{0}(\Omega_{n})}^{2}\to 1$, we can suppose that up to a subsequence $\{\|z_{1,n}\|_{H^{1}_{0}(\Omega_{n})}^{2}\}_{n}$ or $\{\|z_{2,n}\|_{H^{1}_{0}(\Omega_{n})}^{2}\}_{n}$ is uniformly bounded from below by $1/2$. Suppose for instance that $\{\|z_{1,n}\|^{2}_{H^{1}_{0}(\Omega_{n})}\}_{n}$ is bounded from below. Then we test equation (3.15) with $z_{1,n}$, obtaining

[TABLE]

Arguing as in [12, Formula (5.12)], we can easily check that the last two integrals are [math]. Therefore, in this case the first $\liminf$ in the thesis is positive. If $\{\|z_{2,n}\|_{H^{1}_{0}(\Omega_{n})}^{2}\}_{n}$ is bounded from below, in the same way we find that the second $\liminf$ is positive. ∎

We aim to obtain a contradiction with Lemma 3.4. To this end, we fix $\rho>0$ so that $B_{\rho}\subset\subset\Omega$, and we decompose $B_{\rho}\setminus\overline{B_{\varepsilon_{n}}}$ into the union of disjoint annuli as follows:

[TABLE]

with the convention $\delta_{0,n}=\delta_{1,n}^{-1}\rho^{2}$ and $\delta_{k+1,n}=\delta_{k,n}^{-1}\varepsilon_{n}^{2}$. Recall from Remark 2.1 that $\delta_{l+1,n}/\delta_{l,n}\to 0$ as $n\to\infty$. We also set

[TABLE]

and, for every $\ell=1,\dots,k$, we define a cut-off function $\chi_{\ell,n}\in C^{\infty}_{c}(\mathbb{R}^{N})$ with the properties that

[TABLE]

for a positive universal constant $C$. Finally, we define for $\ell=1,\dots,k$ and $i=1,2$ the $\mathcal{D}^{1,2}(\mathbb{R}^{4})$ functions

[TABLE]

naturally extended by [math] in $\mathbb{R}^{4}\setminus\hat{\mathcal{B}}_{\ell,n}$. We have $\hat{z}_{i,n}^{\ell}(x)=\delta_{\ell,n}z_{i,n}(\delta_{\ell,n}x)$ if $x\in\tilde{\mathcal{A}}_{\ell,n}:=\mathcal{A}_{\ell,n}/\delta_{\ell,n}$.

Lemma 3.5**.**

It results that $\hat{z}_{i,n}^{\ell}\to 0$ weakly in $\mathcal{D}^{1,2}(\mathbb{R}^{4})$, and strongly in $L^{q}_{\mathrm{loc}}(\mathbb{R}^{4})$, for every $q\in[2,2^{*})$, for every $i=1,2$, $\ell=1,\dots,k$.

Proof.

We have

[TABLE]

and

[TABLE]

for $x\in\tilde{\mathcal{B}}_{\ell,n}$, that is,

[TABLE]

Notice that $\tilde{\mathcal{B}}_{\ell,n}$ exhausts $\mathbb{R}^{N}$ as $n\to\infty$, by Remark 2.1. Now

[TABLE]

The integral of $|\nabla z_{i,n}|^{2}$ is clearly bounded, since $\|z_{i,n}\|_{H^{1}_{0}(\Omega_{n})}\leqslant 1$. Also, by (3.17),

[TABLE]

and we infer that $\|\hat{z}_{i,n}^{\ell}\|_{\mathcal{D}^{1,2}(\mathbb{R}^{4})}\leqslant C$. Then, up to a subsequence, we have that $\hat{z}_{i,n}^{\ell}\rightharpoonup\hat{z}_{i}^{\ell}$ weakly in $\mathcal{D}^{1,2}$, and $\hat{z}_{i,n}^{\ell}\to\hat{z}_{i}^{\ell}$ strongly in $L^{q}_{\mathrm{loc}}(\mathbb{R}^{4})$ for $q\in[2,2^{*})$. The equation satisfied by the weak limit can be determined using (3.15) and (3.18): for every $\varphi\in C^{\infty}_{c}(\mathbb{R}^{4}\setminus\{0\})$, by combining (3.15) with (3.18) we have that

[TABLE]

The last integral and all the terms involving $h_{i,n}$ and $w_{i,n}$ tend to [math] as $n\to\infty$, exactly as in [12, Formula (5.20)]. Therefore,

[TABLE]

In order to study the behavior of the integrals as $n\to\infty$, it is convenient to observe (see Lemma A.1) that, if $\ell\in I_{1}$, then

[TABLE]

as $n\to\infty$. If instead $\ell\not\in I_{1}$, then we have a similar expansion, but without the term $U_{1,0}(x)/\delta_{\ell,n}$. We focus at first on the first possibility. We have,

[TABLE]

Now, for every $j<\ell$

[TABLE]

and, by Lemma A.3,

[TABLE]

Moreover, for every $j>\ell$, using the fact that $\mathrm{supp\,}\varphi\subset B_{R}\setminus B_{\rho}$ for suitable $0<\rho<R$, we have that

[TABLE]

and that

[TABLE]

The previous estimates yield

[TABLE]

as $n\to\infty$, for every $\varphi\in C^{\infty}_{c}(\mathbb{R}^{4}\setminus\{0\})$.

Notice that, in the above computations, we never used the fact that the indexes $j$ were in $I_{1}$. Therefore, we directly deduce that

[TABLE]

Finally, in an analogue way

[TABLE]

as $n\to\infty$.

Collecting together (3.21), (3.22) and (3.23), and coming back to (3.20), we finally obtain that the weak limit of $\hat{z}_{1,n}^{\ell}$ satisfies

[TABLE]

Let now $\theta_{\rho}\in C^{\infty}(\mathbb{R}^{N})$ be such that $\theta_{\rho}\equiv 1$ in $B_{2\rho}^{c}$, $\theta_{\rho}\equiv 0$ in $B_{\rho}$, and $|\nabla\theta_{\rho}|\leqslant C/\rho$; and let $\varphi\in C^{\infty}_{c}(\mathbb{R}^{4})$; testing the above equation with $\theta_{\rho}\varphi$, and passing to the limit as $\rho\to 0^{+}$, using the fact that $\hat{z}_{1}^{\ell}\in\mathcal{D}^{1,2}(\mathbb{R}^{4})$ (since it is the weak limit of $\mathcal{D}^{1,2}$ functions), we easily deduce that

[TABLE]

In order to show that $\hat{z}_{1}^{\ell}\equiv 0$, recalling that it is symmetric with respect to [math], it is sufficient to verify that $\hat{z}_{1}^{\ell}\perp\psi_{1,0}$. This can be done exactly as in [12, Formula (5.19)], and completes the proof.

It still remains to analyze the case $\ell\not\in I_{1}$. In such a situation we can proceed exactly as before, but this time we end up with

[TABLE]

Since $\beta<0$ and $\hat{z}_{1}^{\ell}\in\mathcal{D}^{1,2}(\mathbb{R}^{4})$, we infer that

[TABLE]

and the conclusion follows also in this case. ∎

Conclusion of the proof of Proposition 3.2.

Using Lemma 3.5, we will obtain a contradiction with Lemma 3.4. Let us consider

[TABLE]

We show that the right hand side tends to [math] as $n\to\infty$. At first, we have

[TABLE]

Now, let $i\neq\ell$. Then we have

[TABLE]

Since $i\neq\ell$, we have that

[TABLE]

and

[TABLE]

Therefore, the fact that

[TABLE]

follows from the integrability of $U_{1,0}^{4}$ on $\mathbb{R}^{N}$. If moreover $i=\ell$, recalling that $\tilde{\mathcal{A}}_{\ell,n}=\mathcal{A}_{\ell,n}/\delta_{\ell,n}$ we have

[TABLE]

as $n\to\infty$, since $U_{1,0}^{2}\in L^{2}(\mathbb{R}^{N})$ and $(\hat{z}_{1,n}^{\ell})^{2}\rightharpoonup 0$ weakly in $L^{2}(\mathbb{R}^{4})$ by Lemma 3.5. By (3.25) and (3.26), we obtain in (3.24) that

[TABLE]

In a completely analogue way, we also have

[TABLE]

Finally,

[TABLE]

as $n\to\infty$. But (3.24), estimates (3.28), (3.29) and (3.30) imply that

[TABLE]

in contradiction with Lemma 3.4. ∎

3.2. Estimates on the reminder term

In this subsection we prove the following

Proposition 3.6**.**

Let $\eta>0$. There exists $\varepsilon_{0}>0$ and $C>0$ such that

[TABLE]

for $i=1,2$, for every $\textbf{d}\in X_{\eta}$ and for every $\varepsilon\in(0,\varepsilon_{0})$.

Proof.

We focus on $i=1$. By continuity of $\Pi_{1}^{\perp}$ and of $\mathcal{I}^{*}$, there exists $C>0$ such that

[TABLE]

We estimate separately $\mathcal{A}$ and $\mathcal{B}$. At first we note that

[TABLE]

Recalling that $f(s)=(s^{+})^{3}$, and using the fact that

[TABLE]

for a positive constant $C_{n}$ depending only on $n$, we obtain

[TABLE]

Let us fix $j\neq h$. Then, by Lemma A.4,

[TABLE]

Recalling (2.5), we see that if $h>j$ and $\mathbf{d}\in X_{\eta}$

[TABLE]

where $C$ denotes a positive constant depending on $\eta$ (but not on $\mathbf{d}$) and we used the fact that $h-j\geqslant 2$ since $j,h\in I_{1}$ with $j\neq h$. The same estimate holds in case $h<j$. Plugging this into (3.34), and coming back to (3.33), we finally conclude that

[TABLE]

Regarding $A_{2}$, we have

[TABLE]

Using the estimate for $U_{i}-P_{\varepsilon}U_{i}$ contained in Lemma A.1, we deduce that

[TABLE]

In a similar way

[TABLE]

Plugging (3.37) and (3.38) into (3.36), and recalling again Remark 2.1, we obtain

[TABLE]

Therefore, (3.32), (3.35) and (3.39) give

[TABLE]

and it remains to estimate $\mathcal{B}$. By Lemma A.4

[TABLE]

and for any $h>j$ we have

[TABLE]

The same estimate holds for $\delta_{j}/\delta_{h}$ in case $j>h$. Therefore, gathering (3.31), (3.40) and (3.41), we obtain the desired result. ∎

3.3. The nonlinear part: end of the proof of Proposition 3.1

In virtue of Proposition 3.2, solving the equation

[TABLE]

reduces to finding a fixed point of the operator

[TABLE]

in the ball

[TABLE]

for some $\rho>0.$ It is quite standard to show that $\mathbf{T}_{\mathbf{d},\varepsilon}:\mathcal{B}_{\rho}\to\mathcal{B}_{\rho}$ is a contraction mapping for $\varepsilon$ small enough. Indeed, Proposition 3.2 together with straightforward computations lead to

[TABLE]

and

[TABLE]

A standard argument also shows that the map $\mathbf{d}\to\boldsymbol{\phi}^{\mathbf{d},\varepsilon}$ is of class $C^{1}.$

Remark 3.7**.**

Suppose that, instead of dealing with the set of odd and even numbers of $\{1,\ldots,k\}$ in the two equation case, we are dealing with system (1.4) with $m$ equations and with a general partition $I_{1},\ldots,I_{m}$ satisfying (1)–(5). Having already splitted the original problem into $2m$ equations (see Remark 2.2), we can repeat the argument used for $m=2$ without substantial changes, using the fact that each set $I_{j}$ does not contain consecutive integers.

4. Expansion of the reduced energy

Recall that the energy funcional is given by

[TABLE]

where $F(s)=(s^{+})^{4}/4$. Recall that we denote $f(s):=F^{\prime}(s)=(s^{+})^{3}$. For every $\eta>0$ small fixed, we introduce the reduced functional $\widetilde{J}_{\varepsilon}:X_{\eta}\to\mathbb{R}$ as being

[TABLE]

This is a $C^{1}$ functional due to Proposition 3.1 and since $\delta_{i}$ depends on $d_{i}$ via (2.5). Finding critical point of $\widetilde{J}_{\varepsilon}$ corresponds to find solutions of our original system, as we prove next.

Lemma 4.1**.**

Given $\varepsilon\in(0,\varepsilon_{0})$ and $\eta>0$ small, let $\mathbf{d}\in X_{\eta}$. We have that

[TABLE]

if, and only if,

[TABLE]

Proof.

To simplify notations, define $V_{i}^{\mathbf{d},\varepsilon}:=\sum_{j\in I_{i}}P_{\varepsilon}U_{\delta_{j}}+\phi_{i}^{\mathbf{d},\varepsilon}$ for $i=1,2$. From (2.5) we see that

[TABLE]

Hence, if $(V_{1}^{\mathbf{d},\varepsilon},V_{2}^{\mathbf{d},\varepsilon})$ solves (2.2) then $J_{\varepsilon}^{\prime}(V_{1}^{\mathbf{d},\varepsilon},V_{2}^{\mathbf{d},\varepsilon})=0$ and so $\tilde{J}_{\varepsilon}^{\prime}(\mathbf{d})=0$. Conversely, assume $\mathbf{d}\in X_{\eta}$ is a solution of $\widetilde{J}_{\varepsilon}^{\prime}(\mathbf{d})=0$. For $i\in\{1,2\}$ and $l\in I_{i}$, recalling that $\psi_{l}:=\partial_{\delta_{l}}U_{\delta_{l}}$, we have from (4.1) that

[TABLE]

From (2.10), Proposition 3.1 and recalling that $K_{i}$ is spanned by $P_{\varepsilon}\psi_{j}$ for $j\in I_{i}$, we deduce the existence of coefficients $c_{i}^{j}=c_{i}^{j}(\varepsilon,\mathbf{d})$, $j\in I_{i}$ such that

[TABLE]

In conclusion, for $i\in\{1,2\}$ and $l\in I_{i}$,

[TABLE]

A straightforward computation shows that

[TABLE]

for some constant $\sigma_{ll}>0$ (see for instance [21, p. 417]). On the other hand, we have $\langle P_{\varepsilon}\psi_{j},\partial_{\delta_{l}}\phi_{n}^{\mathbf{d},\varepsilon}\rangle=\textrm{o}(1)$. Indeed, since $\phi_{n}^{\mathbf{d},\varepsilon}\in K_{n}^{\perp}$ $(n=1,2)$, then $\langle P_{\varepsilon}\psi_{j},\phi_{n}^{\mathbf{d},\varepsilon}\rangle=0$ for every $\mathbf{d}$. Therefore, taking the derivative of the previous identity with respect to $\delta_{l}$ ($l\in I_{n}$), we get $\langle P_{\varepsilon}\psi_{j},\partial\delta_{l}\phi_{n}^{\mathbf{d},\varepsilon}\rangle=-\langle\partial_{\delta_{l}}P_{\varepsilon}\psi_{j},\phi_{n}^{\mathbf{d},\varepsilon}\rangle$. Combining (1.15) with Lemma A.3 we have $\|\delta_{l}\partial_{\delta_{l}}P_{\varepsilon}\phi_{j}\|=\textrm{O}(1)$, while Proposition 3.1 yields $\|\phi_{n}^{\mathbf{d},\varepsilon}\|=\textrm{o}(\delta_{1})$. Therefore, $\langle P_{\varepsilon}\psi_{j},\delta_{l}\partial\delta_{l}\phi_{n}^{\mathbf{d},\varepsilon}\rangle=o(\delta_{1})=o(1)$, as claimed. In conclusion, we end up with a linear system of the form

[TABLE]

which, as $\varepsilon\to 0$, has the unique solution $c_{i}^{j}=0$ for every $i=1,2$, $j\in I_{i}$. Looking back at (4.2) we see that $(V_{1}^{\mathbf{d},\varepsilon},V_{2}^{\mathbf{d},\varepsilon})$ solves (2.2), as we wanted. ∎

We now compute the leading term of the reduced energy. For simplicity, and when there is no risk of confusion, we denote $\phi_{i}=\phi_{i}^{\mathbf{d},\varepsilon}$ and $U_{i}=U_{\delta_{i}}$. We have

[TABLE]

where

[TABLE]

will be an higher order term.

In what follows we show that the reduced energy reads as

[TABLE]

for some constants $c_{1},c_{2},c_{3},c_{4}>0$. This yields the choice of parameters (2.5) (which for convenience of the reader we recall)

[TABLE]

and the existence of towers of bubbles as we want. Observe that, as $\varepsilon\to 0$,

[TABLE]

(see ahead for the details).

Lemma 4.2**.**

Given $i=1,\ldots,k$ we have

[TABLE]

as $\varepsilon\to 0$, uniformly for every $\mathbf{d}\in X_{\eta}$. We recall that $A,B$ and $\Gamma$ are defined in (1.16)–(1.17), while $\tau$ is the Robin function (see the notation section).

Proof.

We reason similarly to [20, Lemma 4.3], to which we refer for more details.

First of all, using (1.6), we have that

[TABLE]

Using a Taylor expansion up to second order, we have that

[TABLE]

for some function $\xi(x)\in[0,1]$. Therefore, we can rewrite (4.9) as

[TABLE]

We now estimate each one of the three terms separately. The first term in (4.10) is, after a change of variables $x=\delta_{i}y$ and recalling that $\Omega_{\varepsilon}=\Omega\setminus B_{\varepsilon}$ and $\varepsilon/\delta_{i}\to 0$,

[TABLE]

As for the second term, we use the fact that

[TABLE]

(by Lemma A.1, which we can apply since $\varepsilon/\delta_{i}\to 0$ as $\varepsilon\to 0$). We have

[TABLE]

where we have used the estimates for the remainder term $R$ contained in Lemma A.1.

As for the last term in (4.10), since $0\leqslant P_{\varepsilon}U_{i}\leqslant U_{i}$ (by the maximum principle) and $\xi(x)\in[0,1]$, we have $0\leqslant U_{i}+\xi(PU_{i}-U_{i})\leqslant U_{i}$. Combining this with Lemma A.2 and since $\varepsilon/\delta_{i}\to 0$,

[TABLE]

The result follows combining (4.10) with (4.13)–(4.17)–(4.19). ∎

Corollary 4.3**.**

The following estimate holds

[TABLE]

as $\varepsilon\to 0$, uniformly for every $\mathbf{d}\in X_{\eta}$.

Proof.

From the previous lemma we see that

[TABLE]

and the first identity of the lemma follows because $d\in N_{\eta}$ and $\delta_{i}=o(\delta_{1})$ for every $i\in\{2,\ldots,k\}$ (recall Remark 2.1), which implies that $\varepsilon^{2}/\delta_{i}^{2}=\textrm{o}(\varepsilon^{2}/\delta_{k}^{2})$ for $i\in\{1,\ldots,k-1\}$ as $\varepsilon\to 0$. The second identity follows directly from the definition of $d_{i}$ (see (2.5)). ∎

Lemma 4.4**.**

Given $i,j\in\{1,\ldots,k\}$ with $i>j$, we have

[TABLE]

as $\varepsilon\to 0$, uniformly for every $\mathbf{d}\in X_{\eta}$.

Proof.

First, we rewrite

[TABLE]

(by Lemma A.1). We estimate the leading term as follows: for $\varepsilon>0$ small and $r>0$ such that $B_{\varepsilon}\subset B_{r}\subset\Omega$ and $\varepsilon<\sqrt{\delta_{i}\delta_{j}}$,

[TABLE]

Asymptotic estimate of $(I)$: scaling $x=\delta_{i}y$,

[TABLE]

We have:

[TABLE]

The first term can be estimated as follows:

[TABLE]

since, as $\delta_{j}/\delta_{i}\to\infty$ ($i>j$) and $\varepsilon/\delta_{i}\to 0$ (recall Remark 2.1):

[TABLE]

and

[TABLE]

As for the second term, because $(\delta_{i}/\delta_{j})|y|\leqslant\sqrt{\delta_{i}/\delta_{j}}\leqslant c$ for $|y|\leqslant\sqrt{\delta_{j}/\delta_{i}}$, as $\varepsilon\to 0$, and recalling the computation done for $(I.a)$, we have

[TABLE]

Combining the expansions of $(I.a)$ and $(I.b)$ with (4.24) yields, in conclusion, that

[TABLE]

Asymptotic estimate of $(II)$: by using this time the scaling $x=\delta_{j}y$ and the fact that

[TABLE]

we have

[TABLE]

Asymptotic estimate of $(III)$:

[TABLE]

By combining the estimates of $(I)$, $(II)$ and $(III)$ we deduce that

[TABLE]

which yields the desired conclusion. ∎

Corollary 4.5**.**

We have, as $\varepsilon\to 0$, uniformly for every $\mathbf{d}\in X_{\eta}$,

[TABLE]

Proof.

The first identity is a simple consequence of the previous lemma together with the fact that $\delta_{l}=\textrm{o}(\delta_{i})$ as $\varepsilon\to 0$, for $l>i$. In fact, since each one of the sets $I_{1}$ and $I_{2}$ do not contain two consecutive integers, and that given $i>j$ with $|i-j|>1$ it holds

[TABLE]

then

[TABLE]

The last identity of the statement is a consequence of the definition of $\delta_{i}$ and the fact that

[TABLE]

Lemma 4.6**.**

We have

[TABLE]

as $\varepsilon\to 0$, uniformly for every $\mathbf{d}\in X_{\eta}$.

Proof.

Recall that $F(s)=(s^{+})^{4}/4$, and we denote $f(s):=F^{\prime}(s)=(s^{+})^{3}$. We have

[TABLE]

Recalling the definition of $P_{\varepsilon}$ from (1.6) and adding and subtracting terms of type $F(\sum_{j\in I_{i}}P_{\varepsilon}U_{j})$ and $f(\sum_{j\in I_{i}}P_{\varepsilon}U_{j})\phi_{i}$, we have

[TABLE]

Moreover,

[TABLE]

Let us rewrite $R(\mathbf{d},\varepsilon)$ as

[TABLE]

Estimates for $a_{1},a_{2}$. First of all, we check that the first two terms satisfy $a_{1},a_{2}=o(\delta_{1}^{2}).$ Indeed, since $0\leqslant P_{\varepsilon}U_{i}\leqslant U_{i}$ (by the maximum principle), $a_{1}+a_{2}$ is controlled by a sum of terms of the form $\int_{\Omega_{\varepsilon}}P_{\varepsilon}U_{i}P_{\varepsilon}U_{j}P_{\varepsilon}U_{h}P_{\varepsilon}U_{l}$ for indices $j,l,h,i$ not all equal at the same time; each term is of higher order with respect to the leading term $\delta_{1}^{2}$, as we will now check. Indeed, if $i\not=j$ we have by Lemma A.4 that

[TABLE]

because we are always in a situation that $i,j$ belong to the same set $I_{h}$, thus $|i-j|>1$, and then by the choices we did in (2.5),

[TABLE]

Moreover, if $i\not=j$, then assuming without loss of generality that $i>j$ with $|i-j|>1$ we have by Lemma A.5

[TABLE]

(note that this term only appears in $a_{2}$). In a similar way, if $i,j,l\in I_{h}$ for some $h\in\{1,2\}$, then $|i-j|,|i-l|,|j-l|>1$ and

[TABLE]

(note that this term only appears in $a_{2}$).Finally, if all the indices are different, then assuming without loss of generality that $i>j>h>l$,

[TABLE]

Estimates for $a_{3}$. Arguing as in the proof of Lemma 7.2 in [12] (see equation (7.6) therein), the term $a_{3}$ is quadratic in $\phi_{1}$ and $\phi_{2}$, and so by Proposition 3.1 it satisfies $a_{2}=o(\delta_{1}^{2})$.

Estimates for $a_{4}$. The first term in $a_{4}$ can be estimated as

[TABLE]

because, by (3.40),

[TABLE]

(this term corresponds to the quantity $A$ defined in the proof of Proposition 3.6). As for the second term in $a_{4}$, given $i,j\in\{1,2\}$ with $i\neq j$ and $h,l\in I_{i}$ with $h<l$, $m\in I_{j}$, by Lemma A.6 we have

[TABLE]

since, for instance when $h<l<m$,

[TABLE]

Estimates for $a_{5}$. Starting from the first term, by Proposition 3.1 and Lemma A.5 we see that, given $i,j\in\{1,2\}$ with $i\neq j$, and $h,l\in I_{i}$ with $h<l$,

[TABLE]

Similarly, the second term in $a_{5}$ is also an $\text{o}(\delta_{1}^{2})$.

Estimates for $a_{6}$. We have, by Proposition 3.1 and Lemma (A.3),

[TABLE]

while the second and third terms in $a_{6}$ are respectively of second and fourth order in $\phi$, thus an $\text{o}(\delta_{1}^{2})$. This ends the proof. ∎

As a direct consequence of (4.7), Corollary 4.3, Corollary 4.5 and Lemma 4.6, we have the following result, which gives us the leading term of the expansion of the reduced energy.

Proposition 4.7**.**

We have

[TABLE]

as $\varepsilon\to 0$, uniformly in $\mathbf{d}\in X_{\eta}$.

Remark 4.8**.**

Suppose that, instead of dealing with the set of odd and even numbers of $\{1,\ldots,k\}$ in the two equation case, we are dealing with system (1.4) with $m$ equations and with a general partition $I_{1},\ldots,I_{m}$ satisfying (1)–(5). Then the reduced energy reads as

[TABLE]

with

[TABLE]

With an analogous proof of the one of Lemma 4.1, we can show that critical points of this functional correspond to solutions of (1.4). The choice of rates is still (2.5) in the general case. As pointed out in the proofs of Corollary 4.5 and Lemma 4.6, besides the exact shape of the rates, the other crucial step is that each set $I_{h}$ does not contain two consecutive integers. Since this property is valid for a general partition (it corresponds to property (5)), it is straightforward to adapt the proofs of these results and show that the quantity

[TABLE]

has the asymptotic expansion (4.26), and that $R$ satisfies (4.29). Combining this with Corollary 4.3 yields that, in the general case, the reduced expression has the exact same expansion, namely (4.31).

5. Proof of the main result

In this section we conclude the proof of Theorem 1.3. Define $\Psi:(\mathbb{R}^{+})^{k}\to\mathbb{R}$ as

[TABLE]

where, since $\beta<0$,

[TABLE]

Lemma 5.1**.**

The function $\Psi$ achieves a unique global minimum at $(x_{1}^{*},\ldots,x_{k}^{*})$, with

[TABLE]

for $i=1,\ldots,k$. In particular, the conclusion of Theorem 1.3 holds true.

Proof.

First of all, observe that

[TABLE]

since $a_{1},a_{2},a_{3}>0$ and $a_{1}x_{1}\to\infty$ if $x_{1}\to\infty$, $x_{i+1}/x_{i}\to\infty$ if $x_{i+1}\to\infty$ and $x_{i}$ is bounded. Moreover, given $\bar{\mathbf{x}}\in(\mathbb{R}^{+}_{0})^{k}$ with $\bar{x}_{i}=0$,

[TABLE]

since in this case at least one of the $k+1$ terms in the expression of $\Psi$ divergences to $\infty$. In conclusion, $\Psi$ admits a global minimum. Let us see next that it is unique, and deduce its expression.

We have

[TABLE]

and

[TABLE]

Hence, at a critical point,

[TABLE]

which yields, by direct substitution,

[TABLE]

The last identity gives

[TABLE]

and the rest of the proof follows. ∎

End of the proof of Theorem 1.3.

From the definition of $\Psi$ and by Proposition 4.7, we have

[TABLE]

where $\textrm{o}(1)\to 0$ as $\varepsilon\to 0$, uniformly in $\mathbf{d}\in X_{\eta}$. Let $d_{i}^{*}:=\sqrt{x_{i}^{*}}$ (cf. Lemma 5.1), and take $\eta>0$ small enough so that $(d_{1}^{*},\ldots,d_{k}^{*})\in X_{\eta}$. Let $K\Subset X_{\eta}$ be a compact set such that $((d_{1}^{*})^{2},\ldots,(d_{k}^{*})^{2})\in\text{int}K$ and

[TABLE]

Then

[TABLE]

Therefore $J_{\varepsilon}|_{K}$ has a minimizer $\mathbf{d}_{\varepsilon}$, which converves to $\mathbf{d}^{*}$ (by the uniqueness stated in Lemma 5.1). Thus $J_{\varepsilon}^{\prime}(\mathbf{d}^{*})=0$. By invoking Lemma 4.1, the proof is finished. ∎

Remark 5.2**.**

The proof of the general case, Theorem 1.1, follows exactly in the same way since, as we commented on Remark 4.8, the reduced functional is the same as in the two equation case.

Appendix A Asymptotic estimates

In this appendix we collect several important asymptotic estimates which are used in the paper.We assume in every statement that $N=4$, that is, $\Omega\subset\mathbb{R}^{4}$.

The following two results are taken from [12], see Lemmas 3.1 and 3.2 therein.

Lemma A.1**.**

Let $a\in\Omega$, $r>0$, and $\tau\in\mathbb{R}^{4}$. Assume that $\xi=a+\delta\tau$, with $\delta=\delta(\varepsilon)\to 0$ and $\varepsilon/\delta\to 0$ as $\varepsilon\to 0^{+}$. Then, for $A=\int_{\mathbb{R}^{N}}U_{1,0}^{3}=\frac{\alpha_{4}}{\gamma_{4}}$ and $R$ defined by

[TABLE]

there exists $C=C(\tau,\textrm{dist}(a,\partial\Omega))>0$ such that, for any $x\in\Omega\setminus B_{r\varepsilon}(a)$,

[TABLE]

Lemma A.2**.**

Under the assumptions and notations of the previous lemma, we have the following estimate:

[TABLE]

The following concerns the asymptotic study of $L^{q}$ norms of the bubble. For the proof, see for instance [20, Lemma A.3] or [21, Lemma A.2].

Lemma A.3**.**

We have, as $\delta\to 0$,

[TABLE]

The following lemmas will be used many times in order to estimate interaction integrals.

Lemma A.4**.**

Let $1<q<2<p$ be such that $p+q=4$. Let $\rho_{1}=\rho_{1}(\varepsilon)>0$, $\rho_{2}=\rho_{2}(\varepsilon)>0$, be such that

[TABLE]

as $\varepsilon\to 0^{+}$. Then

[TABLE]

Proof.

We proceed by direct computations:

[TABLE]

as desired. Analogously

[TABLE]

Lemma A.5**.**

Let $\rho_{1}=\rho_{1}(\varepsilon)>0$, $\rho_{2}=\rho_{2}(\varepsilon)>0$, be such that

[TABLE]

as $\varepsilon\to 0^{+}$. Then

[TABLE]

Proof.

We have

[TABLE]

Since $\Omega$ is bounded, there exists a sufficiently large radius $R>0$ such that

[TABLE]

To sum up

[TABLE]

and the thesis follows. ∎

Lemma A.6**.**

Let $\rho_{1}=\rho_{1}(\varepsilon)>0$, $\rho_{2}=\rho_{2}(\varepsilon)>0$, $\rho_{3}=\rho_{3}(\varepsilon)>0$ be such that

[TABLE]

as $\varepsilon\to 0^{+}$. Then

[TABLE]

as $\varepsilon\to 0^{+}$.

Proof.

Again, by direct computations

[TABLE]