EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 6102 ISSN 1307-5543 – ejpam.com Published by New York Business Global Boundary Blowing Up Solutions for an Elliptic Neumann Problem with Nearly Critical Exponent Rakan Almushahhin1, Mohamed Ben Ayed1,∗ 1 Department of Mathematics, College of Science, Qassim University, Buraydah 51542, Saudi Arabia Abstract. In this paper, we investigate the nonlinear problem (Pε) : −∆u + V (x)u = fu n+2 n−2−ε, u > 0 in Ω and ∂u/∂ν = 0 on ∂Ω, where Ω is a bounded regular domain in Rn, with n ≥ 4, ε is a small positive parameter, V and f are smooth positive functions on Ω. Under certain conditions involving the function f and the mean curvature of the boundary, we construct boundary blowing up solutions of the problem (Pε) which converge weakly to 0 and blow up at some critical points of fb := f|∂Ω. This existence of solutions leads to a multiplicity result for (Pε). The proof of these results involves expanding the gradient of the associated functional and testing the equation with suitable vector fields. This process imposes constraints on the concentration parameters, and a careful analysis of these constraints leads to the conclusions presented. 2020 Mathematics Subject Classifications: 35A15, 35J20, 35J25 Key Words and Phrases: Partial Differential Equations, Neumann elliptic problems, Critical Sobolev exponent 1. Introduction In the last decades, there has been a great deal of interest in studying the following problem (Pλ,q) { −∆u+ λu = uq, u > 0 in Ω, ∂u ∂ν = 0, on ∂Ω, where Ω is a smooth and bounded open set of Rn with n ≥ 3, q > 1 and λ is a positive real number. Problem (Pλ,q) is a well-known example encountered in various applied sciences. For instance, it can be interpreted as the stationary problem arising in a Keller-Segel chemo- taxis model [1, 2] , originally developed to describe cell migration in response to chemical cues. Over time, this model has found broad application in engineering, supporting ad- vances in targeted drug delivery, tissue engineering, microfluidic systems, and the design of bio-inspired robots guided by chemotactic behavior [3–5]. ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.6102 Email addresses: 441112378@qu.edu.sa (R. Almushahhin), M.benayed@qu.edu.sa (M. Ben Ayed) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 2 of 31 From a mathematical standpoint, problem (Pλ,q) is an interesting model due to its solutions often exhibiting the bubbling phenomenon, with the center of the bubble located on the boundary. This boundary bubbling highlights a strong interaction between the geometry of the boundary and the solutions of (Pλ,q). A large body of research has explored problem (Pλ,q) when the exponent q is fixed and λ is treated as the parameter. A noteworthy aspect of the problem (Pλ,q) is the existence of solution families, denoted uλ,q, that exhibit blow up phenomena as λ changes. More precisely, these solutions blow up around certain points within Ω or on its boundary ∂Ω, while remaining negligibly small elsewhere. For the subcritical case, where 1 < q < (n + 2)/(n − 2), the only solution to (Pλ,q) for small λ is the constant one. However, as λ grows larger, non-constant solutions arise, which blow up at one or more points as λ→ ∞ [6]. The least energy solution experiences, for large λ, a blow-up at a boundary point where the mean curvature of the boundary is maximized [6–9]. Numerous works, such as [6, 10–13], have analyzed higher-energy solutions of (Pλ,q) that exhibit this asymptotic profile, whether blow up at boundary or interior points, as λ → ∞. In particular, solutions with any desired number of blow up points, both interior and boundary, have been shown to exist. The case when q = (n + 2)/(n − 2), the critical exponent, differs significantly. For n ∈ {4, 5, 6} and small λ, (Pλ,q) admits non-constant solutions [14–16]. On the other hand, the limiting equation of problem (Pλ,q), which arises when studying the asymptotic behavior of the least-energy solution as λ→ ∞, does not have any solutions. Nevertheless, least-energy solutions uλ,q still exist for large λ, and concentration appears in the form [17, 18] ωaλ,µλ (x), where aλ ∈ ∂Ω behaves as in the subcritical case, converging to the point that maximizes the mean curvature of the boundary. In this context, for any a ∈ Rn and µ ∈ (0,∞), the function ωa,µ represents the standard bubble defined by ωa,µ(x) := β0 µ(n−2)/2 (1 + µ2|x− a|2)(n−2)/2 , where β0 := [n(n− 2)](n−2)/4 (1) which are the only solutions [19] of −∆u = u n+2 n−2 , u > 0 in Rn. Higher-energy solutions of (Pλ,q) with concentration on the boundary as λ→ ∞ have been constructed in several studies, such as [17, 18, 20–27] and the references therein. Unlike the subcritical scenario, at least one blow up point must lie on the boundary [28]. Another interesting avenue of research for problem (Pλ,q) involves studying blow up phenomena by fixing λ while letting the exponent q approaches the critical exponent, i.e., q = n+2 n−2 ± ε, where ε is a small positive parameter. This was first explored by Rey and Wei. For n ≥ 4 and q = n+2 n−2 + ε, they demonstrated the existence of a solution that blows up at a boundary point where the mean curvature is maximized [29]. They also showed R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 3 of 31 the existence of a solution that blows up at a boundary point where the mean curvature is minimized when q = n+2 n−2 − ε and Ω is not convex [29]. In dimension 3, they found a solution with single interior blow-up point [30]. More recently, it was shown that for n ≥ 4 and q = n+2 n−2 + ε, there are no solutions that exhibit blow-up exclusively at interior points when ε is a small positive number [31]. Furthermore, in [32], the authors extended the problem by replacing the constant λ with a function V and studied the problem (PV,ε) { −∆u+ V u = u n+2 n−2 −ε, u > 0 in Ω, ∂u ∂ν = 0, on ∂Ω, where Ω is a smooth, bounded subset of Rn, n ≥ 6, V is a positive C2-function on Ω, and ε is a small positive parameter. They constructed interior bubbling solutions, where the interior blow-up points of these solutions converge, as ε → 0, to the critical points of the function V . More recently, in [33], the authors studied the case where a function f is introduced in front of the nonlinear term. More precisely, they considered the following problem (Pε) { −∆u+ V u = fup−ε, u > 0 in Ω, ∂u ∂ν = 0 on ∂Ω, where Ω is a smooth bounded domain of Rn, n ⩾ 4, V and f are positive C2-functions on Ω, ε is a small positive parameter and p+1 = (2n)/(n−2) is the critical Sobolev exponent for the embedding H1(Ω) ↪→ Lq(Ω). They constructed solutions of (Pε) with multi-blow up points located in the interior. A natural question arises: is it possible to construct solutions with boundary blow-up points? As mentioned earlier, this question was partially addressed by Rey and Wei [29], who constructed a solution with a single boundary blow-up point when Ω is non-convex and the function f is equal to 1. In this paper, our objective is to construct solutions to (Pε) with multiple boundary blow-up points and to present a multiplicity result for this problem. More precisely, the main results of our work are stated as follow: Theorem 1. Let n ≥ 4 and b1, · · · , bN be N non-degenerate critical points of f1 := f|∂Ω. We assume that c5 f(bk) ∂f ∂ν (bk)− (c1 2 − c4 ) H(bk) > 0 ∀ k ∈ {1, · · · , N}, (2) where H is the mean curvature of the boundary ∂Ω and c1, c4 and c5 are defined in Lemmas 4, 6 and 8 respectively. Then, there exists ε0 > 0 such that, for any ε ∈ (0, ε0) and for any subset {bi1 , · · · , biℓ} ⊂ {b1, · · · , bN}, Problem (Pε) has a solution uε which converges weakly to zero and blows up at the points bij ’s with the following properties lim ρ→0 lim ε→0 ∫ Ω∩B(bij ,ρ) fu2n/(n−2) ε = f(bij )Sn for each j ∈ {1, · · · , ℓ}, R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 4 of 31 where Sn is an universal constant defined in (22). More precisely, for any ℓ ≤ N , there exist µ1,ε,..., µℓ,ε having the same order as ε−1/2 for n ≥ 5, and as ε−1/2| ln ε|1/2 for n = 4 and ℓ points aj,ε → bij for all j such that ∣∣∣∣∣∣∣∣uε − ℓ∑ j=1 ωaj,ε,µj,ε ∣∣∣∣∣∣∣∣ H1(Ω) → 0, as ε→ 0 where the function ωa,µ is defined in (1). Theorem 1 allows us to obtain the following multiplicity result for problem (Pε) in relation to the number of non-degenerate critical points of the restriction of the function f on the boundary of Ω. Theorem 2. Let n ≥ 4 and assume that the restriction of f to the boundary has N non- degenerate critical points b1, · · · , bN satisfying assumption (2). Then, for small positive ε, the number of solutions to (Pε) that blow up on the boundary is at least 2N − 1. Remark 1. As examples of functions satisfying the assumption (2), let Ω := B(0, 1) and g be a positive C2-function on Ω such that g|∂Ω has only non-degenerate critical points. Let f(x) := g(x) + γ|x|2. It easy to see that ∂f ∂ν (y) = ∂g ∂ν (y) + 2γ ∀y ∈ ∂Ω and f|∂Ω = g|∂Ω + γ. Since the function H is constant on ∂Ω, it follows that, for γ large, the function f satisfies the assumption (2). The proof of our results relies on certain balancing conditions satisfied by the concen- tration parameters, which are relationships that ensure equilibrium between the various factors influencing the blow-up behavior of the solutions. These conditions are derived by performing an asymptotic expansion of the gradient of the Euler-Lagrange functional associated with the problem and testing the equation with appropriate vector fields. This process leads to constraints on both the concentration points and the corresponding blow- up rates of the solution. Through a careful analysis of these conditions, we derive our results. Note that traditional blow-up analysis methods depend on precise point-wise C0- estimates and the use of Pohozaev identities. In contrast, our method, used in this paper, deviates from these techniques. Bypassing the need for point-wise estimates and Pohozaev identities, our method holds significant promise for handling non-compact variational prob- lems that involve more intricate blow-up behaviors, as the existence of non-simple blow-up points. Moreover, the method developed in this paper is specifically tailored to the varia- tional problem and does not directly extend to non-variational settings. R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 5 of 31 The paper is structured as follows: In Section 2, we introduce the necessary prelim- inaries for studying the problem (Pε). Section 3 focuses on the analysis of the infinite- dimensional part of the solutions. In Section 4, we carry out an asymptotic expansion of the gradient of the Euler-Lagrange functional associated with (Pε). Section 5 presents the proof of our main results and Section 6 explores possible avenues for future research. Finally, the proofs rely on some technical facts, which are provided in the appendix in Section 7 for the reader’s convenience. 2. Preliminaries In this section, we proceed with the parametrization of the variational problem under consideration. Indeed, problem (Pε) is a variational one and its solutions are the positive critical points of the functional Jε(u) := 1 2 ∫ Ω |∇u|2 + 1 2 ∫ Ω V u2 − n− 2 2n− ε(n− 2) ∫ Ω f |u| 2n n−2 −ε, u ∈ H1(Ω). (3) The space H1(Ω) is equipped with the scalar product and its corresponding norm defined by: ⟨u1, u2⟩ := ∫ Ω ∇u1∇u2 + ∫ Ω V u1u2; ∥u∥2 := ∫ Ω |∇u|2 + ∫ Ω V u2. Since V is a bounded positive continuous function on Ω, it follows that this norm is equivalent to the standard norm of H1(Ω). Observe that, if uε is a solution of (Pε), satisfying uε ⇀ 0 (converges weakly to zero), by the concentration compactness principle [34], it follows that uε has to be close to some bubbles as ε → 0, that is, there exist q ∈ N, µ1, · · · , µq −→ ∞ (as ε → 0) and a1, . . . , aq ∈ Ω such that, as ε→ 0,∥∥∥∥∥uε − q∑ i=1 f(ai) (2−n)/4ωai,µi ∥∥∥∥∥→ 0, and µi µj + µj µi + µiµj |ai − aj |2 −→ ∞. In this paper, we want to construct some solutions blowing up at some boundary points. To this aim, we introduce the following set: Let n ⩾ 4, η0 be a small positive real, γ0 be a fixed small positive constant and q ∈ N, we define ϑ (q, γ0, η0) := { (a, µ, α) ∈ (∂Ω)q× ( η−1 0 ,∞ )q × (0,∞)q : |ai − aj | ⩾ γ0 ∀i ̸= j, ε lnµi < η0 and |1− αif(ai) (n−2)/4| ⩽ η0 ∀i } . Furthermore, for a ∈ ∂Ω and µ > η−1 0 , we define Fa,µ := { v ∈ H1(Ω) : ∫ Ω ∇v∇ωa,µ = ∫ Ω ∇v∇∂ωa,µ ∂µ = ∫ Ω ∇v∇∂ωa,µ ∂τj = 0; 1 ≤ j ≤ n− 1 } (4) R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 6 of 31 where the τ ′js, for j = 1, . . . , n − 1, build an orthonormal system of coordinates of the tangent space to ∂Ω at the point a. In addition, for (a, µ, α) ∈ ϑ (q, γ0, η0), we introduce Fa,µ := ⋂ 1≤i≤q Fai,µi . (5) Our aim is to construct solutions u having the form u = ∑q i=1 αiωai,µi + v, with (a, µ, α) ∈ ϑ (q, γ0, η0) and v ∈ Fa,µ. 3. Study of the infinite-dimensional part of the solutions In this section, we take (a, µ, α) ∈ ϑ (q, γ0, η0) and we are going to study the v-part of the solution u. In the sequel, we denote by ũ := q∑ i=1 αiωai,µi , for (a, µ, α) ∈ ϑ (q, γ0, η0) . (6) Furthermore, for (a, µ, α) ∈ ϑ (q, γ0, η0), we denote by Bi := B(ai, γ0/2). It follows that ωai,µi ≤ c µ (n−2)/2 i in Ω \Bi and ũ = αiωai,µi +O (∑ j ̸=i c µ (n−2)/2 j ) in Bi. (7) In the following, we will use the estimate given below, the proof of which is derived by applying Taylor’s expansion. For t1, t2 ∈ R and γ ≥ 2, we have |t1 + t2|γ = |t1|γ + γ |t1|γ−2 t1t2 + 1 2 γ(γ − 1) |t1|γ−2 t22 + { O ( |t1|γ−3 |t2|3 + |t2|γ ) if γ > 3, O (|t2|γ) if γ ⩽ 3. (8) Thus, for u = ũ+v with v ∈ Fa,µ, using Eq. (8), , the expansion of Jε, defined by (3), is as follows Jε(u) = Jε(ũ)− Lε(v) + 1 2 Qε(v) +Rε(v), with (9) Qε(v) := ∥v∥2 − (p− ε) ∫ Ω fũp−ε−1v2, (10) Lε(v) := ∫ Ω fũp−εv, (11) Rε(v) = o ( ∥v∥2 ) , R′ ε(v) = o(∥v∥) and R′′ ε (v) = o(1). Next, we will prove the coercivity of the form Qε. More precisely, we have: Lemma 1. Let (a, µ, α) ∈ ϑ (q, γ0, η0). Then the following fact holds Qε(v) = Q(v) + o (∥∥v2∥∥) where Q(v) := ∥v∥2 − n+ 2 n− 2 q∑ i=1 ∫ Ω ω 4 n−2 ai,µiv 2. Proof. By using the following formula derived from Taylor’s expansion ∣∣∣∑ ti ∣∣∣γ = ∑ |ti|γ +O ∑ j ̸=i |titj |γ/2  ∀ti ∈ R, ∀γ ∈ (0, 2], R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 7 of 31 we derive that∫ Ω fũp−ε−1v2 = q∑ i=1 αp−ε−1 i ∫ Ω fω 4 n−2−ε ai,µi v2 + ∑ j ̸=i O (∫ Ω ( ωai,µi ωaj ,µj ) 2 n−2 v2 ) . Expanding f around ai, we obtain∫ Ω fω 4 n−2−ε ai,µi v2 = f (ai) ∫ Ω ω 4 n−2−ε ai,µi v2 +O (∫ Ω |x− ai|ω 4 n−2−ε ai,µi v2 ) = f (ai) ∫ Ω ω 4 n−2−ε ai,µi v2 + o ( ∥v∥2 ) . Furthermore, since ε lnµi is small, we have ω−ε ai,µi = β−ε 0 µ −εn−2 2 i ( 1 + n− 2 2 ε ln ( 1 + µ2 i |x− ai|2 )) +O ( ε2 ln2 ( 1 + µ2 i |x− ai|2 )) (12) = 1 + o(1). (13) Thus, since ∣∣∣1− αif (ai) (n−2)/4 ∣∣∣ is small, we get α 4 n−2−ε i ∫ Ω fω 4 n−2−ε ai,µi v2 = α 4 n−2 i f (ai) ∫ Ω ω 4 n−2 ai,µiv 2 + o ( ∥v∥2 ) = ∫ Ω ω 4 n−2 ai,µiv 2 + o ( ∥v∥2 ) . This completes the proof of Lemma 1. At this point, we require the following important result regarding the uniform coercivity of the quadratic form Q. Proposition 1. Let (a, µ, α) ∈ ϑ (q, γ0, η0). There exists β3 > 0 (independent of α, µ and a) such that Q(v) ⩾ β3∥v∥2 ∀v ∈ Fa,µ. The proof of this proposition will be presented in Subsection 7.2. Combining Lemma 1 and Proposition 1, we deduce that the quadratic form Qε is coercive, that is Qε(v) ⩾ 1 2 β3∥v∥2 ∀v ∈ Fa,µ. (14) Now, we need to estimate the norm of the linear form Lε. More precisely, we have: Lemma 2. Let a ∈ ∂Ω, µ be a large real satisfying ε lnµ is small. Then, for ψ and v satisfying ψ ∈ { ωa,µ, µ ∂ωa,µ ∂µ , 1 µ ∂ωa,µ ∂a } , v ∈ H1(Ω) with ∫ Ω ∇v · ∇ψ = 0, (15) we have ∣∣∣∣∫ Ω fω 4 n−2−ε a,µ ψv ∣∣∣∣ ⩽ c∥v∥ ( ε+ 1 µ ) . R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 8 of 31 Proof. Observe that, since ε lnµ is small, then Equation (13) holds true. Therefore, using the fact that |ψ| ⩽ cωa,µ, we get∫ Ω fω 4 n−2−ε a,µ ψv = f(a) ∫ Ω ω 4 n−2−ε a,µ ψv +O (∫ Ω |x− a|ω n+2 n−2 a,µ |v| ) = β−ε 0 µ−εn−2 2 f(a) ∫ Ω ω 4 n−2ψv +O ( ε ∫ Ω ω n+2 n−2 a,µ ln ( 1 + µ2|x− a|2 ) |v|+ ∥v∥ µ ) . (16) Easy computation leads to∫ Ω ω 2n n−2 a,µ ( ln ( 1 + µ2|x− a|2 ))γ ≤ c ∀ γ > 0 (17) which implies that ε ∫ Ω ω n+2 n−2 a,µ ln ( 1 + µ2|x− a|2 ) |v| ≤ cε∥v∥. (18) To estimate the first integral in the right hand side of (16), we distinguish three cases: • Case 1. If ψ = ωa,µ, using (15), Holder’s inequality, Lemma 7 and the continuity of the embedding H1(Ω) ↪→ L 2n−2 n−2 (∂Ω), we get∫ Ω ω 4 n−2 a,µ ψv = ∫ Ω ω n+2 n−2 a,µ v = ∫ Ω −∆ωa,µv = ∫ Ω ∇ωa,µ∇v − ∫ ∂Ω ∂ωa,µ ∂ν v = O [∫ ∂Ω |v| 2n−2 n−2 ] n−2 2n−2 [∫ ∂Ω ∣∣∣∣∂ωa,µ ∂ν ∣∣∣∣ 2n−2 n ] n 2n−2  = O ( ∥v∥ µ ) . (19) Combining Eqs. (16), (18) and (19), the proof of the lemma is completed in the case where ψ = ωa,µ. • Case 2. If ψ = µ (∂ωa,µ/∂µ), it holds:∫ Ω ω 4 n−2 a,µ ψv = n− 2 n+ 2 ∫ Ω −∆ ( µ ∂ωa,µ ∂µ ) v = n− 2 n+ 2 (∫ Ω ∇ ( µ ∂ωa,µ ∂µ ) ∇v − ∫ ∂Ω ∂ ∂ν ( µ ∂ωa,µ ∂µ ) v ) = O [∫ ∂Ω |v| 2n−2 n−2 ] n−2 2n−2 [∫ ∂Ω ∣∣∣∣µ∂2ωa,µ ∂ν∂µ ∣∣∣∣ 2n−2 n ] n 2n−2  = O ( ∥v∥ µ ) , by using Claims (i) and (ii) of Lemma 7. This completes the proof of of the lemma in the case where ψ = µ (∂ωa,µ/∂µ). • Case 3. If ψ = µ−1 (∂ωa,µ/∂aj), the proof can be done in the same way. Hence, we omit it. The proof of Lemma 2 is thereby completed. Now, we are ready to present the estimate of the linear form Lε defined in (11). Proposition 2. Let (a, µ, α) ∈ ϑ (q, γ0, η0). Then, we have |Lε(v)| := ∣∣∣∣∫ Ω fũp−εv ∣∣∣∣ ⩽ c∥v∥ ( ε+ ∑ 1 µi ) ∀v ∈ Fa,µ. R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 9 of 31 Proof. We will use the following formula, the proof of which follows from Taylor expansion. For tj > 0, (∑ tj )γ = ∑ tγj + ∑ i ̸=j O ( (titj) γ/2 ) if γ ⩽ 2, O ( tγ−1 i tj ) if γ > 2. (20) Observe that, if n ⩾ 6, it follows that p := (n+2) (n−2) ⩽ 2, and therefore, using (20) and (13), we deduce that Lε(v) = q∑ i=1 αp−ε i ∫ Ω fωp−ε ai,µi v + ∑ i ̸=j O (∫ Ω ( ωai,µiωaj ,µj )p/2 |v|) . Thus, using Lemma 2 and Holder’s inequality, we get |Lε(v)| ⩽ c∥v∥ ( ε+ ∑ 1 µi ) + c∥v∥ (∫ Ω ω n n−2 ai,µiω n n−2 aj ,µj )n+2 2n . But, since |ai − aj | ⩾ γ0, let Bk = B (ak, γ0/2), using (7), it follows that∫ Ω ω n n−2 ai,µiω n n−2 aj ,µj ⩽ c µ n 2 j ∫ Bi ω n n−2 ai,µi + c µ n 2 i ∫ Bj ω n n−2 aj ,µj + c (µiµj) n 2 ∫ Ω\(Bi∪Bj) dx ⩽ c ln (µiµj) (µiµj) n 2 . This completes the proof for n ⩾ 6. However, for n ⩽ 5, we need to estimate∫ Ω ω 4 n−2 ai,µiωaj ,µj |v| ⩽ c∥v∥ (∫ Ω ω 8n (n−2)(n+2) ai,µi ω 2n n+2 aj ,µj )n+2 2n . But, since |ai − aj | ⩾ γ0, using (7), we get∫ Ω ω 8n (n−2)(n+2) ai,µi ω 2n n+2 aj ,µj ⩽ c µ n(n−2) n+2 j ∫ Bi ω 8n (n−2)(n+2) ai,µi + c µ 4n n+2 i ∫ Bj ω 2n n+2 aj ,µj + c µ 4n n+2 i µ n(n−2) (n+2) j ⩽ c (µiµj) n(n−2) (n+2) . This completes the proof of the proposition. Proposition 3. Let (a, µ, α) ∈ ϑ (q, γ0, η0). Then, for ε small, there exists a unique v ∈ Fa,µ verifying ⟨∇Jε(ũ+ v), v⟩ = 0 ∀v ∈ Fa,µ. and ∥v∥ ⩽ c ( ε+ ∑ 1 µi ) . Proof. The proof follows from (9), (14) and Proposition 2 by using the implicit function theorem. 4. Asymptotic expansion of the gradient in the potential sets This section is devoted to the asymptotic expansion of the gradient of the functional Jε defined in (3). To this aim, by easy computation, we see that ⟨∇Jε(u), g⟩ = ⟨u, g⟩ − ∫ Ω f |u|p−ε−1ug ∀u, g ∈ H1(Ω). (21) We start by the expansion with respect to the variable α. R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 10 of 31 Proposition 4. Let (a, µ, α) ∈ ϑ (q, γ0, η0), ũ be defined in (6) and v ∈ Fa,µ where Fa,µ is defined in (5). Then, for ε small and i ∈ {1, . . . , q}, it holds ⟨∇Jε(ũ+ v), ωai,µi ⟩ = αiSn ( 1− µ −εn−2 2 i αp−ε−1 i f (ai) ) +O ( ∥v∥2 + 1 µi + ε ) , with Sn := 1 2 [n(n− 2)]n/2 ∫ Rn 1 (1 + |x|2)n dx. (22) Proof. Since v ∈ Fa,µ, using Lemmas 4 and 5, it follows that ⟨ũ+ v, ωai,µi ⟩ = q∑ j=1 αj 〈 ωaj ,µj ωai,µi 〉 = αi ( Sn +O ( 1 µi )) + q∑ j ̸=i O ( 1 (µiµj) (n−2)/2 ) . (23) Now, observe that |s+ t|γ(s+ t)z = |s|γsz + (γ + 1)|s|γtz +O(|s|γt2 + |t|γ+2) ∀ s, t ∈ R, |z| ≤ |s| and γ > 0. Thus, for each ψi satisfying |ψi| ⩽ cωai,µi , in Bi := B ( ai, γ0/2 ) ∩ Ω, using (7) and (13), it holds |ũ+ v|p−ε−1(ũ+ v)ψi = αp−ε i ωp−ε ai,µi ψi + (p− ε) (αiωai,µi) p−ε−1 ∑ j ̸=i αjωaj ,µj + v ψi +O ( ωp−1 ai,µi [∑ 1 µn−2 j + |v|2 ] + ∑ 1 µn j + |v|p+1−ε ) in Bi. (24) But, in Ω\Bi, we have |ũ+ v|p−ε |ψi| ⩽ c ( |v|p−ε + ∑ ωp aj ,µj ) |ψi| in Ω\Bi. (25) Thus, the integral in Eq. (21) becomes∫ Ω f |ũ+ v|p−ε−1(ũ+ v)ωai,µi = αp−ε i ∫ Bi fωp+1−ε ai,µi + (p− ε)αp−ε−1 i [∑ j ̸=i αj ∫ Bi fωp−ε ai,µi ωaj ,µj + ∫ Bi fωp−ε ai,µi v ] +O ( ∥v∥2 + ∑ 1 µn j + ∑ 1 µn−2 j ∫ Bi ωp−1 ai,µi + ∥v∥p−ε µ n−2 2 i + 1 µ n−2 2 i ∑∫ Ω\Bi ωp aj ,µj ) . (26) Using (13) and Lemma 5, we get∫ Bi ωp−ε ai,µi ωaj ,µj ⩽ c ∫ Bi ωp ai,µi ωaj ,µj ⩽ c (µiµj) (n−2)/2 . (27) In addition, using (13), (17) and Lemma 6, we get∫ Bi fωp+1−ε ai,µi = β−ε 0 µ −ε(n−2)/2 i f(ai) ∫ Ω ωp+1 ai,µi +O (∫ Ω\Bi ωp+1 ai,µi + ε ∫ Ω ωp+1 ai,µi ln ( 1 + µ2 i |x− ai|2 ) + ∫ Bi |x− ai|ωp+1 ai,µi ) R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 11 of 31 = β−ε 0 µ −ε(n−2)/2 i f(ai)Sn +O ( 1 µi + ε ) . (28) Furthermore, taking ψ = ωa,µ in Lemma 2 and using (7), we deduce that∫ Bi fωp−ε ai,µi v = ∫ Ω fωp−ε ai,µi v +O (∫ Ω\Bi ωp ai,µi |v| ) = O ( ∥v∥ [ ε+ 1 µi ]) +O ( ∥v∥ µ (n+2)/2 i ) = O ( ∥v∥ [ ε+ 1 µi ]) . (29) Combining (29), (28) and (27), the equation (26) becomes∫ Bi f |ũ+ v|p−ε−1(ũ+ v)ωai,µi = µ −ε(n−2)/2 i αp−ε i f (ai)Sn +O ( 1 µi + ε+ ∥v∥2 ) . (30) Combining (23) and (30), the proof of Proposition 4 follows. Next, we deal with the expansion with respect to µ. Proposition 5. Let (a, µ, α) ∈ ϑ (q, γ0, η0) and v ∈ Fa,µ. For ε small and i ⩽ q, we have〈 ∇Jε(ũ+ v),µi ∂ωai,µi ∂µi 〉 = n− 2 4 c6µ −εn−2 2 i αp−ε i f (ai) ε+ αi H (ai) µi (c1 2 − µ−ε i n−2 2 αp−ε−1 i f (ai) c4 ) − c5 µi ∂f ∂ν (ai)µ −ε i n−2 2 αp−ε i + O (n=4) ( lnµi µ2 i ) +O ( ∥v∥2 + ε2 + 1 µ2 i + ∑ 1 µn−2 j ) , where O (n=4) appears only if n = 4 and the constants c1, c4, c5 and c6 are defined in (55), (56), (57), (58) respectively. Proof. Observe that, using Lemmas 4, 5, 9 and the fact that v ∈ Fa,µ, we obtain 〈 ũ+ v, µi ∂ωai,µi ∂µi 〉 = q∑ j=1 αj 〈 ωaj ,µj , µi ∂ωai,µi ∂µi 〉 = c1 2 αi H (ai) µi +O ( 1 µ2 i + ∑ 1 µn−2 k ) + O (n=4) ( lnµi µ2 i ) . (31) In addition, using (24) and (25), we get∫ Ω f |ũ+ v|p−ε−1(ũ+ v)µi ∂ωai,µi ∂µi = αp−ε i ∫ Bi fωp−ε ai,µi µi ∂ωai,µi ∂µi + (p− ε)αp−ε−1 i ∫ Bi fωp−ε−1 ai,µi vµi ∂ωai,µi ∂µi +O (∑ j ̸=i 1 µ (n−2)/2 j ∫ Bi ωp ai,µi + ∑ j ̸=i 1 µn−2 j ∫ Bi ωp−1 ai,µi + ∑ 1 µn j + ∥v∥2 + 1 µ (n−2)/2 i ( ∥v∥p−ε + ∑∫ Ω ωp aj ,µj )) . (32) Notice that, the remainder term can be estimated as O ( q∑ j=1 1 µn−2 j + ∥v∥2 ) . R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 12 of 31 Furthermore, taking ψ = µ∂ωa,µ/∂µ in Lemma 2 and using (7), we get∫ Bi fωp−ε−1 ai,µi vµi ∂ωai,µi ∂µi = ∫ Ω fωp−ε−1 ai,µi vµi ∂ωai,µi ∂µi +O (∫ Ω\Bi ωp ai,µi |v| ) = O ( ∥v∥ [ ε+ 1 µi ]) +O ( ∥v∥ µ (n+2)/2 i ) = O ( ∥v∥ [ ε+ 1 µi ]) . (33) To complete the estimate of (32), using (12), we write∫ Bi fωp−ε ai,µi µi ∂ωai,µi ∂µi = β−ε 0 µ −εn−2 2 i [∫ Bi fωp ai,µi µi ∂ωai,µi ∂µi + n− 2 2 ε× × ∫ Bi fωp ai,µi µi ∂ωai,µi ∂µi ln ( 1 + µ2 i |x− ai|2 )] +O ( ε2 ∫ Bi ωp+1 ai,µi ln2 ( 1 + µ2 i |x− ai|2 )) . (34) Expanding f around ai and using Lemmas 6 and 8, we obtain∫ Bi fωp ai,µi µi ∂ωai,µi ∂µi = f(ai) ∫ Bi ωp ai,µi µi ∂ωai,µi ∂µi +∇f (ai) ∫ Bi (x− ai)ω p ai,µi µi ∂ωai,µi ∂µi +O (∫ Bi |x− ai|2 ωp ai,µi ) = f (ai) ( c4 H (ai) µi +O ( 1 µ2 i )) + ∂f ∂ν (ai) ( c5 µi +O ( 1 µ2 i )) +O ( 1 µ2 i ) . For the second integral in the right hand side of (34), expanding f around ai and using Lemma 8, we get ∫ Bi fωp ai,µi µi ∂ωai,µi ∂µi ln ( 1 + µ2 i |x− ai|2 ) = f (ai) ( −c6 2 +O ( 1 µi )) +O ( 1 µi ) . Thus, using (17), (34) becomes∫ Bi fωp−ε ai,µi µi ∂ωai,µi ∂µi =β−ε 0 µ −εn−2 2 i [ c4f(ai) H (ai) µi + c5 µi ∂f ∂ν (ai)− n− 2 4 c6f (ai) ε ] +O ( ε2 + 1 µ2 i ) . (35) This completes the proof of (32) and we get (by combining (35) and (33))∫ Ω f |ũ+ v|p−ε−1(ũ+ v)µi ∂ωai,µi ∂µi = β−ε 0 µ −εn−2 2 i [ c4f (ai) H (ai) µi + c5 µi ∂f ∂ν (ai) −n− 2 4 c6f (ai) ε ] +O ( ∥v∥2 + ε2 + 1 µi 2 + ∑ 1 µn−2 j ) . (36) Thus, Combining (36) and (31), the proof of Proposition 5 follows. We end this section by expanding the gradient of Jε with respect to the concentration point a. R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 13 of 31 Proposition 6. Let (a, µ, α) ∈ ϑ (q, γ0, η0) and v ∈ Fa,µ. Let i ∈ {1, . . . , q}, we denote by ( σk i ) , where 1 ⩽ k ⩽ n− 1, an orthonormal system of coordinates on the tangent space to the boundary ∂Ω at ai. It holds〈 ∇Jε(ũ+ v), 1 µi ∂ωai,µi ∂σk i 〉 = −c7µ −εn−2 2 i αp−ε i 1 µi ∂f ∂σk i (ai) +O ( ε2 + 1 µ2 i + ∑ 1 µn−2 j + ∥v∥2 ) where c7 is defined in (59). Proof. To simply the presentation, without loss of generality, we will assume that ai = 0 and the normal exterior vector νai = −en. By this choose, we deduce that the tangent space to the boundary ∂Ω at the point ai = 0 is Rn−1 × {0}. Let k ∈ {1, · · · , n− 1}, using Lemmas 4, 5, the proof is similar to the proof of Proposition 5. Here, we will give a sketch and precise the argument of the estimate of some integrals. Using Lemma 9 and the fact that v ∈ Fa,µ, we deduce that〈 ũ+ v, 1 µi ∂ωai,µi ∂ai,k 〉 = q∑ j=1 αj 〈 ωaj ,µj , 1 µi ∂ωai,µi ∂ai,k 〉 = O ( 1 µ2 i + ∑ 1 µn−1 j ) . For the other part of the gradient, using (24) and (25), we deduce that (32), (33) and (34) hold true by taking 1 µi ∂ωai,µi ∂ai,k instead of µi ∂ωai,µi ∂µi . Now, using Lemmas 6 and 8, it holds∫ Bi fωp ai,µi 1 µi ∂ωai,µi ∂ai,k = f (ai) ∫ Bi ωp ai,µi 1 µi ∂ωai,µi ∂ai,k +∇f (ai) ∫ Bi (x− ai)ω p ai,µi 1 µi ∂ωai,µi ∂ai,k +O (∫ Bi |x− ai|2ωp+1 ai,µi ) = c7 µi ∂f ∂xi (ai) +O ( 1 µ2 i ) ,∫ Bi fωp ai,µi 1 µi ∂ωai,µi ∂ai,k ln ( 1 + µ2 i |x− ai|2 ) = O ( 1 µi ) , by expanding f around ai. Hence, we obtain∫ Bi fωp−ε ai,µi 1 µi ∂ωai,µi ∂ai,k = β−ε 0 µ −ε(n−2)/2 i c7 µi ∂f ∂xk (ai) +O ( 1 µ2 i ) +O ( ε2 ) . This completes the proof of Proposition 6. 5. Proof of Theorems 1 and 2 Since Theorem 2 is a direct consequence of Theorem 1, it is sufficient to prove the latter. Adopting the proof strategy from [35], let N ∈ N, b1, . . . , bN be as defined in Theorem 1 and ε > 0 be small. We consider the set Dε,N := { (a, µ, α, v) ∈ (∂Ω)N × (0,∞)N × (0,∞)N ×H1(Ω) : |ai − bi| < √ ε ; M−1 1 ⩽ µiε ⩽M1; ∣∣∣1− αif (ai) (n−2)/4 ∣∣∣ < ε ln2 ε, v ∈ Fa,µ and ∥v∥ < √ ε } where M1 is a fixed large constant. Let J̃ε be the function defined by J̃ε : Dε,N −→ R ; Λ := (a, µ, α, v) 7−→ J̃ε(Λ) := Jε ( N∑ i=1 αiwai,µi + v ) . R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 14 of 31 There exists a biunivoque relation between the critical points of J̃ε and the ones of Jε. Proposition 7. Let Λ := (a, µ, α, v) ∈ Dε,N . u = ∑N i=1 αiwai,µi + v is a critical point of Jε if and only if Λ is a critical point of J̃ε, i.,e., there exists (γ, η, σ) ∈ ( Rn−1 )N × RN × RN such that the following system is satisfied: (Ak) ∂J̃ε ∂αk (Λ) = 0 ∀ k, (V ) ∂J̃ε ∂v (Λ) = ∑N k=1 ( ηk ∂g1,k ∂v (Λ) + σk ∂g2,k ∂v (Λ) + ∑n−1 j=1 γk,j ∂g3,k,j ∂v (Λ) ) , (Mk) ∂J̃ε ∂µk (Λ) = σk ∫ Ω ∇v · ∇ ( µk ∂2wak,µk ∂µ2 k ) + ∑n−1 j=1 γk,j ∫ Ω ∇v · ∇ ( 1 µk ∂2wak,µk ∂µk∂τk,j ) ,∀ k, (Tk) ∂J̃ε ∂τk,l (Λ) = σk ∫ Ω ∇v · ∇ ( µk ∂2wak,µk ∂µk∂τk,l ) + ∑n−1 j=1 γk,j ∫ Ω ∇v · ∇ ( 1 µk ∂2wak,µk ∂τk,j∂τk,l ) , ∀ k, ∀ l. (37) Proof. Observe that Dε,N is not an open set in (∂Ω)N × (0,∞)N × (0,∞)N ×H1(Ω) since the elements Λ of Dε,N have to satisfy the following orthogonality constraints: g1,k(Λ) := ∫ Ω ∇v · ∇wak,µk = 0, k ∈ {1, · · · , N}, g2,k(Λ) := ∫ Ω ∇v · ∇ ( µk ∂wak,µk ∂µk ) = 0, k ∈ {1, . . . , N}, g3,k,j(Λ) := ∫ Ω ∇v · ∇ ( 1 µk ∂wak,µk ∂τk,j ) = 0, k ∈ {1, . . . , N}, j ∈ {1, . . . , n− 1}. Therefore, it follows from the multiplier Lagrange theorem that Λ is a critical point of J̃ε in Dε,N if and only if there exist some constants γ ∈ ( Rn−1 )N , η ∈ RN and σ ∈ RN such that ∇J̃ε(Λ) = N∑ k=1 ( ηk∇g1,k(Λ) + σk∇g2,k(Λ) + ∑ 1≤j≤n−1 γk,j∇g3,k,j(Λ) ) . (38) Notice that ∇J̃ε(Λ) = ((∂J̃ε(Λ) ∂τk,1 ) k≤N , · · · , ( ∂J̃ε(Λ) ∂τk,n−1 ) k≤N , (∂J̃ε(Λ) ∂µk ) k≤N , (∂J̃ε(Λ) ∂αk ) k≤N , ∂J̃ε(Λ) ∂v ) . (39) Combining (38), (39) and the fact that the functions g1,k, g2,k and g3,k are independent of the variable α, we easily derive the result. To prove Theorem 1, we observe that Proposition 7 implies that it is sufficient to study the system (37) and demonstrate that (37) has a solution. First, for Λ ∈ Dε,N , let u = ∑N i=1 αiwaiµi +v, the definition of J̃ε implies that, for each k ∈ {1, . . . , N}, ∂J̃ε ∂v (Λ) = ∇Jε(u), ∂J̃ε ∂µk (Λ) = 〈 ∇Jε(u), αk ∂wak,µk ∂µk 〉 , ∂J̃ε ∂αk (Λ) = ⟨∇Jε(u), wak,µk ⟩ , ∂J̃ ∂τk,j (Λ) = 〈 ∇Jε(u), αk ∂wak,µk ∂τk,j 〉 , 1 ≤ j ≤ n− 1. (40) Second, using Proposition 3, for each (a, µ, α, 0) ∈ Dε,N , there exists v := v̄ε,a,µ,α ∈ Fa,µ such that〈 ∇Jε ( N∑ i=1 αiwaiµi + v ) , h 〉 = 0 ∀h ∈ Fa,µ and ∥v∥ ⩽ c ( ε+ ∑ 1 µi ) . (41) R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 15 of 31 Therefore, Eqs. (40) and (41) imply the existence of η ∈ RN , σ ∈ RN and γ ∈ ( Rn−1 )N such that ∇Jε(ū) = ∇J̃ε(Λ̄) = N∑ k=1 ηkwak,µk + σkµk ∂wak,µk ∂µk + n−1∑ j=1 γk,j 1 µk ∂wak,µk ∂τk,j  (42) where Λ̄ = (a, µ, α, v̄) and ū = ∑N i=1 αiwai,µi + v̄. Thus, (42) implies that the second equation of (37) is satisfied for each (a, µ, α, v̄) ∈ Dε,N . Hence, it remains to solve the three other equations. To do so, we start by giving the estimate of the multiplier Lagrange coefficients (η, σ, γ). Lemma 3. The multiplier Lagrange coefficients (η, σ, γ) found in (42) satisfy: |ηk| ⩽ cε ln2 ε; |σk| ⩽ cε; |γk,j | ⩽ cε3/2, for each k ∈ {1, . . . , N} and each j ∈ {1, . . . , n− 1}. Proof. Using Λ̄ := (a, µ, α, v̄) ∈ Dε,N and Propositions 4, 5 and 6, it follows that ⟨∇Jε(ū), wak,µk ⟩ = O ( ε+ ∣∣∣1− α 4/(n−2) i f (ai) ∣∣∣+ ε |ln ε| ) = O ( ε| ln ε|2 ) ,〈 ∇Jε(ū), µk ∂wak,µk ∂µk 〉 = O(ε),〈 ∇Jε(ū), 1 µk ∂wak,µk ∂τk,j 〉 = O ( ε |ak − bk|+ ε2 ) = O ( ε3/2 ) . Furthermore, for ψk, ψl ∈ N⋃ i=1 { wai,µi , µi ∂wai,µi ∂µi , 1 µi ∂wai,µi ∂τi,j , j ∈ {1, . . . , n− 1} } , using Lemmas 4, 5 and 9, we deduce that ⟨ψk, ψl⟩ = { c+O(ε) if k = l, O(ε) if k ̸= l, for some positive constant c. Thus, the scalar products of (42) with wai,µi , µi ∂wai,µi ∂µi and 1 µi ∂wai,µi ∂τi,j , respectively, give the following quasi-diagonal system: cηi + ∑ k O (ε (|γk|+ |σk|+ |ηk|)) = O ( ε ln2 ε ) , cσi + ∑ k O (ε (|γk|+ |σk|+ |ηk|)) = O(ε), cγij +O (ε (|γk|+ |σk|+ |ηk|)) = O ( ε3/2 ) , which implies the result. Now, we are ready to solve the equations (Ak,Mk, Tk) for k ∈ {1, . . . , N} defined in Proposition 7. R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 16 of 31 First, using Lemma 3 and Eq. (41), for Λ̄ = (a, µ, α, v̄), the equations (Ak,Mk, Tk) in the system (37) are equivalent to: ∂J̃ε ∂αk (Λ̄) = 0, ∀k ∈ {1, . . . , N}, µk ∂J̃ε ∂µk (Λ̄) = O (( |σk|+ ∑ |γk,j | ) ∥v̄∥ ) = O ( ε2 ) ,∀k ∈ {1, . . . ;N}, (43) 1 µk ∂J̃ε ∂τk,j (Λ̄) = O (( |σk|+ ∑ |γk,j | ) ∥v̄∥ ) = O ( ε2 ) , ∀k ∈ {1, . . . , N},∀j ∈ {1, . . . , n− 1}. Second, using (41) and Propositions 4, 5 and 6, we derive that the system (43) is equivalent to 1− µ −ε(n−2)/2 k αp−1−ε k f (ak) = O(ε), ∀k, n− 2 4 c6ε+ (c1 2 − c4 ) H (ak) µk − c5 µk 1 f (ak) ∂f ∂ν (ak) = { O ( ε2 ) if n ⩾ 5, O ( ε2| ln ε| ) if n = 4, ∀k, 1 µk ∇f1 (ak) = O ( ε2 ) , ∀k. (44) At this step, to solve the system (44), it is better to take a change of variables to obtain an easier system to solve. Notice that αk ∈ (0,∞) and µk ∈ (0,∞), however ak ∈ ∂Ω and therefore, we need to be move careful in the change of variables for ak. To be more precise, let y ∈ ∂Ω and( e′1, . . . , e ′ n−1,−νy ) be an orthonormal basis of Rn. In this basis, the tangent space to ∂Ω at y is Rn−1 × {0}. Written y = (y′, yn) ∈ Rn−1 × R, since Ω is a C2-domain, there exist ρ > 0 (small) and a C2-function g : Bn−1(0, ρ) ⊂ Rn−1 → R such that: • g(0) = 0, ∇g(0) = 0 and therefore |g (z′)| ⩽ c |z′|2 ∀z′, • Ω ∩Bn (y, ρ) = { (y′ + z′, yn + zn) ∈ Rn−1 × R : |(z′, zn)| < ρ and zn > g (z′) } , • ∂Ω ∩Bn (y, ρ) = { (y′ + z′, yn + zn) ∈ Rn−1 × R : |(z′, zn)| < ρ and zn = g (z′) } . Furthermore, assume that y is a critical point of f1 := f|∂Ω (the restriction of f on the boundary), for a ∈ ∂Ω ∩B(y, ρ), written a as a := (a′, an) = (y′ + z′, yn + g (z′)) with z′ ∈ Bn−1(0, ρ), (45) then it holds that ∇T f(a) = ∇f1(a) = D2f1(y) ((z ′, 0) , ·) +O ( |z′|2 ) . (46) Recall that (a, µ, α, 0) ∈ Dε,N which implies that ak is close to bk (which is a critical point of f1) and αkf (ak) n−2/4 is close to 1 for each k. Hence, let us consider the following change of variables: ρk :=1− αkf (bk) n−2/4 , k ∈ {1, . . . , N}, (47) 1 µk := [ c5 f (bk) ∂f ∂ν (bk)− (c1 2 − c4 ) H (bk) ]−1 n− 2 4 c6ε (1 + λk) , (48) ak := (b′k + z′k, (bk)n + g (z′k)) , k ∈ {1, . . . , N}, (49) by using the notation of (45). Using this change of variables, we get: f (ak) = f1 (ak) = f1 (bk) +O ( |ak − bk|2 ) = f1 (bk) +O ( |z′k| 2 ) , R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 17 of 31 H (ak) = H (bk) +O (|z′k|) , 1 f(ak) ∂f ∂ν (ak) = 1 f (bk) ∂f ∂ν (bk) +O (|z′k|) , 1− µ −ε(n−2)/2 k α 4 n−2−ε k f (ak) = 1− α 4 n−2 k f1 (bk) +O ( |z′|2 + ε |ln ε| ) = 4 n− 2 ρk +O ( ρ2k + |z′|2 + ε |ln ε| ) , (50) n− 2 4 c6ε+ 1 µk [(c1 2 − c4 ) H (ak)− c5 f (ak) ∂f ∂ν (ak) ] = n− 2 4 c6ε− 1 µk [ c5 f (bk) ∂f ∂ν (bk)− (c1 2 − c4 ) H (bk) ] +O ( 1 µk |z′| ) = −n− 2 4 c6ελk +O (ε |z′|) . (51) Thus, using (46), (50) and (51), the system (44) becomes equivalent to: ρk = O ( ρ2k + |z′|2 + ε |ln ε| ) , k ∈ {1, . . . , N}, λk = O ( |z′|+ (if n ⩾ 5) ε+ (if n = 4) ε |ln ε| ) D2f1 (bk) ((z ′, 0) , ·) = O ( ε+ |z′|2 ) , k ∈ {1, . . . , N}. (52) Since D2f1 (bk) is assumed to be non-degenerate, the last equation in (52) implies that |z′| ⩽ c ( ε+ |z′|2 ) , and therefore, the system (52) can be rewritten as ρk = O ( ρ2k + |z′|2 + ε |ln ε| ) , k ∈ {1, . . . , N}, λk = O ( |z′|2 + (if n ⩾ 5) ε+ (if n = 4) |ε ln ε| ) , D2f1 (bk) = O ( ε+ |z′|2 ) , k ∈ {1, . . . , N}. (53) Since, D2f1 (bk) is non-degenerate, using Brouwer’s fixed point theorem, we deduce that (53) has a solution ( ρε, λε, (z′) ε) . Furthermore, it holds, for each k ∈ {1, . . . , N}, ρεk = O (ε |ln ε|) ; λεk = O ( (if n ⩾ 5) ε+ (if n = 4) |ε ln ε| ) ; (z′k) ε = O(ε). Taking αε k, µ ε k and aεk by using the equations (47), (48) and (49) and taking uε = ∑N k=1 α ε kwaε k,µ ε k + v̄ε, we deduce that uε is a critical point of Iε and therefore it satisfies{ (−∆+ V )uε = f |uε| 4 n−2−εuε in Ω, ∂uε/∂ν = 0 on ∂Ω. Finally, we have to prove that uε > 0. To this aim, let u−ε := max(0,−uε), it follows that 0 ≤ u−ε ≤ |vε|. Furthermore, multiplying the previous equation by u−ε and integrating over Ω, we obtain ∥u−ε ∥2 = ∫ Ω ∇uε∇u−ε + ∫ Ω V uεu − ε R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 18 of 31 = ∫ Ω (−∆+ V )uεu − ε = ∫ Ω f |uε|p−ε−1uεu − ε = ∫ Ω f(u−ε ) p−ε+1 (54) which implies that ∥u−ε ∥2 = o(1) (since 0 ≤ u−ε ≤ |vε|). Now, using the Holder’s inequality, we obtain ∥u−ε ∥2 ≤ c∥u−ε ∥p+1−ε. Thus, u−ε has to be zero and therefore, by the maximum principle, we derive that uε > 0 in Ω. Hence uε is a solution of Problem (Pε). This completes the proof of Theorem 1. 6. Conclusion By expanding the gradient of the associated functional and testing the equation with appro- priate vector fields, we were able to construct boundary blow-up solutions for the problem (Pε), which exhibit isolated bubbles. This construction exploits the structure of the problem, using asymptotic analysis to capture the intricate behavior of the concentration points and the cor- responding blow-up rates of the solution as the perturbation parameter ε approaches zero. By carefully analyzing the interaction between the nonlinearities of the equation and the boundary conditions, we establish a connection between the number of isolated bubbles and the topology of the problem. This approach ultimately leads to a multiplicity result, demonstrating that the number of boundary blow-up solutions is closely related to the number of non-degenerate critical points of the restriction of the function f on the boundary of the domain Ω. This result provides a deeper understanding of the solution structure, offering insights into bifurcation behavior and the stability of solutions as the boundary conditions are varied. Nevertheless, several promising avenues for further research and open questions remain: (i) Do boundary clustered bubble solutions exist for the problem? (ii) Can we provide a complete description of the asymptotic profile of the boundary blowing up solutions? (iii) What occurs if the critical points of the restriction f1 of the function f on the boundary are degenerate? In particular, what occurs when f1 satisfies certain flatness conditions? (iv) Is it possible to get the same results presented in this paper when the solutions do not converge weakly to zero? 7. Appendix In this section, we gather estimates for several integrals, which are crucial for refining the expansion of the gradient of the Euler-Lagrange functional Jε. Additionally, we prove the coercivity of the quadratic form defined by (10). 7.1. Useful estimates of some integrals We start by the following lemma which is extracted from [27] (see equations (D.6), (D.7) and (D.8)). Lemma 4. [27] Let n ⩾ 3, a ∈ ∂Ω and µ be a large real. We have (i) ∫ Ω |∇ωa,µ|2 = Sn − c1 H(a) µ +O ( 1 µ2 ) , R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 19 of 31 (ii) ∫ Ω ∇ωa,µ∇ ( µ ∂ωa,µ ∂µ ) = c1 2 H(a) µ +O ( 1 µ2 ) , (iii) ∫ Ω ∇ωa,µ∇ ( 1 µ ∂ωa,µ ∂τj ) = O ( 1 µ2 ) ∀j ∈ {1, . . . , n− 1}, where τ ′js, for j = 1, . . . , n − 1, build an orthonormal system of coordinates on the tangent space to ∂Ω at the point a ∈ ∂Ω, the constant Sn is defined in (22) and the constant c1 is defined by c1 := [n(n− 2)] (n−2)/2 (n− 2)2 4 meas ( Sn−2 ) Γ (n+3 2 ) Γ ( n−3 2 ) Γ(n) . (55) We notice that, in this paper we use ωa,µ = β0Ua,µ where β0 = [n(n − 2)](n−2)/4 and Ua,µ is the function used in [27]. For this reason, there is some changes in the constants found in Lemma 4 and the following lemmas with the corresponding results in [27]. The second lemma deals with some integrals involving the bubbles. Lemma 5. Let n ⩾ 4, a ∈ ∂Ω and µ be a large real. It holds: (i) ∫ Ω ω2 a,µ ⩽ c { µ−2 if n ⩾ 5, µ−2 lnµ if n = 4, (ii) ∣∣∣∣∫ Ω ωa,µµ ∂ωa,µ ∂µ ∣∣∣∣ ⩽ c { µ−2 if n ⩾ 5, µ−2 lnµ if n = 4, (iii) ∣∣∣∣∫ Ω ωa,µ 1 µ ∂ωa,µ ∂a ∣∣∣∣ ⩽ c µ3 . Proof. Notice that, since Ω is bounded, there exists R > 0 such that Ω ⊂ B(a,R). Claim (i) follows by standard computations. Concerning Claim (ii), it follows from the first one and the fact that µ ∣∣∣∂ωa,µ ∂µ ∣∣∣ ⩽ cωa,µ. Finally, for Claim (iii), observe that 1 µ ∣∣∣∣∂ωa,µ ∂a ∣∣∣∣ ⩽ 1 µ|x− a| ωa,µ. Hence, the result follows by standard computations. The next lemma is extracted from [27] (see the equations (D.17), (D.18) and (D.19)). Lemma 6. [27] Let n ⩾ 4, a ∈ ∂Ω and µ be a large real. There hold: (i) ∫ Ω ω 2n n−2 a,µ = Sn − 2n n− 2 c4 H(a) µ +O ( 1 µ2 ) , (ii) ∫ Ω ω n+2 n−2 a,µ µ ∂ωa,µ ∂µ = c4 H(a) µ +O ( 1 µ2 ) , (iii) ∫ Ω ω n+2 n−2 a,µ 1 µ ∂ωa,µ ∂τj = O ( 1 µ2 ) ∀j ∈ {1, . . . , n− 1}, where c4 := n− 2 2n [n(n− 2)] n/2 1 4 meas ( Sn−2 ) Γ (n+1 2 ) Γ ( n−1 2 ) Γ(n) . (56) R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 20 of 31 We also have the following estimates: Lemma 7. Let a ∈ ∂Ω and µ be a large real. It holds: (i) µ ∣∣∣∣∂2ωa,µ ∂ν∂µ ∣∣∣∣ ⩽ c ∣∣∣∣∂ωa,µ ∂ν ∣∣∣∣, (ii) (∫ ∂Ω ∣∣∣∂ωa,µ ∂ν ∣∣∣(2n−2)/n)n/(2n−2) ≤ c µ , (iii) (∫ ∂Ω ∣∣∣∂2ωa,µ ∂ν∂a ∣∣∣(2n−2)/n)n/(2n−2) ⩽ c. Proof. Claims (ii) and (iii) are extracted from [27] (See the equations (D.49) and (D.50)). Concerning Claim (i), it follows easily. We end this subsection by the following two lemmas: Lemma 8. Let a ∈ ∂Ω and µ be a large real. It holds: (i) ∫ Ω (x− a) · τjω n+2 n−2 a,µ µ ∂ωa,µ ∂µ = O ( 1 µ2 ) ∀j ∈ {1, . . . , n− 1}, (ii) ∫ Ω (x− a) · νaω n+2 n−2 a,µ µ ∂ωa,µ ∂µ = c5 µ +O ( 1 µ2 ) , (iii) ∫ Ω ω n+2 n−2 a,µ µ ∂ωa,µ ∂µ ln ( 1 + µ2 |x− a|2 ) = −c6 2 +O ( 1 µ ) , (iv) ∫ Ω (x− a) · τkω n+2 n−2 a,µ 1 µ ∂ωa,µ ∂τj = { O ( µ−2 ) if k ̸= j, c7 µ +O ( 1 µ2 ) if k = j, for each j ∈ {1, . . . , n− 1} and k ∈ {1, . . . , n}, (v) ∫ Ω ω n+2 n−2 a,µ 1 µ ∂ωa,µ ∂τj ln ( 1 + µ2 |x− a|2 ) = O ( 1 µ2 ) ∀j ∈ {1, . . . , n− 1}, where c5 := [n(n− 2] n/2 n− 2 2 ∫ Rn + xn |x|2 − 1 (1 + |x|2)n+1 dx > 0, (57) c6 := n− 2 2 [n(n− 2)]n/2 ∫ Rn |x|2 − 1 (1 + |x|2)n+1 ln ( 1 + |x|2 ) dx > 0, (58) c7 := n− 2 2n [n(n− 2)]n/2 ∫ Rn |x|2 (1 + |x|2)n+1 dx. (59) Proof. Without loss of generality, we can assume that a = 0 and νa = −en. (60) Since we assumed that Ω is smooth, there exit ρ > 0 (we take it small) and a function φ : Bn−1(0, ρ) ⊂ Rn−1 −→ R such that φ(0) = 0, φ′(0) = 0 and Ω ∩Bn(0, ρ) = {x = (x′, xn) ∈ Bn−1(0, ρ)× R : xn > φ (x′)} . R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 21 of 31 Since φ′(0) = 0, it is easy to see that φ (x′) = O ( |x′|2 ) ∀x′ ∈ Bn−1(0, ρ). (61) Observe that ∫ Ω\B(0,ρ) |x|ω 2n n−2 0,µ ⩽ ∫ Ω\B(0,ρ) |x| µn|x|2n dx ⩽ c µn . (62) To estimate the integral over Ω ∩B(0, ρ), we introduce the following sets B+(0, ρ) := {x = (x′, xn) ∈ B(0, ρ) : xn > 0} , Ω1 := {x = (x′, xn) ∈ B(0, ρ) : 0 < xn < φ (x′)} , Ω2 := {x = (x′, xn) ∈ B(0, ρ) : φ (x′) < xn < 0} , and we have ∫ Ω∩B(0,ρ) . . . = ∫ B+(0,ρ) . . .− ∫ Ω1 . . .+ ∫ Ω2 . . . . (63) Proof of (i): Let j ∈ {1, · · · , n− 1}. By (60), it follows that (x− a) · τj = xj . By oddness of the function, it is easy to get that the first integral is zero. Concerning the other ones, using (61), we derive that |xn| ⩽ φ (x′) = O ( |x′|2 ) ∀ (x′, xn) ∈ Ωi, i = 1, 2. (64) Furthermore, it is easy to see that 1 + µ2|x|2 ⩾ 1 + µ2 |x′|2. Thus, we obtain, for i ∈ {1, 2},∣∣∣∣∫ Ωi . . . ∣∣∣∣ ⩽ ∫ Ωi |x′|ω 2n n−2 a,µ ⩽ c ∫ Ωi µn |x′| (1 + µ2|x′|2 dx′dxn ⩽ c ∫ Bn−1(0,ρ) µn |x′|3( 1 + µ2 |x′|2 )n dx′ ⩽ c µ2 · (65) Hence, Eqs. (62), (63) and (65) end the proof of Claim (i). Proof of (ii): From (60), we deduce that (x− a) · νa = −xn. As in the proof of Claim (i), for i ∈ {1, 2}, we have (using (64))∣∣∣∣∫ Ωi xnω n+2 n−2 a,µ µ ∂ωa,µ ∂µ ∣∣∣∣ ⩽ c ∫ Ωi µn|xn|( 1 + µ2 |x′|2 )n dx′dxn ⩽ c ∫ Bn−1(0,ρ) µn |x′|4( 1 + µ2 |x′|2 )n dx′ ⩽ c µ3 . (66) For the integral over B+(0, ρ), it holds∫ B+(0,ρ) −xnω n+2 n−2 0,µ µ ∂ω0,µ ∂µ = ∫ B+(0,ρ) −xn ( n− 2 2 ) ω 2n n−2 0,µ 1− µ2|x|2 1 + µ2|x|2 dx = −β 2n n−2 0 n− 2 2 ∫ B+(0,ρ) µnxn 1− µ2|x|2 (1 + µ2|x|2)n+1 dx = −β 2n n−2 0 n− 2 2 1 µ ∫ Rn + xn 1− |x|2 (1 + |x|2)n+1 dx+O ( 1 µn ) . (67) Combining Eqs. (62), (66) and (67), the proof of Claim (ii) follows. Proof of (iii): Following the proof of the previous claims, we need to estimate:∣∣∣∣∣ ∫ Ω\B(0,ρ) . . . ∣∣∣∣∣ ⩽ ∫ Ω\B(0,ρ) ω 2n n−2 a,µ ln ( 1 + µ2 |x− a|2 ) ⩽ c lnµ µn , (68) R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 22 of 31∫ B+(0,ρ) . . . = β 2n n−2 0 n− 2 2 ∫ B+(0,ρ) µn ( 1− µ2|x|2 ) (1 + µ2|x|2)n+1 ln ( 1 + µ2|x|2 ) dx = β 2n n−2 0 n− 2 2 ∫ (B+(0,λρ) 1− |x|2 (1 + |x|2)n+1 ln ( 1 + |x|2 ) dx = −1 2 β 2n n−2 0 n− 2 2 ∫ Rn |x|2 − 1 (1 + |x|2)n+1 ln ( 1 + |x|2 ) dx+O ( lnµ µn ) . (69) For the integrals over Ωi, i = 1, 2, note that, using Eq. (64), we have |xn| ≤ c |x′|2, which implies that 1 + µ2|x|2 = 1 + µ2 |x′|2 + µ2x2n ⩽ 1 + µ2 |x′|2 ( 1 + c |x′|2 ) ⩽ 2(1 + µ2 |x′|2), (70) since |x′| < ρ which is small. Thus we obtain ∣∣∣∣∫ Ωi . . . ∣∣∣∣ ⩽ c ∫ Ωi µn ln ( 1 + µ2|x|2 ) (1 + µ2|x|2)n dx ≤ c ∫ Bn−1(0,ρ) µn |x′|2 ln ( 1 + µ2 |x′|2 ) ( 1 + µ2 |x′|2 )n dx′ ⩽ c µ ∫ Rn−1 |x′|2 ln ( 1 + |x′|2 ) ( 1 + |x′|2 )n dx′ ⩽ c µ . (71) Hence, (68), (69) and (71) imply the proof of Cham (iii). Proof of (iv): Note that, by (60), it follows that (x− a) · τk = xk and ∂ωa,µ ∂τj = ∂ωa,µ ∂aj . As before, we compute: ∣∣∣∣∣ ∫ Ω\B(0,ρ) . . . ∣∣∣∣∣ ⩽ ∫ Ω\B(0,ρ) |x| 1 µ|x| ω 2n n−2 0,µ ⩽ c µn+1 , (72) where we have uses the fact that ∣∣∣∂ωa,µ ∂a ∣∣∣ ⩽ c ωa,µ |x−a| . Concerning the integral over Ω ∩B(0, ρ), using (63), we need to compute:∫ B+(0,ρ) · · · = (n− 2)β 2n n−2 0 ∫ B+(0,ρ) xk µn+1xj (1 + µ2|x|2)n+1 dx = 0, k ̸= j, (73) (by oddness with respect to the variable xj). However, if k = j, we obtain∫ B+(0,ρ) · · · = (n− 2)β 2n n−2 0 ∫ B+(0,ρ) µn+1x2j (1 + µ2|x|2)n+1 dx = 1 2 n− 2 µ β 2n n−2 ∫ B(0,µρ) x2j (1 + |x|2)n+1 dx = 1 2µ n− 2 n β 2n n−2 0 ∫ B(0,µρ) |x|2 (1 + |x|2)n+1 dx = c7 µ +O ( 1 µn+1 ) . (74) It remains the integrals over Ωi, i = 1, 2. Using (70), it holds∣∣∣∣∫ Ωi . . . ∣∣∣∣ ⩽ c ∫ Ωi |xk| 1 µ|x| ω 2n n−2 a,µ ⩽ c ∫ Ωi µn−1 (1 + µ2|x′|2)n dx′dxn R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 23 of 31 ⩽ c ∫ Bn−1(0,ρ) µn−1 |x′|2( 1 + µ2 |x′|2 )n ⩽ c µ2 . (75) Thus, Combining (72) - (75), the proof of Claim (iv) follows. Proof of (v): It can be done in the same way than the proof of Claims (iii) and (iv). Hence, we omit it. Lemma 9. Let a1, a2 ∈ ∂Ω with |a1 − a2| ⩾ c > 0 and µ1, µ2 be large reals. We have: (i) ∫ Ω |∇ωa1,µ1 | |∇ωa2,µ2 | ⩽ c (µ1µ2)(n−2)/2 ⩽ c µn−2 1 + c µn−2 2 , (ii) ∫ Ω |∇ωa1,µ1 | ∣∣∣∣∇(µ2 ∂ωa2,µ2 ∂µ2 )∣∣∣∣ ⩽ c (µ1µ2)(n−2)/2 ⩽ c µn−2 1 + c µn−2 2 , (iii) ∫ Ω |∇ωa1,µ1 | ∣∣∣∣∇( 1 µ2 ∂ωa2,µ2 ∂a2 )∣∣∣∣ ⩽ c µ (n−2)/2 1 lnµ2 µ n/2 2 , (iv) ∫ Ω ωa1,µ1ωa2,µ2 ⩽ c (µ1µ2)(n−2)/2 ⩽ c µn−2 1 + c µn−2 2 , (v) ∫ Ω ωa1,µ1 ∣∣∣∣µ2 ∂ωa2,µ2 ∂µ2 ∣∣∣∣ ⩽ c (µ1µ2)(n−2)/2 ⩽ c µn−2 1 + c µn−2 2 , (vi) ∫ Ω ωa1,µ1 ∣∣∣∣ 1µ2 ∂ωa2,µ2 ∂a2 ∣∣∣∣ ⩽ c µ (n−2)/2 1 µ n/2 2 ⩽ c µn−1 1 + c µn−1 2 , (vii) ∫ Ω ω n+2 n−2 a1,µ1ωa2,µ2 ⩽ c (µ1µ2)(n−2)/2 ⩽ c µn−2 1 + c µn−2 2 , (viii) ∫ Ω ω n+2 n−2 a1,µ1µ2 ∣∣∣∣∂ωa2,µ2 ∂µ2 ∣∣∣∣ ⩽ c (µ1µ2)(n−2)/2 ⩽ c µn−2 1 + c µn−2 2 , (ix) ∫ Ω ω n+2 n−2 a1,µ1 1 µ2 ∣∣∣∣∂ωa2,µ2 ∂a2 ∣∣∣∣ ⩽ c µ (n−2)/2 1 µ n/2 2 ⩽ c µn−1 1 + c µn−1 2 . Proof. We will focus on the proof of the first one and the other proofs can be done in the same way. Note that ∣∣∇ωai,µi| ∣∣ ⩽ c µ n+2 2 |x− ai|( 1 + µ2 i |x− ai|2 )n/2 ⩽ c µ (n−2)/2 i |x− ai|n−1 . Thus, let ρ := |a1 − a2| /2, it holds :∫ Ω |∇ωa1,µ1 | |∇ωa2,µ2 | ⩽ 1 µ (n−2)/2 1 µ (n−2)/2 2 ∑ i=1,2 ∫ B(ai,ρ) dx |x− ai|n−1 + ∫ Ω\∪B(ai,ρ) 1dx  . ⩽ c (µ1µ2) (n−2)/2 ⩽ c ( 1 µn−2 1 + 1 µn−2 2 ) . Hence, the proof of Claim (i) is completed. R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 24 of 31 7.2. Coercivity of the quadratic form The goal of this subsection is to prove Proposition 1. To this aim, for µ > 0 and x = (x1, . . . , xn) ∈ Rn, we denote by ψ1(x) := ω0,µ(x) := β0 µ(n−2)/2 (1 + µ2|x|2)(n−2)/2 , ψ2(x) := µ ∂ω0,µ ∂µ (x) = n− 2 2 β0 µ(n−2)/2 ( 1− µ2|x|2 ) (1 + µ2|x|2)n/2 , ψj(x) := (n− 2)β0 µn/2xj−2 (1 + µ2|x|2)n/2 , for j ∈ {3, · · · , n+ 2}. (76) We begin by the following lemma: Lemma 10. Let ρ > 0 be a small radius and B+ ρ := { x := (x′, xn) ∈ Rn−1 × R : |x| < ρ and xn > 0 } . For µ large and γ̄ > 0, let us define Q+(v) := ∫ B+ ρ |∇v|2 + γ̄ ∫ B+ ρ v2 − n+ 2 n− 2 ∫ B+ ρ ω 4 n−2 0,µ v2 Then there exists a constant β1 > 0 such that Q+(v) ⩾ β1 (∫ B+ ρ |∇v|2 + γ̄ ∫ B+ ρ v2 ) ∀v ∈ E+ µ , where E+ µ := { v ∈ H1 ( B+ ρ ) : ∫ B+ ρ ∇v · ∇ψj = 0 ∀j ∈ {1, . . . , n+ 1} } . Proof. Let us introduce the function ṽ defined on B(0, ρ) by for y := (y′, yn) ∈ B(0, ρ), ṽ(y) := { v(y) if yn > 0, v (y′,−yn) if yn < 0. Easy Computations imply that ṽ ∈ H1 (B(0, ρ)) , 2Q+(v) = Q̃+(ṽ) := ∫ B(0,ρ) |∇ṽ|2 + γ̄ ∫ B(0,ρ) (ṽ)2 − n+ 2 n− 2 ∫ B(0,ρ) ω 4 n−2 0,µ ṽ2. (77) Notice that the function Q̃+ is a positive definite quadratic form on the space E0,µ := { v ∈ H1 (B(0, ρ)) : ∫ B(0,ρ) ∇v∇ψj = 0 ∀j = 1, . . . , n+ 2 } , (See Proposition 1 of [32] and equation (19) by taking Ω = B (0, 1), K = γ̄ and N = 1). This implies that there exists a constant β0 > 0 such that Q̃+(w) ⩾ β0∥w∥H1(B(0,ρ)) ∀w ∈ E0,µ. (78) R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 25 of 31 In the following, we will prove that ṽ ∈ E0,µ. For this aim, for 1 ⩽ j ⩽ n+ 1, we compute∫ B(0,ρ) ∇ṽ∇ψj = 2 ∫ B+ ρ ∇v∇ψj = 0, since v ∈ E+ µ . Now, for j = n+ 2, observe that (by easy computations) −∆ψn+2 = n+ 2 n− 2 ω0,µψn+2 in B(0, ρ) ; ∂ψn+2 ∂ν = c(ρ, µ)xn on ∂B(0, ρ). Thus, by oddness (with respect the variable xn), we obtain∫ B(0,ρ) ∇ṽ∇ψn+2 = ∫ B(0,ρ) −∆ψn+2ṽ + ∫ ∂B(0,ρ) ∂ψn+2 ∂ν ṽ = 0, Hence, v ∈ E0,µ and the assumptions of Proposition 1 of [32] are satisfied. Combining (77) and (78) (by taking w = ṽ), we get Q+(v) ⩾ (β0/2) ∥ṽ∥2H1(B(0,ρ)) ⩾ β0∥v∥2H1(B+ ρ ) . We remark that∫ B+ ρ |∇v|2 + γ̄ ∫ B+ ρ v2 ⩽ ∥v∥2 H1(B+ ρ ) ⩽ 1 γ̄ (∫ B+ ρ |∇v|2 + γ̄ ∫ B+ ρ v2 ) if γ̄ ⩽ 1, 1 γ̄ (∫ B+ ρ |∇v|2 + γ̃ ∫ B+ ρ v2 ) ⩽ ∥v∥2 H1(B+ ρ ) ⩽ ∫ B+ ρ |∇v|2 + γ̄ ∫ B+ ρ v2 if γ̄ > 1. (79) The proof of the lemma is thereby completed. Notice that, for a ∈ ∂Ω, a neighborhood of a in Ω is not necessary a half ball. For this reason, we need to take a general case. Lemma 11. Let a ∈ ∂Ω, µ be a large real and ρ be a small radius. Let Qa,ρ(v) := ∫ B(a,ρ)∩Ω |∇v|2 + γ̄ ∫ B(a,ρ)∩Ω v2 − n+ 2 n− 2 ∫ B(a,ρ)∩Ω ω 4 n−2 a,µ v2. Then, there exists a constant β2 > 0 such that Qa,ρ(v) ⩾ β2 (∫ B(a,ρ)∩Ω |∇v|2 + γ̄ ∫ B(a,ρ)∩Ω v2 ) + o ( ∥v∥2H1(Ω) ) ∀v ∈ Fa,µ, where Fa,µ is defined in (4). Proof. Let (e1, . . . , en) be the canonical basis of Rn. Without loss of generality, we can assume that a = 0 and νa = −en (which implies that the tangent space to ∂Ω at a = 0 is Rn−1 × {0} and a basis of this tangent space is (e1, . . . , en−1)). Since ρ is small and Ω is a regular domain, there exists a smooth function f : Rn−1 −→ R, satisfying f(0) = 0, ∇f(0) = 0 and Ω ∩B(0, ρ) = { x := (x′, xn) ∈ Rn−1 × R : |x| < ρ, xn > f (x′) } . R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 26 of 31 Now, we define φ : Ω ∩B(0, ρ) −→ Rn−1 × R, φ (x′, xn) = (x′, xn − f (x′)) . (80) From (80), we remark that there exists a neighborhood V of 0 in B (0, ρ) such that φ induces a diffeomorphism between V ∩ Ω and B+ := {x = (x′, xn) ∈ B (0, ρ/2) : xn > 0}, that is φ(V ∩ Ω) = B+. (81) In addition, we have B (0, ρ/4) ⊂ V. Furthermore, from the definition of φ in (80), we deduce that ∂φ ∂xi (x) = ei − ∂f ∂xi (x′) en for 1 ⩽ i ⩽ n− 1 and ∂φ ∂xn (x) = en, (82) which implies that the Jacobian of φ at each point x is 1 (|Jacφ| = 1). Now, let us define the function v1 by v1 : B+ −→ R, v1 := v ◦ φ−1. (83) Using (82), easy computations imply that |∇v(x)|2 = |(∇v1) (φ(x))|2 +O ( ρ |(∇v1) (φ(x))|2 ) , which implies that, by using (81),∫ V∩Ω |∇v(x)|2 dx = ∫ B+ |∇v1(z)|2 dz +O ( ρ ∫ B+ |∇v1(z)|2 dz ) , and therefore∫ V∩Ω |∇v|2 + γ̄ ∫ V∩Ω |v|2 = ∫ B+ |∇v1|2 + γ̄ ∫ B+ (v1) 2 +O ( ρ ∫ B+ |∇v1|2 ) . (84) Concerning the last integral in the definition of Qa,ρ, we have∫ V∩Ω ω 4 n−2 0,µ (x)v2(x)dx = ∫ V∩Ω ω 4 n−2 a,µ (x)v21(φ(x))dx = ∫ B+ ω 4 n−2 0,µ ( φ−1(z) ) v21(z)dz. (85) Observe that equation (C. 25) of [27] gives us( 1 + µ2|x|2 )γ = ( 1 + µ2 ∣∣φ−1(z) ∣∣2)γ = ( 1 + µ2|z|2 )γ +O (( 1 + µ2|z|2 )γ−1 µ2|z|2ρ ) = ( 1 + µ2|z|2 )γ +O (( 1 + µ2|z|2 )γ ρ ) . (86) Hence we obtain∫ B+ ω 4 n−2 0,µ ( φ−1(z) ) v21(z)dz = ∫ B+ ω 4 n−2 0,µ v21(z)dz +O ( ρ ∫ B+ ω 4 n−2 0,µ v21(z)dz ) = ∫ B+ ω 4 n−2 0,µ v21(z)dz +O ( ρ ∥v1∥2L2n/(n−2)(B+) ) . (87) Combining (84), (85) and (87), we get Qa,ρ(v) = ∫ (B(0,ρ)∩Ω)\V ( |∇v|2 + γ̄v2 ) − n+ 2 n− 2 ∫ (B(0,ρ)∩Ω)\V ω 4 n−2 0,µ v2 R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 27 of 31 + ∫ B+ |∇v1|2 + γ̄ ∫ B+ (v1) 2 − n+ 2 n− 2 ∫ B+ ω 4 n−2 0,µ v21 +O ( ρ ∥v1∥H1(B+) ) , (88) Observe that, since B(0, ρ/4) ⊂ V, we deduce that∫ (B(0,ρ)∩Ω\V) ω 4 n−2 0,µ v2 ⩽ c∥v∥2L2(B(0,ρ)) (∫ Rn\B(0,ρ/4) ω 2n n−2 0,µ )2/n ⩽ c (µρ)2 ∥v∥2L2(B(0,R)∩Ω). Thus (88) becomes Qa,ρ(v) = ∫ (B(0,ρ)∩Ω)\V ( |∇v|2 + γ̄v2 ) +Q+ (v1) +O ( ρ∥v∥H1(B(0,ρ)∩Ω) ) . (89) At this step, we need to apply Lemma 10 to the function v1, defined by (83), but v1 /∈ E+ µ . For this reason, we decompose v1 as follows: v1 = n+1∑ j=1 σjψj + v⊥1 with v⊥1 ∈ E+ µ , where the ψj ’s are defined in (76). Now, we need to estimate the parameters σj ’s. Observe that, on one hand we have:∫ B+ ∇v1∇ψ1 = σ1 ∫ B+ |∇ψ1|2 + ∑ j ̸=1 ∫ B+ ∇ψ1∇ψj = cσ1 + o (∑ |σj | ) . (90) On the other hand, we have: ∫ B+ ∇v1∇ψ1 = ∫ B+ ∇v1∇ω0,µ = ∫ B+ (−∆ω0,µ) v1 + ∫ ∂B+ ( ∂ ∂ν ω0,µ ) v1 = ∫ B+ ω n+2 n−2 0,µ v1 + ∫ ∂B+ ( ∂ ∂ν ω0,µ ) v1. (91) Let Γ1 := {x = (x′, xn) : |x| = ρ and xn ⩾ 0} and Γ2 := {x = (x, 0) : |x| ⩽ ρ}. It is easy to see that ∂B+ = Γ1UΓ2, ∂ ∂ν ω0,µ = 0 on Γ2 and ∂ ∂ν ω0,µ = O ( 1 µ(n−2)/2 ) on Γ1. Thus ∣∣∣∣∫ ∂B+ ( ∂ ∂ν ω0,µ ) v1 ∣∣∣∣ ⩽ c µ(n−2)/2 ∫ Γ1 |v1| ⩽ c µ(n−2)/2 ∥v1∥H1(B+) . (92) For the other integral, we get∫ B+ ω n+2 n−2 0,µ v1 = ∫ B+ ω n+2 n−2 0,µ (z)v ( φ−1(z) ) dz = ∫ V ∩Ω ω n+2 n−2 0,µ (φ(x))v(x)dx. (93) Using (86), (93) becomes∫ B+ ω n+2 n−2 0,µ v1 = ∫ V∩Ω ω n+2 n−2 0,µ (x)v(x)dx+O ( ρ ∫ V∩Ω ω n+2 n−2 0,µ (x) |v(x)| dx ) = ∫ Ω ω n+2 n−2 0,µ v − ∫ Ω\V ω n+2 n−2 0,µ v +O ( ρ∥v∥L2n/n−2(V ∩Ω) ) = ∫ Ω ∇ω0,µ∇v − ∫ ∂Ω ( ∂ ∂ν ω0,µ ) v +O ( ∥v∥L2n/n−2(Ω) [ ρ+ 1 (µρ)(n+2)/2 ]) . R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 28 of 31 Using the fact that v ∈ Fa,µ and equation (19), we obtain∫ B+ ω n+2 n−2 0,µ v1 = O ( ∥v∥L2n/n−2(Ω) [ 1 µ + ρ+ 1 (µρ)(n+2)/2 ]) = o ( ∥v∥L2n/(n−2)(Ω) ) . (94) Combining (90), (91), (92) and (94), we get σ1 = o (∑ |σj | ) + o ( ∥v∥H1(Ω) ) . In the same way, we get the estimate of σi for i ⩾ 2 and therefore we obtain σi = o (∑ |σj | ) + o ( ∥v∥H1(Ω) ) ∀i = 1, . . . , n+ 1, which implies that σi = o ( ∥v∥H1(Ω) ) ∀i = 1, . . . , n+ 1. Hence we deduce that v1 − v⊥1 = o (∥v∥ω0,µ) and ∇ ( v1 − v⊥1 ) = O ( ∥v∥ ∑ |∇ψj | ) . (95) This implies that (by using v⊥2 ∈ E+ µ ) Q+ (v1) = ∫ B+ |∇v1|2 + γ̄ ∫ B+ v21 − n+ 2 n− 2 ∫ B+ ω 4 n−2 0,µ v21 = ∫ B+ ∣∣∇v⊥1 ∣∣2 + ∫ B+ ∣∣∇ (v1 − v⊥1 )∣∣2 + γ̄ ∫ B+ ( v⊥1 )2 + 2γ̄ ∫ B+ ( v⊥1 ) ( v1 − v⊥1 ) + γ̄ ∫ B+ ( v1 − v⊥1 )2 − n+ 2 n− 2 ∫ B+ ω 4 n−2 0,µ [( v⊥1 )2 + 2v⊥1 ( v1 − v⊥1 ) + ( v1 − v⊥1 )2] = Q+ ( v⊥1 ) + o ( ∥v⊥1 ∥2 + ∥v∥2 ) ⩾ 1 2 β1 (∫ B+ ∣∣∇v⊥1 ∣∣2 + γ̄ ∫ B+ ( v⊥1 )2) + o ( ∥v∥2 ) , (96) by using Lemma 10, Eq. (79) and the fact that v⊥1 ∈ E+ µ . Combining (96), (95), (89) and (84), the proof of Lemma 11 follows. Now, we are ready to prove Proposition 1. Proof of Proposition 1 Let ρ be a small radius and let Bi := B (ai, ρ)∩Ω. Since |ai − aj | ⩾ c > 0 for i ̸= j, it follows that Bi ∩Bj = ∅ for each i ̸= j. Thus we get Q(v) = q∑ i=1 (∫ Bi |∇v|2 + ∫ Bi V v2 − n+ 2 n− 2 ∫ Bi ω 4 n−2 ai,µiv 2 ) + ∫ Ω\(UBi) |∇v|2 + ∫ Ω\(UBi) V v2 − n+ 2 n− 2 q∑ i=1 ∫ Ω\Bi ω 4 n−2 ai,µiv 2. Observe that, for each i ∈ {1, . . . , q}, we have ∫ Ω\Bi ω 4 n−2 ai,µiv 2 ⩽ (∫ Ω\Bi v 2n n−2 )n−2 n (∫ Ω\Bi ω 2n n−2 ai,µi )2/n ⩽ c (µiρ) 2 ∥v∥ 2 H1(Ω). (97) R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 29 of 31 In addition, let γ̄ = minV > 0, using (97) and Lemma 11, we derive that Q(v) ⩾ q∑ i=1 Qai,ρ(v) + ∫ Ω\(∪Bi) |∇v|2 + ∫ Ω\(∪Bi) V v2 + ∑ O ( ∥v∥2 (µiρ) 2 ) ⩾ q∑ i=1 β2 (∫ Bi |∇v|2 + γ̄ ∫ Bi v2 ) + ∫ Ω\∪Bi |∇v|2 + ∫ Ω\∪Bi V v2 + ∑ O ( ∥v∥2 (µiρ) 2 ) ⩾ β3∥v∥2, for some positive constant β3 (since ρ is fixed and the µi’s are large). This completes the proof. Acknowledgements The authors gratefully acknowledge Qassim University, represented by the Deanship of Grad- uate Studies and Scientific Research, on the financial support for this research under the number (QU-J-PG-2-2025-53906) during the academic year 1446 AH / 2024 AD. Author Contributions: R.A and M.B.A.: conceptualization, methodology, investigation, writing original draft, writing-review and editing. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by the Deanship of Scientific Research, Qassim University, grant number project QU-J-PG-2-2025-53906. Data Availability Statement: No data to report in this manuscript. Conflicts of Interest: The authors declare no conflict of interest. References [1] E.F. Keller and L.A. Segel. Initiation of slime mold aggregation viewed as an instability. J. Theor. Biol., 26:399–415, 1970. [2] R. Schaaf. Stationary solutions of chemotaxis systems. Trans. Amer. Math. Soc., 292:531–556, 1985. [3] T. Hillen and K.J. Painter. A user’s guide to pde models for chemotaxis. J. Math. Biol., 58:183–217, 2009. [4] M. Bezerra, C. Cuevas, C. Silva, and H. Soto. On the fractional doubly parabolic keller-segel system modelling chemotaxis. Science China Mathematics, 65:1827–1874, 2022. [5] L. Almeida, F. Bubba, B. Perthame, and C. Pouchol. Energy and implicit discretization of the fokker-planck and keller-segel type equations. Networks and Heterogeneous Media, 14:23–41, 2019. [6] C.S. Lin, W.M. Ni, and I. Takagi. Large amplitude stationary solutions to a chemotaxis system. J. Differential Equations, 72:1–27, 1988. [7] W.M. Ni and I. Takagi. On the shape of least-energy solutions to a semi-linear neumann problem. Comm. Pure Appl. Math., 44:819–851, 1991. [8] W.M. Ni and I. Takagi. Locating the peaks of least-energy solutions to a semi-linear neumann problem. Duke Math. J., 70:247–281, 1993. R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 30 of 31 [9] M. Del Pino and P. Felmer. Spike-layered solutions of singularly perturbed elliptic problems in a degenerate setting. Indiana Univ. Math. J., 48:883–898, 1999. [10] E. N. Dancer and S. Yan. Multipeak solutions for a singularly perturbed neumann problem. Pacific J. Math., 189(2):241–262, 1999. [11] M. Del Pino, P. Felmer, and J. Wei. On the role of mean curvature in some singularly perturbed neumann problems. SIAM J. Math. Anal., 31:63–79, 2000. [12] M. Grossi, A. Pistoia, and J. Wei. Existence of multipeak solutions for a semilinear neumann problem via nonsmooth critical point theory. Calc. Var. PDE, 11(2):143–175, 2000. [13] C. Gui and J. Wei. Multiple interior peak solutions for some singularly perturbed neumann problems. J. Differential Equations, 158:1–27, 1999. [14] J. Wei, B. Xu, and W. Yang. On lin-ni’s conjecture in dimensions four and six. Science in China: Mathematics, 49(2):281–306, 2019. [15] Adimurthi and S.L. Yadava. Existence and nonexistence of positive radial solutions of neu- mann problems with critical sobolev exponents. Arch. Rat. Mech. Anal., 115:275–296, 1991. [16] O. Rey and J. Wei. Arbitrary number of positive solutions for elliptic problem with critical nonlinearity. J. Eur. Math. Soc., 7:449–476, 2005. [17] Adimurthi, F. Pacella, and S.L. Yadava. Interaction between the geometry of the boundary and positive solutions of a semilinear neumann problem with critical nonlinearity. J. Funct. Anal., 113:318–350, 1993. [18] W.M. Ni, X.B. Pan, and I. Takagi. Singular behavior of least-energy solutions of a semi-linear neumann problem involving critical sobolev exponents. Duke Math. J., 67:1–20, 1992. [19] L. Caffarelli, B. Gidas, and J. Spruck. Asymptotic symmetry and local behavior of semilinear elliptic equations with critical sobolev growth. Comm. Pure Appl. Math., 42:271–297, 1989. [20] Adimurthi and G. Mancini. The neumann problem for elliptic equations with critical nonlin- earity. In A tribute in honour of G. Prodi, pages 9–25. Scuola Norm. Sup. Pisa, 1991. [21] O. Druet, F. Robert, and J. Wei. The lin-nis problem for mean convex domains. Mem. Amer. Math. Soc., 218(1027), 2012. [22] N. Ghoussoub and C. Gui. Multi-peak solutions for a semilinear neumann problem involving the critical sobolev exponent. Math. Z., 229:443–474, 1998. [23] C. Gui and C.-S. Lin. Estimates for boundary-bubbling solutions to an elliptic neumann problem. J. Reine Angew. Math., 546:201–235, 2002. [24] L. Wang, J. Wei, and S. Yan. A neumann problem with critical exponent in nonconvex domains and lin-ni’s conjecture. Trans. Amer. Math. Soc., 362(9):4581–4615, 2010. [25] Z.-Q. Wang. Construction of multi-peaked solutions for a nonlinear neumann problem with critical exponent in symmetric domains. Nonlinear Anal., 27(11):1281–1306, 1996. [26] J. Wei and S. Yan. Arbitrary many boundary peak solutions for an elliptic neumann problem with critical growth. J. Math. Pures Appl. (9), 88(4):350–378, 2007. [27] O. Rey. Boundary effect for an elliptic neumann problem with critical nonlinearity. Comm. Partial Differential Equations, 22:1055–1139, 1997. [28] O. Rey. The question of interior blow-up points for an elliptic neumann problem: The critical case. J. Math. Pures Appl., 81:655–696, 2002. [29] O. Rey and J. Wei. Blow-up solutions for an elliptic neumann problem with sub- or super- critical nonlinearity, ii: n ≥ 4. Ann. Inst. H. Poincaré, Anal. non-lin, 22(4):459–484, 2005. [30] O. Rey and J. Wei. Blow-up solutions for an elliptic neumann problem with sub- or super- critical nonlinearity, i: n = 3. J. Funct. Anal., 212:472–499, 2004. [31] M. Ben Ayed and K. El Mehdi. Non-existence of interior bubbling solutions for slightly supercritical elliptic problems. Boundary Value Problems, 2023(90), 2023. [32] K. El Mehdi and F. Mohamed Salem. Interior bubbling solutions for an elliptic equation with R. Almushahhin, M. Ben Ayed / Eur. J. Pure Appl. Math, 18 (2) (2025), 6102 31 of 31 slightly subcritical nonlinearity. Mathematics, 11(6):1471, 2023. [33] M. Ben Ayed, K. El Mehdi, and F. Mohamed Salem. Interior multi-peak solution for a slightly subcritical nonlinear neumann equation. Symmetry, 16(291), 2024. [34] M. Struwe. A global compactness result for elliptic boundary value problems involving limiting nonlinearities. Math. Z., 187:511–517, 1984. [35] A. Bahri, Y.Y. Li, and O. Rey. On a variational problem with lack of compactness: the topo- logical effect of the critical points at infinity. Calculus of Variations and Partial Differential Equations, 3:67–94, 1995.