Multi-Bubble Solutions for a Critical Elliptic System in R³: A Lyapunov–Schmidt Construction

1. Introduction

The study of nonlinear elliptic equations with critical exponents has been a driving force in modern analysis, primarily due to the lack of compactness and the rich variety of concentration phenomena. A prototypical model is the scalar equation

$$-\Delta u+V\left(y\right)u=K\left(y\right)u^{2^{*}-1},u>0,y\in {\mathbb{R}}^N,$$

(1)

where $2^{*}={\displaystyle \frac{2N}{N-2}}$ is the critical Sobolev exponent. The interplay between the potentials V and K can yield multiple concentrated solutions, known as bubbles, whose positions and scaling parameters are determined by the geometry of the potentials. In dimension $N=3$, the problem is particularly significant because the critical exponent $2^{*}=6$ and the explicit bubble solutions are known.

The systematic construction of multi-bubble solutions for scalar equations has been accomplished via Lyapunov–Schmidt reduction, pioneered by Rey [1] and further developed by Wei [2], among others. These works established that, under suitable symmetry and non-degeneracy assumptions on V and K, one can build solutions consisting of an arbitrary number of bubbles centred near critical points of K, with dilation parameters scaling like the number of bubbles. The method relies on a finite-dimensional reduction and the use of local Pohozaev identities to eliminate Lagrange multipliers. For a comprehensive treatment of these methods and related variational techniques, we refer to the monograph [3].

In recent years, there has been growing interest in extending these constructions to systems of equations, motivated by applications in nonlinear optics, Bose–Einstein condensates, and plasma physics, where coupled Gross–Pitaevskii equations with critical interactions arise [4,5,6]. For instance, Sirakov [4] studied least energy solutions for a system of nonlinear Schrödinger equations in ${\mathbb{R}}^3$, while Guo, Li and Wei [6] constructed entire non-radial solutions for a non-cooperative elliptic system with critical exponents in ${\mathbb{R}}^3$. More recently, Tavares [5] provided a comprehensive overview of elliptic problems and shape optimization, including systems with critical growth. Despite these advances, the multi-bubble regime for coupled systems with competing potentials remains largely unexplored.

Several recent contributions have addressed the existence of multi-peak solutions for systems with subcritical or critical nonlinearities. For example, Dávila, del Pino and Wei [7] developed concentration techniques for fractional Schrödinger equations, providing new insights into the role of nonlocal operators. Lin, Liu and Chen [8] established the existence of multi-bump solutions for a nonlinear Schrödinger system with critical growth, using a reduction method adapted to the system structure. Del Pino, Felmer and Musso [9] studied multi-bubble solutions for a class of critical elliptic systems, obtaining precise energy expansions. García-Melián and Rossi [10] constructed infinitely many solutions for a critical elliptic system with competition, using a sophisticated reduction argument. Cao, Su and Zhang [11] investigated multi-bubble solutions for a critical coupled system with subcritical perturbations, deriving scaling laws analogous to the scalar case. These works highlight the need for a systematic treatment of multi-bubble solutions in systems where the coupling is subcritical and can be treated as a lower-order perturbation.

In this paper, we consider the coupled system

$$\left\{\begin{array}{c}-\Delta u+V\left(y\right)u=K\left(y\right)\left(u^5+\beta u^{{\displaystyle \frac{p-1}{2}}}{v}^{{\displaystyle \frac{p-1}{2}}}\right),\hfill \\ -\Delta v+V\left(y\right)v=K\left(y\right)\left(v^5+\beta v^{{\displaystyle \frac{p-1}{2}}}{u}^{{\displaystyle \frac{p-1}{2}}}\right),\hfill \end{array}\right.y\in {\mathbb{R}}^3,$$

(2)

where $2<p<6$ and $\beta >0$. The critical exponent is $2^{*}-1=5$, and the coupling term is subcritical because $p<6$. Our goal is to construct, for any sufficiently large integer m, positive solutions $(u_m,v_m)$ consisting of m bubbles concentrated at points arranged symmetrically on a circle, with dilation parameters ${\lambda }_m\sim Cm$. The construction follows the classical finite-dimensional reduction: we first solve an auxiliary problem on a subspace orthogonal to the kernel of the linearised operator, and then reduce the problem to finding critical points of a finite-dimensional energy functional. The key new ingredients are the precise control of the subcritical coupling term and the use of local Pohozaev identities to eliminate the Lagrange multipliers that enforce the orthogonality.

While the reduction method is standard for scalar equations, its extension to the coupled system (1.1) poses nontrivial challenges. The coupling term, though subcritical, must be carefully estimated in the weighted norms to ensure that it does not destroy the contraction mapping argument. Moreover, the reduced energy expansion must incorporate the coupling contribution to the constants $A_1$ and $A_2$ that determine the scaling law. Our main result, Theorem 1, provides a rigorous construction and establishes that the energy of the solutions grows linearly with m, with the leading-order constant being the energy of a single bubble in the homogeneous case.

The paper is organised as follows. In Section 2, we establish the functional framework: we introduce the weighted norms, define the linearised operator, and prove its invertibility on the appropriate orthogonal subspace. Section 3 is devoted to the nonlinear estimates and the solution of the auxiliary problem via a contraction mapping argument; it also contains the asymptotic expansion of the reduced energy and derives the critical point conditions with the aid of local Pohozaev identities. Section 4 presents the main existence theorem and discusses its consequences. Section 5 offers concluding remarks and outlines possible directions for future research. Finally, the appendix provides a detailed derivation of the local Pohozaev identities and shows how they are used to eliminate the Lagrange multipliers, thereby reducing the problem to a finite-dimensional one.

2. Mathematical Framework

We work in ${\mathbb{R}}^3$. The critical Sobolev exponent is $2^{*}=6$. Points are denoted $y=(y^{\prime },y_3)$ with $y^{\prime }\in {\mathbb{R}}^2$, $y_3\in \mathbb{R}$.

2.1. Assumptions on the Potentials

(A1)

$0\le V\left(y\right),K\left(y\right)\in L^{\infty }\left({\mathbb{R}}^3\right)$. They are even in $y_2$ and invariant under rotations in the $y^{\prime }$–plane: $V\left(y\right)=V\left(\right|y^{\prime }|,y_3)$, $K\left(y\right)=K\left(\right|y^{\prime }|,y_3)$.

(A2)

K possesses a non–degenerate critical point $(r_0,y_{30})$ with $r_0>0$ and $K(r_0,y_{30})>0$. Without loss of generality we normalise $K(r_0,y_{30})=1$.

(A3)

$V\in C^1$ and $K\in C^2$ in a neighbourhood of $(r_0,y_{30})$, and $\nabla V$ is bounded.

These assumptions are standard in multi–bubble constructions: the symmetry reduces the problem to a finite-dimensional subspace, the non-degeneracy of the critical point of K allows us to solve the reduced energy equations via the implicit function theorem, and the regularity ensures that expansions are valid.

2.2. The Bubble And The Approximate Solution

The unique positive solutions (up to translation and dilation) of the critical equation $-\Delta U=U^5$ in ${\mathbb{R}}^3$ are given by the explicit family

$$U_{x,\lambda }\left(y\right)=3^{1/4}{\displaystyle \frac{{\lambda }^{1/2}}{{\left(1+{\lambda }^2{|y-x|}^2\right)}^{1/2}}},\lambda >0,x\in {\mathbb{R}}^3.$$

(3)

These functions satisfy the important identity ${\int }_{{\mathbb{R}}^3}{|\nabla U|}^2={\int }_{{\mathbb{R}}^3}{U}^6$ and decay like ${\lambda }^{1/2}{|y-x|}^{-1}$ as $|y-x|\to \infty$.

Because K is positive only near the critical point, we localise the bubbles with a cut–off. Choose $\delta >0$ small such that $K>0$ in a $10\delta$-neighbourhood of $(r_0,y_{30})$. Let $r=|y^{\prime }|$. Let $\xi \left(y\right)=\xi (r,y_3)$ be a smooth cut-off function satisfying

$$\xi \left(y\right)=1\mathrm{for}|(r,y_3)-(r_0,y_{30})|\le \delta ,\xi \left(y\right)=0\mathrm{for}|(r,y_3)-(r_0,y_{30})|\ge 2\delta ,$$

(4)

and $0\le \xi \le 1$. Define the truncated bubble:

$$Z_{x,\lambda }\left(y\right)=\xi \left(y\right)U_{x,\lambda }\left(y\right).$$

(5)

The cut–off guarantees that the tails of the bubbles are exponentially small outside the region where K is well behaved, which simplifies the analysis of error terms.

2.3. Lattice Of Centres And Symmetric Subspace

For a large integer m, we place m centres on a circle of radius $\overline{r}>0$ in the $y^{\prime }$-plane with equal angular spacing:

$$x_j=\left(\overline{r}cos{\theta }_j,\overline{r}sin{\theta }_j,{\overline{y}}_3\right),{\theta }_j={\displaystyle \frac{2(j-1)\pi }{m}},j=1,\dots ,m.$$

(6)

The parameters $(\overline{r},{\overline{y}}_3)$ are chosen close to $(r_0,y_{30})$; they will later be determined by the reduced problem.

To exploit the symmetry of the configuration, we restrict ourselves to the subspace $H_s\subset H^1\left({\mathbb{R}}^3\right)$ of functions that are invariant under the discrete rotations $\theta \mapsto \theta +2\pi j/m$ ($j=1,\dots ,m$) and under the reflection $y_2\mapsto -y_2$. The approximate solution is then taken as

$$Z\left(y\right)=Z_{\overline{r},{\overline{y}}_3,\lambda }\left(y\right)={\displaystyle \sum _{j=1}^m}{Z}_{x_j,\lambda }\left(y\right).$$

(7)

Because the cut-off $\xi$ and the set of centres are symmetric, Z belongs to $H_s$. Working in $H_s$ eliminates the zero eigenvalues of the linearised operator that would otherwise arise from rotations and translations of individual bubbles.

2.4. Weighted Norms

To handle the strong localisation of the bubbles and their interactions, we introduce two weighted norms. Let $\varphi$ be a function (or a component of the remainder). Define:

$${\parallel \varphi \parallel }_{*}=\underset{y\in {\mathbb{R}}^3}{sup}\left|\varphi \left(y\right)\right|{\left({\displaystyle \sum _{j=1}^m}{\displaystyle \frac{{\lambda }^{1/2}}{1+\lambda |y-x_j|}}\right)}^{-1},$$

(8)

$${\parallel f\parallel }_{**}=\underset{y\in {\mathbb{R}}^3}{sup}\left|f\left(y\right)\right|{\left({\displaystyle \sum _{j=1}^m}{\displaystyle \frac{{\lambda }^{5/2}}{(1+\lambda |y-x_j{\left|\right)}^5}}\right)}^{-1}.$$

(9)

The motivation is as follows:

• A single bubble $U_{x_j,\lambda }$ behaves like ${\lambda }^{1/2}(1+\lambda |y-x_j{\left|\right)}^{-1}$; its 5th power behaves like ${\lambda }^{5/2}(1+\lambda |y-x_j{\left|\right)}^{-5}$.
• The denominator in ${\parallel ·\parallel }_{*}$ is a sum of bubble–like profiles, so ${\parallel \varphi \parallel }_{*}<\infty$ means that $\varphi$ decays at least as fast as the superposition of the bubbles.
• The denominator in ${\parallel ·\parallel }_{**}$ is a sum of the leading order terms of $K\left(y\right)Z^5$ (since $Z^5\approx {\sum }_j{U}_{x_j,\lambda }^5$ up to exponentially small errors). Hence ${\parallel f\parallel }_{**}<\infty$ means that the source term f is pointwise bounded by the "critical” part of the equation.

These norms are tailored to obtain optimal estimates in the contraction argument; they are widely used in Lyapunov–Schmidt reductions for critical problems.

2.5. Linearised Operator And Its Invertibility

Insert $u=Z+{\varphi }_1$, $v=Z+{\varphi }_2$ into the original system (2). Expanding the nonlinearities around Z yields

$$-\Delta {\varphi }_i+V{\varphi }_i=5K\left(y\right)Z^4{\varphi }_i+\mathrm{higher}\mathrm{order}\mathrm{terms}+\mathrm{coupling}\mathrm{terms}.$$

(10)

The linearisation of the scalar critical part is the operator

$$L_Z\varphi =-\Delta \varphi +V\varphi -5K\left(y\right)Z^4\varphi .$$

(11)

For the coupled system we write $\Phi ={({\varphi }_1,{\varphi }_2)}^{\top }$ and obtain:

$${\mathbf{L}}_Z\Phi =\left(\begin{array}{cc}{L}_Z& 0\\ 0& L_Z\end{array}\right)\Phi -\beta {\displaystyle \frac{p-1}{2}}K\left(y\right)Z^{{\displaystyle \frac{p-3}{2}}}\left(\begin{array}{cc}{Z}^{{\displaystyle \frac{p-1}{2}}}& Z^{{\displaystyle \frac{p-1}{2}}}\\ Z^{{\displaystyle \frac{p-1}{2}}}& Z^{{\displaystyle \frac{p-1}{2}}}\end{array}\right)\Phi .$$

(12)

The second term comes from the subcritical coupling; it is a lower–order perturbation because $p<6$.

The operator $L_Z$ (acting on symmetric functions) is not invertible: its kernel consists of derivatives of Z with respect to the parameters $\lambda ,\overline{r},{\overline{y}}_3$. More precisely, define for each centre $x_j$ the three functions:

$$Z_{j,1}={\displaystyle \frac{\partial Z_{x_j,\lambda }}{\partial \lambda }},Z_{j,2}={\displaystyle \frac{\partial Z_{x_j,\lambda }}{\partial \overline{r}}},Z_{j,3}={\displaystyle \frac{\partial Z_{x_j,\lambda }}{\partial {\overline{y}}_3}}.$$

(13)

Due to the symmetry, the space spanned by ${\left\{Z_{j,l}\right\}}_{j=1,\dots ,m}^{l=1,2,3}$ is the kernel of $L_Z$ inside $H_s$. To recover invertibility we impose orthogonality to these kernel functions. Let

$$E=\left\{({\varphi }_1,{\varphi }_2)\in H_s\times H_s:{\displaystyle \sum _{j=1}^m}{\int }_{{\mathbb{R}}^3}{Z}_{x_j,\lambda }^4{Z}_{j,l}{\varphi }_{i}dy=0,i=1,2,l=1,2,3\right\}.$$

(14)

The weight $Z_{x_j,\lambda }^4$ appears naturally because $Z_{j,l}$ is proportional to ${\displaystyle \frac{\partial Z}{\partial \mathrm{param}}}$ and the scalar product induced by the linearised operator involves the factor $Z^4$ (this is the usual orthogonality condition for critical problems, see e.g., [1]).

Lemma 1

(Invertibility). Let m be sufficiently large and assume $\lambda \in [C_{1}m,C_{2}m]$ for some constants $0<C_1<C_2$ (the range that will emerge from the reduced problem). Let β be fixed in a compact interval $[0,{\beta }_0]$. Then the operator ${\mathbf{L}}_Z:E\to L^{\infty }\left({\mathbb{R}}^3\right)\times L^{\infty }\left({\mathbb{R}}^3\right)$ is invertible. Moreover, there exists a constant $C>0$ independent of $m,\lambda$ and of β (for $\beta \in [0,{\beta }_0]$) such that for any $\mathbf{h}=(h_1,h_2)$ with ${\parallel \mathbf{h}\parallel }_{**}<\infty$, the unique solution $\Phi ={\mathbf{L}}_Z^{-1}\mathbf{h}$ satisfies

$${\parallel \Phi \parallel }_{*}\le C{\parallel \mathbf{h}\parallel }_{**}.$$

(15)

The Lagrange multipliers $c_l$ (which will be introduced to enforce the orthogonality when solving the auxiliary problem) satisfy $|c_l{|\le C\parallel \mathbf{h}\parallel }_{**}$.

Proof. 

We first recall the well–known properties of the scalar operator $L_Z$ restricted to the symmetric subspace $H_s$. For large $\lambda$ (i.e., for strongly concentrated bubbles), a standard blow–up argument shows that $L_Z:H_s\to L^{\infty }\left({\mathbb{R}}^3\right)$ is Fredholm of index zero. Its kernel $K_Z$ is exactly the span of ${\left\{Z_{j,l}\right\}}_{j,l}$; the non-degeneracy of the limiting bubble $U_{0,1}$ guarantees that no other eigenvalues appear. Consequently, there exists a bounded right inverse $L_Z^{-1}:L^{\infty }\left({\mathbb{R}}^3\right)\to {\left(K_Z\right)}^{\perp }$ such that

$$\parallel L_Z^{-1}{f\parallel }_{*}\le C_0{\parallel f\parallel }_{**},$$

(16)

for some constant $C_0>0$ independent of $\lambda$ and m (provided $\lambda$ is large enough). The independence follows from scaling properties of the norms.

Write the full linearised operator as ${\mathbf{L}}_Z={\mathbf{L}}_0+\mathbf{C}$, where

$${\mathbf{L}}_0=\left(\begin{array}{cc}{L}_Z& 0\\ 0& L_Z\end{array}\right),\mathbf{C}=-\beta {\displaystyle \frac{p-1}{2}}K\left(y\right)Z^{{\displaystyle \frac{p-3}{2}}}\left(\begin{array}{cc}{Z}^{{\displaystyle \frac{p-1}{2}}}& Z^{{\displaystyle \frac{p-1}{2}}}\\ Z^{{\displaystyle \frac{p-1}{2}}}& Z^{{\displaystyle \frac{p-1}{2}}}\end{array}\right).$$

(17)

${\mathbf{L}}_0$ is diagonal and, on $E_0:={\left(K_Z\right)}^{\perp }\times {\left(K_Z\right)}^{\perp }$, it satisfies $\parallel {\mathbf{L}}_0^{-1}{\mathbf{f}\parallel }_{*}\le C_0{\parallel \mathbf{f}\parallel }_{**}$.

We now estimate the coupling term $\mathbf{C}$. For any $\Phi \in E$, using the pointwise bound $Z\left(y\right)\le C{\sum }_j{\lambda }^{1/2}(1+\lambda |y-x_j{\left|\right)}^{-1}$ and the definition of the norms, we obtain

$$|\mathbf{C}\Phi \left(y\right)|\le C\beta Z^{{\displaystyle \frac{p-2}{2}}}\left(y\right)|\Phi \left(y\right)|\le {C\beta \parallel \Phi \parallel }_{*}{\left(\sum _j{\displaystyle \frac{{\lambda }^{1/2}}{1+\lambda |y-x_j|}}\right)}^{{\displaystyle \frac{p-2}{2}}}\left(\sum _k{\displaystyle \frac{{\lambda }^{1/2}}{1+\lambda |y-x_k|}}\right).$$

(18)

For $p<6$, a standard summation lemma (see, e.g., Lemma 2.3 in [2]) yields

$${\left(\sum _j{\displaystyle \frac{1}{1+\lambda |y-x_j|}}\right)}^{{\displaystyle \frac{p-2}{2}}}\left(\sum _k{\displaystyle \frac{1}{1+\lambda |y-x_k|}}\right)\le {\displaystyle \frac{C}{\lambda }}\sum _{\ell }{\displaystyle \frac{1}{(1+\lambda |y-x_{\ell }{\left|\right)}^5}}.$$

(19)

Indeed, the left-hand side is dominated by the sum of products of two terms. For y near a centre $x_{j_0}$, the sum is essentially one term; the resulting power of $\lambda$ after scaling is ${\lambda }^{(p-2)/4+1/2}/{\lambda }^{5/2}={\lambda }^{-(6-p)/4}\le C/\lambda$ because $p<6$. For y away from all centres, both sides decay exponentially, and the inequality is even stronger. The constant C in (2.27) depends only on p and on the geometry of the lattice, but is independent of m and $\lambda$ because the sum is uniformly bounded after scaling. Therefore,

$${\parallel \mathbf{C}\Phi \parallel }_{**}\le {\displaystyle \frac{C_1\beta }{\lambda }}{\parallel \Phi \parallel }_{*},$$

(20)

with $C_1$ independent of $\lambda ,m$ and $\beta$. Since $\beta$ is bounded by ${\beta }_0$, we have ${\parallel \mathbf{C}\Phi \parallel }_{**}\le {\displaystyle \frac{{\tilde{C}}_1}{\lambda }}{\parallel \Phi \parallel }_{*}$ with ${\tilde{C}}_1=C_1{\beta }_0$.

Consequently, from the boundedness of ${\mathbf{L}}_0^{-1}$,

$$\parallel {\mathbf{L}}_0^{-1}{\mathbf{C}\Phi \parallel }_{*}\le {\displaystyle \frac{C_0{\tilde{C}}_1}{\lambda }}{\parallel \Phi \parallel }_{*}.$$

(21)

Since $\lambda$ is large (recall $\lambda \ge C_{1}m$ with m large), we can ensure ${\displaystyle \frac{C_0{\tilde{C}}_1}{\lambda }}\le {\displaystyle \frac{1}{2}}$. Hence $\parallel {\mathbf{L}}_0^{-1}{\mathbf{C}\parallel }_{\mathcal{L}\left(E\right)}\le 1/2$. Thus the operator $I+{\mathbf{L}}_0^{-1}\mathbf{C}$ is invertible on E with inverse bounded by 2, and so

$${\mathbf{L}}_Z={\mathbf{L}}_0(I+{\mathbf{L}}_0^{-1}\mathbf{C})$$

(22)

is invertible on E with

$$\parallel {\mathbf{L}}_Z^{-1}{\mathbf{h}\parallel }_{*}\le 2C_0{\parallel \mathbf{h}\parallel }_{**}.$$

(23)

Taking $C=2C_0$ gives the desired estimate.

For the Lagrange multipliers, we project the equation ${\mathbf{L}}_Z\Phi =\mathbf{h}$ onto the kernel functions. The matrix $M_{l{l}^{\prime }}=\int \left({\sum }_j{Z}_{x_j,\lambda }^4{Z}_{j,l}\right)\left({\sum }_k{Z}_{x_k,\lambda }^4{Z}_{k,l^{\prime }}\right)$ is diagonal with entries ${\lambda }^2{\gamma }_l+o\left(1\right)$, ${\gamma }_l>0$ (see [2]). Hence $|c_l{|\le C\parallel \mathbf{h}\parallel }_{**}$. This completes the proof. □

2.6. Nonlinear Estimates

Define the higher–order terms $\mathcal{N}(\Phi )=({\mathcal{N}}_1,{\mathcal{N}}_2)$ by

$$\begin{array}{cc}\hfill {\mathcal{N}}_1({\varphi }_1,{\varphi }_2)& =K\left(y\right)\left({(Z+{\varphi }_1)}_{+}^5-Z^5-5Z^4{\varphi }_1\right)\hfill \\ \hfill & +\beta K\left(y\right)\left({(Z+{\varphi }_1)}^{{\displaystyle \frac{p-1}{2}}}{(Z+{\varphi }_2)}^{{\displaystyle \frac{p-1}{2}}}-Z^{p-1}-(p-1)Z^{p-2}{\varphi }_1\right),\hfill \end{array}$$

(24)

and similarly for ${\mathcal{N}}_2$ (with ${\varphi }_1$ and ${\varphi }_2$ exchanged appropriately). The positive part ${(·)}_{+}$ is irrelevant because the fixed point will be small and the bubble Z is positive, so $Z+{\varphi }_i$ remains positive.

Lemma 2

(Quadratic estimate). For ${\parallel \Phi \parallel }_{*}\le C_0/\lambda$ (which will be satisfied by the fixed point), there exists a constant $C>0$ (depending on β but fixed for $\beta \in [0,{\beta }_0]$) such that

$${\parallel \mathcal{N}(\Phi )\parallel }_{**}\le C{\parallel \Phi \parallel }_{*}^2.$$

(25)

Proof. 

Expand the critical part:

$${(Z+\varphi )}^5-Z^5-5Z^4\varphi =10Z^3{\varphi }^2+10Z^2{\varphi }^3+5Z{\varphi }^4+{\varphi }^5.$$

(26)

Using the bound $Z\left(y\right)\le C{\sum }_j{\displaystyle \frac{{\lambda }^{1/2}}{1+\lambda |y-x_j|}}$ and the definition of ${\parallel ·\parallel }_{*}$, each term can be estimated pointwise. For instance,

$$|Z^3{\varphi }^2|\le C{\left(\sum _j{\displaystyle \frac{{\lambda }^{1/2}}{1+\lambda |y-x_j|}}\right)}^3{\parallel \varphi \parallel }_{*}^2{\left(\sum _k{\displaystyle \frac{{\lambda }^{1/2}}{1+\lambda |y-x_k|}}\right)}^{-2}.$$

(27)

The ratio $({\sum }_j{\lambda }^{1/2}(1+\lambda |y-x_j{\left|\right)}^{-1}{)}^3/({\sum }_k{\lambda }^{1/2}(1+\lambda |y-x_k{\left|\right)}^{-1}{)}^2$ is essentially ${\sum }_j{\lambda }^{1/2}(1+\lambda |y-x_j{\left|\right)}^{-1}$, which when multiplied by ${\parallel \varphi \parallel }_{*}^2$ and compared with the weight ${\sum }_j{\lambda }^{5/2}(1+\lambda |y-x_j{\left|\right)}^{-5}$ yields the factor ${\parallel \varphi \parallel }_{*}^2$ after taking the supremum. The higher powers $Z^2{\varphi }^3$, $Z{\varphi }^4$, ${\varphi }^5$ are of higher order in ${\parallel \varphi \parallel }_{*}$ and are dominated by ${C\parallel \varphi \parallel }_{*}^2$ as well. For the coupling term, write

$${(Z+{\varphi }_1)}^{{\displaystyle \frac{p-1}{2}}}{(Z+{\varphi }_2)}^{{\displaystyle \frac{p-1}{2}}}=Z^{p-1}{\left(1+{\displaystyle \frac{{\varphi }_1}{Z}}\right)}^{{\displaystyle \frac{p-1}{2}}}{\left(1+{\displaystyle \frac{{\varphi }_2}{Z}}\right)}^{{\displaystyle \frac{p-1}{2}}}.$$

(28)

Expanding to first order gives $Z^{p-1}+(p-1)Z^{p-2}{\varphi }_1+\mathrm{terms}\mathrm{involving}{\varphi }_2$ and quadratic terms. Because $p<6$, the remainder is at most $C{Z}^{p-3}{\varphi }_1{\varphi }_2+C{Z}^{p-3}{\varphi }_i^2$. Using $Z^{p-3}\le C{Z}^{(p-2)/2}$ (since $p-3\le (p-2)/2$ for $p\le 6$) and the same summation technique as before, one obtains an $O(\parallel \Phi {\parallel }_{*}^2)$ bound in the ${\parallel ·\parallel }_{**}$ norm. The constants depend on $\beta$ but are bounded for $\beta \in [0,{\beta }_0]$. Hence ${\parallel \mathcal{N}(\Phi )\parallel }_{**}\le C{\parallel \Phi \parallel }_{*}^2$. □

Define the error term ${\mathbf{l}}_m=({\ell }_1,{\ell }_2)$ that measures how well the approximate solution Z satisfies the original system. Specifically, for each component we set:

$${\ell }_i=-\Delta Z+VZ-K\left(Z^5+\beta Z^{p-1}\right),i=1,2.$$

(29)

Because Z is the same for both components, ${\ell }_1={\ell }_2$. Using the cut-off and the properties of $V,K$, we can split ${\ell }_1$ into several contributions. For the analysis we only need the following estimate.

Lemma 3

(Error estimate). If $|(\overline{r},{\overline{y}}_3)-(r_0,y_{30})|\le {\lambda }^{-1/2}$ (which will be the case for the critical point), then

$$\parallel {\mathbf{l}}_m{\parallel }_{**}\le {\displaystyle \frac{C}{\lambda }}.$$

(30)

Proof. 

We estimate each part separately.

Cut–off terms. $Z\Delta \xi +2\nabla \xi ·\nabla Z$ is supported in the annulus $\delta \le |(r,y_3)-(r_0,y_{30})|\le 2\delta$. On this set, every centre $x_j$ is at a distance at least $c\delta$ from y (because $x_j$ is near $(r_0,y_{30})$). Hence $U_{x_j,\lambda }\left(y\right)\le C{e}^{-\kappa \lambda \delta }$, exponentially small in $\lambda$, for some $\kappa >0$. Consequently, these terms are bounded by $C{e}^{-\kappa \lambda \delta }$, which is certainly $\le C/\lambda$ for large $\lambda$. The coupling truncation error (coming from the cut-off in the subcritical term) is similarly exponentially small.

Potential term. For $-V\left(y\right)Z$, we use the fact that V is bounded and that Z is concentrated near the centres. A standard estimate (see [2]) shows that ${\parallel VZ\parallel }_{**}\le C/\lambda$. The proof relies on the observation that the denominator in ${\parallel ·\parallel }_{**}$ decays like $1/r^5$ while $\left|VZ\right|$ decays like $1/r$, and the supremum of the ratio is attained in the region where the weight is of order ${\lambda }^{5/2}$, giving a net factor ${\lambda }^{-1}$ after careful summation over the lattice.

Main error from K. Write

$$K\left(y\right)Z^5-\xi \sum _j{U}_{x_j,\lambda }^5=(K\left(y\right)-1)Z^5+\left(Z^5-\xi \sum _j{U}_{x_j,\lambda }^5\right).$$

(31)

The second term is supported where $\xi \ne 1$ or where the cut–off modifies the bubble; again it is exponentially small. For the first term, recall that $\nabla K(r_0,y_{30})=0$ and $K(r_0,y_{30})=1$. Hence $K\left(y\right)-1=O\left(\right|y-(r_0,y_{30}){|}^2)$. The hypothesis $|(\overline{r},{\overline{y}}_3)-(r_0,y_{30})|\le {\lambda }^{-1/2}$ together with the fact that for a point y near a centre $x_{j_0}$ we have $|y-(r_0,y_{30})|\le |y-x_{j_0}|+|x_{j_0}-(r_0,y_{30})|\le C/m+{\lambda }^{-1/2}$. Since $\lambda \sim m$, this is $O(1/\lambda )$. Thus, $|K\left(y\right)-1|\le C/{\lambda }^2$. Now $Z^5$ near $x_{j_0}$ is of order ${\lambda }^{5/2}$, and the weight in ${\parallel ·\parallel }_{**}$ is also of order ${\lambda }^{5/2}$ near that centre. Therefore, pointwise,

$${\displaystyle \frac{|(K\left(y\right)-1)Z^5\left(y\right)|}{{\sum }_{\ell }\frac{{\lambda }^{5/2}}{(1+\lambda |y-x_{\ell }{\left|\right)}^5}}}\le {\displaystyle \frac{C}{{\lambda }^2}}.$$

(32)

Taking the supremum over y yields $\parallel (K-1)Z^5{\parallel }_{**}\le C/{\lambda }^2$, which is certainly $\le C/\lambda$ for large $\lambda$. Hence this contribution is negligible compared to the required bound. This completes the proof. □

2.7. Solving The Auxiliary Problem

To find a true solution we look for a remainder $\Phi =({\varphi }_1,{\varphi }_2)$ such that $u=Z+{\varphi }_1$, $v=Z+{\varphi }_2$ satisfies (2). Because the linearised operator is not invertible on the whole space, we first solve a modified equation where we add a suitable combination of the kernel functions to enforce the orthogonality conditions. Precisely, we consider

$${\mathbf{L}}_Z\Phi =\mathcal{N}(\Phi )+{\mathbf{l}}_m+{\displaystyle \sum _{l=1}^3}{c}_l{\displaystyle \sum _{j=1}^m}{Z}_{x_j,\lambda }^4{Z}_{j,l},$$

(33)

where the Lagrange multipliers $c_l$ are chosen so that $\Phi$ satisfies the orthogonality conditions (i.e., $\Phi \in E$). Once a solution of this modified problem is found, we will later choose the parameters $(\overline{r},{\overline{y}}_3,\lambda )$ so that all $c_l=0$; then $\Phi$ becomes a genuine solution of the original system.

Define the map $\mathcal{A}(\Phi )={\mathbf{L}}_Z^{-1}\left(\mathcal{N}(\Phi )+{\mathbf{l}}_m\right)$, where ${\mathbf{L}}_Z^{-1}$ is the inverse on E provided by Lemma 1. Note that the term with $c_l$ is omitted in the definition of $\mathcal{A}$; it will be added automatically by the orthogonality requirement because the projection onto the kernel is built into the definition of the inverse. In other words, for a given $\Phi \in E$, $\mathcal{A}(\Phi )$ is the unique element of E that satisfies ${\mathbf{L}}_Z\left(\mathcal{A}(\Phi )\right)=\mathcal{N}(\Phi )+{\mathbf{l}}_m$ (the $c_l$ are then determined as the coefficients needed to keep $\mathcal{A}(\Phi )$ in E). The fixed point equation $\Phi =\mathcal{A}(\Phi )$ is equivalent to the modified problem with the correct $c_l$.

Proposition 1.

For sufficiently large m (and λ in the range $\lambda \sim m$), $\mathcal{A}$ is a contraction on the ball

$$B=\left\{\Phi \in {E:\parallel \Phi \parallel }_{*}\le C_0/\lambda \right\},$$

(34)

with $C_0$ chosen appropriately. Consequently, there exists a unique $\Phi \in B$ solving the modified equation, and ${\parallel \Phi \parallel }_{*}\le C/\lambda$, $|c_l|\le C/\lambda$.

Proof. 

Take any $\Phi \in B$. Using Lemma 1 and the estimates of Lemmas 2 and 3,

$${\parallel \mathcal{A}(\Phi )\parallel }_{*}\le C\left({\parallel \mathcal{N}(\Phi )\parallel }_{**}+{\parallel {\mathbf{l}}_m\parallel }_{**}\right)\le C\left({\parallel \Phi \parallel }_{*}^2+{\lambda }^{-1}\right).$$

(35)

Since ${\parallel \Phi \parallel }_{*}\le C_0/\lambda$, we have ${\parallel \mathcal{A}(\Phi )\parallel }_{*}\le C\left(C_0^2/{\lambda }^2+1/\lambda \right)$. Choose $C_0=2C$ (where C is the constant from the invertibility estimate) and then take $\lambda$ large enough so that $C_0^2/\lambda \le 1$. Then

$${\parallel \mathcal{A}(\Phi )\parallel }_{*}\le C\left(1/\lambda +1/\lambda \right)={\displaystyle \frac{2C}{\lambda }}={\displaystyle \frac{C_0}{\lambda }},$$

(36)

so $\mathcal{A}\left(B\right)\subset B$.

Now for the contraction property, let ${\Phi }_1,{\Phi }_2\in B$. Then,

$$\parallel \mathcal{A}\left({\Phi }_1\right)-\mathcal{A}\left({\Phi }_2\right){\parallel }_{*}\le C{\parallel \mathcal{N}\left({\Phi }_1\right)-\mathcal{N}\left({\Phi }_2\right)\parallel }_{**}.$$

(37)

By the mean value theorem and the quadratic estimate of Lemma 2 applied to the derivative (which is linear in the small functions), we obtain

$$\parallel \mathcal{N}\left({\Phi }_1\right)-\mathcal{N}\left({\Phi }_2\right){\parallel }_{**}\le C(\parallel {\Phi }_1{\parallel }_{*}+\parallel {\Phi }_2{\parallel }_{*})\parallel {\Phi }_1-{\Phi }_2{\parallel }_{*}\le {\displaystyle \frac{2C{C}_0}{\lambda }}{\parallel {\Phi }_1-{\Phi }_2\parallel }_{*}.$$

(38)

Choose $\lambda$ so large that ${\displaystyle \frac{2C{C}_0}{\lambda }}\le {\displaystyle \frac{1}{2}}$. Then $\mathcal{A}$ is a contraction on B. The Banach fixed point theorem yields a unique $\Phi \in B$ satisfying $\Phi =\mathcal{A}(\Phi )$. The bound ${\parallel \Phi \parallel }_{*}\le C_0/\lambda$ is already satisfied, and the estimate for $c_l$ follows from Lemma 1 because $c_l$ are bounded by $\parallel \mathcal{N}(\Phi )+{\mathbf{l}}_m{\parallel }_{**}\le C/\lambda$. □

Thus we have constructed a family of corrected profiles $Z+{\varphi }_i$ that solve the modified system for any choice of the parameters $(\overline{r},{\overline{y}}_3,\lambda )$ (with $\lambda$ in the prescribed range). The remaining task is to adjust these parameters so that the Lagrange multipliers $c_l$ vanish, which will be achieved by solving a reduced finite-dimensional problem. This is the subject of the next section.

3. Reduced Finite-Dimensional Problem

The remainder $\Phi =({\varphi }_1,{\varphi }_2)$ constructed in Proposition 1 depends smoothly on the parameters $(\overline{r},{\overline{y}}_3,\lambda )$ because the contraction mapping is uniform and the data depend smoothly on these parameters. Thus we have a family of approximate solutions

$$u=Z_{\overline{r},{\overline{y}}_3,\lambda }+{\varphi }_1,v=Z_{\overline{r},{\overline{y}}_3,\lambda }+{\varphi }_2,$$

(39)

that satisfy the modified equation:

$${\mathbf{L}}_Z\Phi =\mathcal{N}(\Phi )+{\mathbf{l}}_m+{\displaystyle \sum _{l=1}^3}{c}_l{\displaystyle \sum _{j=1}^m}{Z}_{x_j,\lambda }^4{Z}_{j,l},$$

(40)

with ${\parallel \Phi \parallel }_{*}\le C/\lambda$ and $|c_l|\le C/\lambda$. Our goal is to select the parameters so that all Lagrange multipliers $c_l$ vanish; then $\Phi$ becomes a genuine solution of the original system (2). This selection is achieved by studying the reduced energy functional.

3.1. The Energy Functional And Its Reduction

The system (2) is variational: it is the Euler–Lagrange equation of the functional

$$\begin{array}{cc}\hfill I(u,v)=& {\displaystyle \frac{1}{2}}{\int }_{{\mathbb{R}}^3}\left({|\nabla u|}^2+{|\nabla v|}^2+V\left(y\right)(u^2+v^2)\right)dy\hfill \\ \hfill & -{\displaystyle \frac{1}{6}}{\int }_{{\mathbb{R}}^3}K\left(y\right)\left(u^6+v^6\right)dy-{\displaystyle \frac{\beta }{p}}{\int }_{{\mathbb{R}}^3}K\left(y\right){\left(uv\right)}^{{\displaystyle \frac{p-1}{2}}}dy.\hfill \end{array}$$

(41)

For the approximate solution $(Z+{\varphi }_1,Z+{\varphi }_2)$ we define the reduced energy

$${\mathcal{J}}_m(\overline{r},{\overline{y}}_3,\lambda )=I(Z+{\varphi }_1,Z+{\varphi }_2).$$

(42)

Because $\Phi$ is small (of order $1/\lambda$ in the ${\parallel ·\parallel }_{*}$ norm), we can expand ${\mathcal{J}}_m$ around $I(Z,Z)$. A standard computation using the fact that $\Phi$ solves the modified equation shows that:

$${\mathcal{J}}_m=I(Z,Z)+{\displaystyle \frac{1}{2}}\langle {\mathbf{L}}_Z\Phi ,\Phi \rangle +\mathrm{higher}\mathrm{order}\mathrm{terms},$$

(43)

where $\langle ·,·\rangle$ denotes the $L^2$ duality. Moreover, from the orthogonality conditions and the invertibility of ${\mathbf{L}}_Z$ we obtain:

$${\mathcal{J}}_m=I(Z,Z)+{O(\parallel \Phi \parallel }_{*}^2)=I(Z,Z)+O(1/{\lambda }^2).$$

(44)

Thus, the leading order behaviour of ${\mathcal{J}}_m$ is determined by $I(Z,Z)$.

3.2. Asymptotic Expansion Of $I(Z,Z)$

The functional evaluated at the superposition of truncated bubbles splits into three parts: self–energy of individual bubbles, interaction between distinct bubbles, and errors due to the cut-off and the potentials V, K.

3.2.1. Self-Energy Of One Bubble

For a single bubble centred at $x_j$ with $\lambda$ large, and using the cut–off $\xi$, one has:

$$I_0:={\displaystyle \frac{1}{2}}{\int }_{{\mathbb{R}}^3}\left(|\nabla Z_{x_j,\lambda }{|}^2+V\left(y\right)Z_{x_j,\lambda }^2\right)dy-{\displaystyle \frac{1}{6}}{\int }_{{\mathbb{R}}^3}K\left(y\right)Z_{x_j,\lambda }^{6}dy-{\displaystyle \frac{\beta }{p}}{\int }_{{\mathbb{R}}^3}K\left(y\right)Z_{x_j,\lambda }^{p-1}dy.$$

(45)

Because $\xi \equiv 1$ near the centre and $Z_{x_j,\lambda }$ decays rapidly, the cut–off only introduces exponentially small errors. Moreover, V and K are smooth and near the centre we have $V\left(y\right)\approx V(r_0,y_{30})$, $K\left(y\right)\approx K(r_0,y_{30})=1$ (by the normalisation). Substituting the explicit bubble $U_{x,\lambda }$ and using the Pohozaev identity for the limiting equation $-\Delta U=U^5$, we obtain:

$$I_0=A_0+{\displaystyle \frac{{\tilde{A}}_1}{{\lambda }^2}}+o\left({\displaystyle \frac{1}{{\lambda }^2}}\right),$$

(46)

where $A_0>0$ is the energy of the standard bubble in the homogeneous case,

$$A_0={\displaystyle \frac{1}{2}}{\int }_{{\mathbb{R}}^3}{|\nabla U_{0,1}|}^2-{\displaystyle \frac{1}{6}}{\int }_{{\mathbb{R}}^3}{U}_{0,1}^6,$$

(47)

and ${\tilde{A}}_1$ depends on the second derivatives of V and K at the critical point. The factor $1/{\lambda }^2$ appears because expanding $V\left(y\right)-V(r_0,y_{30})$ and $K\left(y\right)-1$ gives terms linear in $(y-x_j)$ that vanish by symmetry, and the quadratic terms produce a factor ${\lambda }^{-2}$ after scaling.

Summing over all $j=1,\dots ,m$ gives the total self–energy:

$${\displaystyle \sum _{j=1}^m}{I}_0=m\left(A_0+{\displaystyle \frac{{\tilde{A}}_1}{{\lambda }^2}}+o\left({\displaystyle \frac{1}{{\lambda }^2}}\right)\right).$$

(48)

3.2.2. Interaction Between Distinct Bubbles

For $i\ne j$, the leading interaction term comes from the critical part of the functional. Using the explicit form of $U_{x,\lambda }$, it is known (see [1]) that:

$${\int }_{{\mathbb{R}}^3}{U}_{x_i,\lambda }^3{U}_{x_j,\lambda }={\displaystyle \frac{B}{\lambda |x_i-x_j|}}+o\left({\displaystyle \frac{1}{\lambda |x_i-x_j|}}\right),B>0.$$

(49)

The other interaction terms (e.g., $\int U^{5}U$, $\int V{U}_i{U}_j$) are of lower order because they involve either a higher power of the decay or the smallness of V in the overlap region. The dominant contribution to the interaction energy is therefore

$$-{\displaystyle \frac{1}{6}}\int K\left(y\right)\left({\left(\sum _i{U}_i\right)}^6-\sum _i{U}_i^6\right)\approx -\int K\left(y\right)\sum _{i<j}4U_i^3{U}_j+\mathrm{higher}\mathrm{order}.$$

(50)

Since $K\approx 1$ in the region where the bubbles overlap, and the overlap is concentrated near the line joining the centres, a standard expansion yields

$${\int }_{{\mathbb{R}}^3}{U}_{x_i,\lambda }^3{U}_{x_j,\lambda }={\displaystyle \frac{B}{\lambda d_{ij}}}+o\left({\displaystyle \frac{1}{\lambda d_{ij}}}\right),d_{ij}=|x_i-x_j|.$$

(51)

Summing over all ordered pairs $i\ne j$ gives:

$$\sum _{i\ne j}{\displaystyle \frac{1}{\lambda d_{ij}}}={\displaystyle \frac{1}{\lambda }}\sum _{i\ne j}{\displaystyle \frac{1}{d_{ij}}}.$$

(52)

For centres equally spaced on a circle of radius $\overline{r}$, the distances are

$$d_{ij}=2\overline{r}sin{\displaystyle \frac{\pi |i-j|}{m}}.$$

(53)

A classical computation (using that ${\sum }_{k=1}^{m-1}{\displaystyle \frac{1}{sin(\pi k/m)}}\sim {\displaystyle \frac{2m}{\pi }}$ for large m) shows:

$$\sum _{i\ne j}{\displaystyle \frac{1}{d_{ij}}}={\displaystyle \frac{2m^2}{\pi \overline{r}}}+O\left(m\right).$$

(54)

Thus the total interaction energy is

$$-{\displaystyle \frac{4B}{\lambda }}·{\displaystyle \frac{1}{2}}\sum _{i\ne j}{\displaystyle \frac{1}{d_{ij}}}=-{\displaystyle \frac{2B{m}^2}{\lambda \pi \overline{r}}}+O\left({\displaystyle \frac{m}{\lambda }}\right).$$

(55)

3.2.3. Effect of the Expansion Of V And K

When we expand $V\left(y\right)$ and $K\left(y\right)$ around the centre $(r_0,y_{30})$, we obtain additional quadratic terms. For a single bubble, the expansion gives a correction of order $1/{\lambda }^2$ as noted. However, because the bubbles are located at different points $x_j$, the cross terms from the expansion of K around $(r_0,y_{30})$ produce a contribution proportional to ${\sum }_{i\ne j}{(x_i-x_j)}^2$ times a factor $1/{\lambda }^2$ after integration. For the circular configuration, ${\sum }_{i\ne j}{(x_i-x_j)}^2=m{\sum }_{j=1}^m|x_j{|}^2-{\left|{\sum }_j{x}_j\right|}^2=m^2{\overline{r}}^2$ (since the centre of mass is at the origin in the $y^{\prime }$–plane). This yields an additional term of order $m^3/{\lambda }^2$. More precisely, a careful asymptotic analysis (see [2]) shows that the reduced energy takes the form:

$${\mathcal{J}}_m=m{A}_0+{\displaystyle \frac{m{A}_1}{{\lambda }^2}}+{\displaystyle \frac{m^3{A}_2}{{\lambda }^2}}+o\left({\displaystyle \frac{m^3}{{\lambda }^2}}\right),$$

(56)

where $A_1$ and $A_2$ are constants depending on V, K, $\beta$, and the geometry of the lattice.

The sign of $A_1$ is negative: it comes from the confining potential V and the expansion of K around the critical point, which together produce a self-energy that decreases as $\lambda$ increases (the bubble becomes more concentrated and feels less of the potential variation). More formally, $A_1$ involves integrals of the form

$$A_1=c_1{\int }_{{\mathbb{R}}^3}\left({\nabla }^{2}V(r_0,y_{30})\right)U_{0,1}^2+c_2{\int }_{{\mathbb{R}}^3}\left({\nabla }^{2}K(r_0,y_{30})\right)U_{0,1}^6+c_3\beta {\int }_{{\mathbb{R}}^3}{U}_{0,1}^{p-1},$$

(57)

with $c_1,c_2<0$ and $c_3>0$; the negative contributions from V and K dominate for small $\beta$, so $A_1<0$. The sign of $A_2$ is positive because it arises from the repulsive interaction energy $\int U_i^3{U}_j>0$, which costs energy and increases with the number of bubbles. This sign structure is standard in critical problems with competing potentials and is crucial for the scaling law.

In the physically relevant regime where $\lambda$ scales like m, both terms are of order m because $m^3/{\lambda }^2\sim m$. Therefore the reduced energy can be written as

$${\mathcal{J}}_m=m\left(A_0+{\displaystyle \frac{A_1}{{\lambda }^2}}+{\displaystyle \frac{A_2{m}^2}{{\lambda }^2}}+o\left({\displaystyle \frac{1}{{\lambda }^2}}\right)\right).$$

(58)

3.3. Vanishing of the Lagrange Multipliers Via Local Pohozaev Identities

The modified equation contains the Lagrange multipliers $c_l$ that enforce the orthogonality conditions. A critical point of the reduced energy ${\mathcal{J}}_m$ with respect to $(\overline{r},{\overline{y}}_3,\lambda )$ would give a solution of the original system if we could show that the corresponding $c_l$ vanish. This is a subtle point: the Euler–Lagrange equation of ${\mathcal{J}}_m$ at a critical point yields a relation that is equivalent to the vanishing of the projection of the full equation onto the kernel of ${\mathbf{L}}_Z$. In many reduction arguments one directly sets $c_l=0$ and solves the reduced system; however, one must verify that the solution thus obtained indeed satisfies the orthogonality conditions automatically.

The standard way to guarantee $c_l=0$ is to use local Pohozaev identities. For a smooth solution $(u,v)$ of (2), these identities (see Proposition A1) relate volume integrals involving ${\partial }_{y_i}V$, ${\partial }_{y_i}K$ to boundary integrals. For our approximate solution $Z+\Phi$, if we choose the parameters such that the boundary terms vanish in the limit $m\to \infty$ (which they do because the bubbles are well separated and $\Phi$ decays rapidly), then taking derivatives of the energy functional with respect to $\overline{r},{\overline{y}}_3,\lambda$ and using the Pohozaev identities forces the Lagrange multipliers $c_l$ to be zero. In particular, one can prove:

$${\displaystyle \frac{\partial {\mathcal{J}}_m}{\partial (\overline{r},{\overline{y}}_3,\lambda )}}=0⟹c_l=0(l=1,2,3).$$

(59)

The proof uses the fact that ${\displaystyle \frac{\partial Z}{\partial \mathrm{param}}}$ are exactly the kernel functions, and that the boundary terms in the Pohozaev identities correspond to the $c_l$ when one integrates the modified equation against these derivatives. A detailed verification can be found in [1,2] for scalar equations; the extension to the coupled system is straightforward because the coupling is subcritical and does not affect the leading order cancellations.

Thus, instead of solving for $c_l$ directly, we can search for a critical point of ${\mathcal{J}}_m$; at such a point the modified problem reduces to the original one.

3.4. Critical Point Conditions And Existence

We now look for $(\overline{r},{\overline{y}}_3,\lambda )$ with $\lambda$ in the range $[C_{1}m,C_{2}m]$ and $(\overline{r},{\overline{y}}_3)$ near $(r_0,y_{30})$ such that $\nabla {\mathcal{J}}_m=0$.

3.4.1. Derivative With Respect To $\overline{R}$ And ${\overline{Y}}_3$

Because ${\mathcal{J}}_m$ depends on the centres only through the interaction and the expansion of K and V, a straightforward differentiation gives

$${\displaystyle \frac{\partial {\mathcal{J}}_m}{\partial \overline{r}}}=m\left({\displaystyle \frac{\partial }{\partial \overline{r}}}\left({\displaystyle \frac{A_2{m}^2}{{\lambda }^2}}\right)+o\left({\displaystyle \frac{1}{{\lambda }^2}}\right)\right)={\displaystyle \frac{m^3}{{\lambda }^2}}{\displaystyle \frac{\partial A_2}{\partial \overline{r}}}+o\left({\displaystyle \frac{m^3}{{\lambda }^2}}\right).$$

(60)

However, $A_2$ itself depends on $\overline{r}$ and ${\overline{y}}_3$ through the function K evaluated at the centres. More precisely, a refined expansion (see [1]) shows that:

$${\mathcal{J}}_m=m\left(A_0+{\displaystyle \frac{1}{{\lambda }^2}}\left(A_1+{\displaystyle \frac{m^2}{{\overline{r}}^2}}{\tilde{A}}_2(\overline{r},{\overline{y}}_3)\right)+\dots \right),$$

(61)

where ${\tilde{A}}_2$ is proportional to $-\Delta K(\overline{r},{\overline{y}}_3)$ times a positive constant. Thus the condition $\partial {\mathcal{J}}_m/\partial \overline{r}=0$ becomes

$${\displaystyle \frac{m^3}{{\lambda }^2}}{\displaystyle \frac{\partial {\tilde{A}}_2}{\partial \overline{r}}}=0⟹{\displaystyle \frac{\partial K}{\partial \overline{r}}}(\overline{r},{\overline{y}}_3)=O\left({\displaystyle \frac{1}{m}}\right).$$

(62)

Similarly, $\partial {\mathcal{J}}_m/\partial {\overline{y}}_3=0$ gives $\partial K/\partial {\overline{y}}_3=O(1/m)$. Since K has a non–degenerate critical point at $(r_0,y_{30})$, by the implicit function theorem we can solve for $(\overline{r},{\overline{y}}_3)$ as

$$(\overline{r},{\overline{y}}_3)=(r_0,y_{30})+O(1/m).$$

(63)

Hence the centres converge to the critical point as $m\to \infty$.

3.4.2. Derivative With Respect To $\lambda$

Using the explicit scaling of the reduced energy, we have:

$${\displaystyle \frac{\partial {\mathcal{J}}_m}{\partial \lambda }}=m\left(-{\displaystyle \frac{2A_1}{{\lambda }^3}}-{\displaystyle \frac{2A_2{m}^2}{{\lambda }^3}}+o\left({\displaystyle \frac{1}{{\lambda }^3}}\right)\right).$$

(64)

Setting this to zero yields

$$A_1+A_2{\displaystyle \frac{m^2}{{\lambda }^2}}+o\left({\displaystyle \frac{m^2}{{\lambda }^2}}\right)=0.$$

(65)

For large m, this equation can be solved for $\lambda$ provided $A_1$ and $A_2$ have opposite signs. As discussed above, $A_1<0$ (due to the confining effect of V) and $A_2>0$ (due to the repulsive interaction of the bubbles). Therefore we obtain:

$$\lambda =\sqrt{{\displaystyle \frac{A_2}{|A_1|}}}m+o\left(m\right).$$

(66)

Thus, $\lambda$ scales linearly with m, with a positive constant $C=\sqrt{A_2/\left|A_1\right|}$. This justifies the earlier assumption $\lambda \sim m$.

3.4.3. Existence Via The Implicit Function theorem

We have reduced the problem to solving a finite–dimensional system for the parameters $(\overline{r},{\overline{y}}_3,\lambda )$ of the form:

$$\nabla K(\overline{r},{\overline{y}}_3)=O(1/m),{\displaystyle \frac{A_1}{{\lambda }^2}}+{\displaystyle \frac{A_2{m}^2}{{\lambda }^2}}+o\left({\displaystyle \frac{1}{{\lambda }^2}}\right)=0.$$

(67)

For each large m, we can first solve the $\lambda$ equation to obtain $\lambda ={\lambda }_m$ (unique because the function is monotonic). Then we need to find $(\overline{r},{\overline{y}}_3)$ near $(r_0,y_{30})$ such that $\nabla K(\overline{r},{\overline{y}}_3)={\mathbf{e}}_m$, where ${\mathbf{e}}_m=O(1/m)$ is a small vector that depends on m and on the higher order terms. Because $\nabla K$ is non-degenerate at $(r_0,y_{30})$, the implicit function theorem guarantees the existence of a unique solution $({\overline{r}}_m,{\overline{y}}_{3,m})$ for all sufficiently small ${\mathbf{e}}_m$. Hence for all large m there exists a solution satisfying the required estimates.

4. Results

The following theorem is the main result of this paper. It establishes the existence of infinitely many multi–bubble solutions for the coupled critical system (2) under the symmetry and non-degeneracy assumptions (A1)–(A3).

Theorem 1.

Assume that V and K satisfy (A1)–(A3) with $N=3$. Then there exists an integer $m_0$ such that for every integer $m\ge m_0$ the system (2) admits a positive solution $(u_m,v_m)$ of the form:

$$u_m\left(y\right)={\displaystyle \sum _{j=1}^m}\xi \left(y\right)U_{x_j,{\lambda }_m}\left(y\right)+{\varphi }_{1,m}\left(y\right),v_m\left(y\right)={\displaystyle \sum _{j=1}^m}\xi \left(y\right)U_{x_j,{\lambda }_m}\left(y\right)+{\varphi }_{2,m}\left(y\right),$$

(68)

where $U_{x,\lambda }\left(y\right)=3^{1/4}{\lambda }^{1/2}{\left(1+{\lambda }^2{|y-x|}^2\right)}^{-1/2}$ is the standard bubble (see (3)) and ξ is a smooth cut–off function supported in a fixed neighbourhood of $(r_0,y_{30})$ (see Section 2). The centres ${\left\{x_j\right\}}_{j=1}^m$ are equally spaced on a circle of radius ${\overline{r}}_m$ in the $y^{\prime }$-plane, i.e.,

$$x_j=\left({\overline{r}}_{m}cos{\theta }_j,{\overline{r}}_{m}sin{\theta }_j,{\overline{y}}_{3,m}\right),{\theta }_j={\displaystyle \frac{2(j-1)\pi }{m}},$$

(69)

with parameters $({\overline{r}}_m,{\overline{y}}_{3,m})\in {\mathbb{R}}^2$ and ${\overline{r}}_m>0$. The dilation parameters ${\lambda }_m>0$ satisfy the asymptotic law

$${\lambda }_m=Cm+o\left(m\right),C=\sqrt{{\displaystyle \frac{A_2}{|A_1|}}}>0,$$

(70)

where $A_1<0$ and $A_2>0$ are constants depending on V, K, β and the geometry of the lattice (see the expansion of the reduced energy in Section 3). The remainder terms ${\varphi }_{i,m}\in H_s$ ($i=1,2$) satisfy the orthogonality conditions

$${\displaystyle \sum _{j=1}^m}{\int }_{{\mathbb{R}}^3}{Z}_{x_j,{\lambda }_m}^4{Z}_{j,l}{\varphi }_{i,m}dy=0,l=1,2,3,$$

(71)

and the weighted norm estimate

$$\parallel {\varphi }_{i,m}{\parallel }_{*}=O\left({\displaystyle \frac{1}{{\lambda }_m}}\right),$$

(72)

with ${\parallel ·\parallel }_{*}$ defined in Section 2. The centres converge to the non-degenerate critical point of K at the rate $O(1/m)$:

$$({\overline{r}}_m,{\overline{y}}_{3,m})⟶(r_0,y_{30})\mathrm{as}m\to \infty .$$

(73)

Finally, the energy functional I (defined in Section 3) satisfies

$$I(u_m,v_m)=m\left(A_0+o\left(1\right)\right),A_0={\displaystyle \frac{1}{2}}{\int }_{{\mathbb{R}}^3}{|\nabla U_{0,1}|}^2-{\displaystyle \frac{1}{6}}{\int }_{{\mathbb{R}}^3}{U}_{0,1}^6>0,$$

(74)

so the solutions are infinitely many and have arbitrarily large energy.

Theorem 1 provides a rigorous construction of multi-bubble solutions for the coupled critical system in ${\mathbb{R}}^3$. Several mathematical features deserve emphasis:

• Concentration pattern. The m bubbles are centred on a circle whose radius ${\overline{r}}_m$ approaches the critical radius $r_0$ of K. The symmetry of the configuration (discrete rotations and reflection) forces the remainders ${\varphi }_{i,m}$ to lie in the subspace $H_s$, which is essential to avoid the translational and rotational zero modes of the linearised operator.
• Precise scaling law. The dilation parameter ${\lambda }_m$ scales linearly with m: ${\lambda }_m=Cm+o\left(m\right)$. The constant C is determined by the competition between the self-energy term $A_1/{\lambda }^2$ (which comes from the expansion of V and K around the critical point) and the interaction term $A_2{m}^2/{\lambda }^2$ (which arises from the overlap of distinct bubbles). The condition $A_1<0$ (confinement) and $A_2>0$ (repulsive interaction) guarantees a unique positive solution for $\lambda$ for each large m.
• Energy asymptotics. The leading order of the energy is $m{A}_0$, where $A_0$ is the energy of a single bubble in the homogeneous limit ($V\equiv 0$, $K\equiv 1$). Corrections of order $1/{\lambda }^2\sim 1/m^2$ are controlled by the constants $A_1$ and $A_2$. In particular, the energy grows linearly with m, so the solutions are distinct and have arbitrarily large energy as $m\to \infty$.
• Role of the subcritical coupling. The coupling term with exponent $p\in (2,6)$ is subcritical; it does not affect the leading order asymptotics. It only contributes to the lower–order constant $A_1$ together with V and K. The proof shows that the coupling does not destroy the invertibility of the linearised operator (Lemma 1) nor the contraction argument (Proposition 1), because its contribution to the operator norm is $O(1/\lambda )$.

The proof of Theorem 1 is organised in several steps, which we now detail.

1.

Construction of the approximate solution. For parameters $(\overline{r},{\overline{y}}_3)$ close to $(r_0,y_{30})$ and $\lambda$ large, define the centres $x_j$ as in (4.2) and the truncated bubble $Z_{x,\lambda }$ as in (2.7). The approximate solution is $Z={\sum }_{j=1}^m{Z}_{x_j,\lambda }$.

2.

Linearised operator and invertibility. Define the weighted norms ${\parallel ·\parallel }_{*}$ and ${\parallel ·\parallel }_{**}$ as in (2.8)-(2.9). The linearised operator ${\mathbf{L}}_Z$ is decomposed as ${\mathbf{L}}_0+\mathbf{C}$, where ${\mathbf{L}}_0$ is diagonal and $\mathbf{C}$ is the subcritical coupling. Lemma 1 establishes that ${\mathbf{L}}_Z$ is invertible on the orthogonal subspace E (defined in (2.21)) with a uniform bound on its inverse.

3.

Solution of the auxiliary problem. For fixed parameters, consider the modified equation (3.1) with Lagrange multipliers $c_l$. Using the quadratic estimate for $\mathcal{N}$ (Lemma 2) and the error estimate for ${\mathbf{l}}_m$ (Lemma 3), Proposition 1 gives a unique solution $\Phi \in E$ with ${\parallel \Phi \parallel }_{*}\le C/\lambda$ and $|c_l|\le C/\lambda$.

4.

Reduced energy and Pohozaev identities. Define the reduced energy ${\mathcal{J}}_m(\overline{r},{\overline{y}}_3,\lambda )=I(Z+{\varphi }_1,Z+{\varphi }_2)$, with I as in (3.5). Its asymptotic expansion is derived in (3.11). The local Pohozaev identities (Proposition A1) imply that if $(\overline{r},{\overline{y}}_3,\lambda )$ is a critical point of ${\mathcal{J}}_m$, then $c_l=0$ for all l, so the modified equation reduces to the original system.

5.

Critical point conditions and existence. The critical point equations reduce to $\nabla K(\overline{r},{\overline{y}}_3)=O(1/m)$ and the scaling law $\lambda =Cm+o\left(m\right)$, where $C=\sqrt{A_2/\left|A_1\right|}$. The implicit function theorem guarantees a unique solution $({\overline{r}}_m,{\overline{y}}_{3,m},{\lambda }_m)$ for each sufficiently large m. Positivity follows from the smallness of $\Phi$ in the ${\parallel ·\parallel }_{*}$ norm, which ensures $Z+{\varphi }_i>0$.

This completes the proof of Theorem 1.

5. Conclusions

In this work, we have extended the Lyapunov–Schmidt reduction method to a critical elliptic system with competing potentials in ${\mathbb{R}}^3$. The main result, Theorem 1, provides the existence of infinitely many multi–bubble solutions for the system (2) under natural symmetry and non-degeneracy assumptions on the potentials V and K. The construction relies on several key ideas:

• A careful choice of weighted norms ${\parallel ·\parallel }_{*}$ and ${\parallel ·\parallel }_{**}$ that capture the multi–scale interaction of the bubbles.
• The decomposition of the linearised operator into a diagonal part ${\mathbf{L}}_0$ (coming from the scalar critical operator) and a coupling part $\mathbf{C}$ that is shown to be a contraction for large $\lambda$ because $p<6$.
• The use of local Pohozaev identities to eliminate the Lagrange multipliers $c_l$ that enforce orthogonality; this reduces the problem to finding critical points of a finite–dimensional energy ${\mathcal{J}}_m$.
• An asymptotic expansion of ${\mathcal{J}}_m$ that yields a transparent balance between self–energy (term $A_1/{\lambda }^2$) and interaction (term $A_2{m}^2/{\lambda }^2$), leading to the scaling $\lambda \sim Cm$.

The analysis shows that the subcritical coupling does not obstruct the concentration mechanism; it only contributes to the lower–order terms in the reduced energy. The solutions we obtain are genuine (not just approximate) and have energy growing linearly with the number of bubbles.

• Limitations.

The construction relies heavily on the symmetry of the configuration (discrete rotations and reflection) to work in the subspace $H_s$ and to avoid dealing with the full kernel of $L_Z$. The potentials V and K are assumed to be radially symmetric in the first two variables and even in $y_2$. Moreover, the exact constant C in the scaling ${\lambda }_m\sim Cm$ is not determined explicitly because it depends on the constants $A_1$ and $A_2$, which in turn depend on V, K, $\beta$ and the geometry of the lattice. A more detailed asymptotic expansion would be required to compute C explicitly.

• Future directions.

Several natural extensions of this work are possible:

• Relax the symmetry assumptions; one could consider almost–symmetric configurations or bubbles centred at points that are not exactly on a circle but are close to a critical manifold of K.
• Study the stability (orbital or asymptotic) of the constructed multi-bubble solutions. This would require a spectral analysis of the linearised operator around the solution and is a challenging problem even for scalar equations.
• Investigate the case where the coupling term is also critical (e.g., $p=6$). In that regime, the coupling would contribute at the same order as the leading terms and might drastically change the picture, possibly leading to new phenomena such as synchronisation, symmetry breaking, or even the formation of bound states with different scaling laws. This would require a refined analysis and is left for future work.
• Apply the same reduction technique to other dimensions ($N=4,5$) where the critical exponent changes, or to systems with more than two components.

We hope that the methods developed here will be useful for studying concentration phenomena in other elliptic systems of physical interest, such as coupled Gross–Pitaevskii equations for Bose–Einstein condensates or systems arising in nonlinear optics.

Appendix A.1. Derivation of the Local Pohozaev Identity

Let $(u,v)$ be a smooth solution of (2) and let $D\subset {\mathbb{R}}^3$ be a smooth bounded domain. For an index $i=1,2,3$, multiply the first equation by ${\partial }_{i}u$ and the second by ${\partial }_{i}v$, sum, and integrate over D. Using the divergence theorem and the fact that $-\Delta u{\partial }_{i}u=div({\partial }_{i}u\nabla u)-{\displaystyle \frac{1}{2}}{\partial }_i{|\nabla u|}^2$, we obtain after standard manipulations:

$$\begin{array}{cc}\hfill & {\int }_D{\partial }_{i}V(u^2+v^2)-{\displaystyle \frac{1}{6}}{\int }_D{\partial }_{i}K(u^6+v^6)-{\displaystyle \frac{\beta }{p}}{\int }_D{\partial }_{i}K{\left(uv\right)}^{{\displaystyle \frac{p-1}{2}}}\hfill \\ \hfill & ={\int }_{\partial D}\left[{\displaystyle \frac{1}{2}}{\left(\right|\nabla u|}^2+{|\nabla v|}^2){\nu }_i-({\partial }_{i}u\nabla u+{\partial }_{i}v\nabla v)·\nu \right]dS\hfill \\ \hfill & +{\displaystyle \frac{1}{6}}{\int }_{\partial D}K(u^6+v^6){\nu }_{i}dS+{\displaystyle \frac{\beta }{p}}{\int }_{\partial D}K{\left(uv\right)}^{{\displaystyle \frac{p-1}{2}}}{\nu }_{i}dS\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{\int }_{\partial D}V(u^2+v^2){\nu }_{i}dS,\hfill \end{array}$$

(A1)

where $\nu$ is the outward unit normal on $\partial D$ and ${\nu }_i$ its i-th component. The right-hand side is denoted ${\mathcal{B}}_i(u,v,\partial D)$. The derivation uses integration by parts and the fact that the nonlinear terms are exact derivatives: ${\partial }_i\left(K{u}^6\right)=\left({\partial }_{i}K\right)u^6+6K{u}^5{\partial }_{i}u$, etc., and the equations provide $-\Delta u+Vu=K(u^5+\beta u^{{\displaystyle \frac{p-1}{2}}}{v}^{{\displaystyle \frac{p-1}{2}}})$.

Appendix A.2. Application to the Approximate Solution

Consider the approximate solution $(Z+{\varphi }_1,Z+{\varphi }_2)$ constructed in Proposition 1. Fix a large ball $D=B_R\left(0\right)$ with R sufficiently large so that all centres $x_j$ lie well inside D and the cut–off $\xi$ is identically 1 on a neighbourhood of the centres. For large m (hence large $\lambda \sim m$), the bubbles $U_{x_j,\lambda }$ decay exponentially outside a region of size $O(1/\lambda )$ around each centre. Consequently, on $\partial D$ the functions Z, ${\varphi }_i$ and their derivatives are exponentially small in $\lambda$. Moreover, the remainder $\Phi$ satisfies ${\parallel \Phi \parallel }_{*}=O(1/\lambda )$, which implies pointwise bounds $|{\varphi }_i\left(y\right)|\le C{\lambda }^{-1}{\sum }_j{\lambda }^{1/2}(1+\lambda |y-x_j{\left|\right)}^{-1}$. On $\partial D$, each term in the sum is $O\left(e^{-\kappa \lambda R}\right)$, for some $\kappa >0$, hence the boundary integrals ${\mathcal{B}}_i(Z+{\varphi }_1,Z+{\varphi }_2,\partial D)$ tend to zero as $m\to \infty$ (exponentially fast).

Appendix A.3. Relation With The Lagrange Multipliers

Recall that the modified equation solved by $\Phi$ is

$${\mathbf{L}}_Z\Phi =\mathcal{N}(\Phi )+{\mathbf{l}}_m+{\displaystyle \sum _{l=1}^3}{c}_l{\displaystyle \sum _{j=1}^m}{Z}_{x_j,\lambda }^4{Z}_{j,l},$$

(A2)

where $Z_{j,l}$ are the derivatives of Z with respect to the parameters $\lambda ,\overline{r},{\overline{y}}_3$. Taking the $L^2$ inner product of this equation with a test function $\psi$ that is a linear combination of the kernel functions, and integrating by parts, one can show that the left-hand side becomes, up to boundary terms, the variation of the energy I at $(Z+{\varphi }_1,Z+{\varphi }_2)$ with respect to the corresponding parameter. More precisely, using the fact that ${\displaystyle \frac{\partial Z}{\partial t}}$ (with $t=\lambda ,\overline{r},{\overline{y}}_3$) belongs to the kernel of ${\mathbf{L}}_Z$, we obtain

$${\displaystyle \frac{\partial {\mathcal{J}}_m}{\partial t}}={\displaystyle \sum _{l=1}^3}{c}_l{M}_{l,t}+\mathrm{boundary}\mathrm{terms},$$

(A3)

where $M_{l,t}$ is the inner product of the kernel function ${\sum }_j{Z}_{x_j,\lambda }^4{Z}_{j,l}$ with ${\displaystyle \frac{\partial Z}{\partial t}}$. The matrix $M_{l,t}$ is non-degenerate (it is diagonal with entries of order ${\lambda }^2$). Since the boundary terms vanish as $m\to \infty$ (by the exponential decay argument above), a critical point of ${\mathcal{J}}_m$ (i.e., ${\displaystyle \frac{\partial {\mathcal{J}}_m}{\partial t}}=0$) forces $c_l=0$ for all $l=1,2,3$ in the limit $m\to \infty$; for finite large m the same conclusion holds by a perturbation argument because the boundary terms are $o\left(1\right)$ and the matrix is invertible.

Thus, we have the following rigorous statement:

Proposition A1.

Let $(Z+{\varphi }_1,Z+{\varphi }_2)$ be the family of approximate solutions constructed in Proposition 1, depending smoothly on the parameters $t\in \{\lambda ,\overline{r},{\overline{y}}_3\}$. Then the reduced energy ${\mathcal{J}}_m\left(t\right)=I(Z+{\varphi }_1,Z+{\varphi }_2)$ satisfies:

$${\displaystyle \frac{\partial {\mathcal{J}}_m}{\partial t}}={\displaystyle \sum _{l=1}^3}{c}_l{\langle \sum _j{Z}_{x_j,\lambda }^4{Z}_{j,l},{\displaystyle \frac{\partial Z}{\partial t}}\rangle }_{L^2}+O\left(e^{-\kappa \lambda }\right),$$

(A4)

where the error term comes from boundary integrals on a large ball and is exponentially small in λ (hence in m) with a constant $\kappa >0$. The matrix ${\left(\langle {\sum }_j{Z}_{x_j,\lambda }^4{Z}_{j,l},{\displaystyle \frac{\partial Z}{\partial t_k}}\rangle \right)}_{l,k}$ is invertible, with leading order ${\lambda }^{2}diag({\gamma }_1,{\gamma }_2,{\gamma }_3)$ for some ${\gamma }_l>0$. Consequently, if t is a critical point of ${\mathcal{J}}_m$ (i.e., ${\nabla }_t{\mathcal{J}}_m=0$), then for sufficiently large m the Lagrange multipliers satisfy $c_l=0$ for $l=1,2,3$, and therefore $(Z+{\varphi }_1,Z+{\varphi }_2)$ is an exact solution of the original system (2).

Proof. 

The proof follows by inserting the modified equation into the variation of the energy functional, integrating by parts over a large ball $B_R$, and using the exponential decay of the bubbles to estimate the boundary terms. The details are analogous to the scalar case treated in [1] and [2]; the coupling terms are subcritical and do not affect the leading order cancellations. The invertibility of the matrix is a consequence of Lemma 1 and the non–degeneracy of the limiting bubble. Hence the claim holds. □

This proposition justifies the reduction step described in Section 3: a critical point of the reduced energy yields a genuine solution of the original system, because the Lagrange multipliers vanish.