# One-bubble nodal blow-up for asymptotically critical stationary Schrödinger-type equations

Frédéric Robert, Bruno Premoselli

<figure><img src="image_1.png" /></figure>

## To cite this version:

Frédéric Robert, Bruno Premoselli. One-bubble nodal blow-up for asymptotically critical stationary Schrödinger-type equations. 2024. hal-04558683

HAL Id: hal-04558683

https://hal.science/hal-04558683

Preprint submitted on 25 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# One-bubble nodal blow-up for asymptotically critical stationary Schrödinger-type equations

Bruno Premoselli, Frédéric Robert

► To cite this version:

Bruno Premoselli, Frédéric Robert. One-bubble nodal blow-up for asymptotically critical stationary Schrödinger-type equations. Journal of Functional Analysis, 2025, 288 (6), pp.110808. 10.1016/j.jfa.2024.110808. hal-04558683v2

HAL Id: hal-04558683

https://hal.science/hal-04558683v2

Submitted on 8 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ONE-BUBBLE NODAL BLOW-UP FOR ASYMPTOTICALLY CRITICAL STATIONARY SCHRÖDINGER-TYPE EQUATIONS

## BRUNO PREMOSELLI AND FRÉDÉRIC ROBERT

**ABSTRACT.** We investigate in this work families $(u_\epsilon)_{\epsilon>0}$ of sign-changing blowing-up solutions of asymptotically critical stationary nonlinear Schrödinger equations of the following type:

$$ \Delta_g u_\epsilon + h_\epsilon u_\epsilon = |u_\epsilon|^{p_\epsilon - 2} u_\epsilon $$

in a closed manifold $(M, g)$, where $h_\epsilon$ converges in $C^1(M)$. Assuming that $(u_\epsilon)_{\epsilon>0}$ blows-up as a single sign-changing bubble, we obtain necessary conditions for blow-up that constrain the localisation of blow-up points and exhibit a strong interaction between $h$, the geometry of $(M, g)$ and the bubble itself. These conditions are new and are a consequence of the sign-changing nature of $u_\epsilon$.

## 1. INTRODUCTION

### 1.1. Statement of the main results.
Let $(M^n, g)$, $n \ge 3$ be a smooth, connected and closed manifold, where closed means compact without boundary. We study in this paper sign-changing solutions $u \in C^2(M)$ of the equation

$$ (1.1) \qquad \Delta_g u + hu = |u|^{p-2}u \quad \text{in } M $$

where $\Delta_g = -\text{div}_g(\nabla \cdot)$ is the Laplace-Beltrami operator, $h \in C^1(M)$ and $2 < p \le 2^*$ with $2^* = \frac{2n}{n-2}$. When $h \equiv \frac{n-2}{4(n-1)}S_g$, where $S_g$ is the scalar curvature of $(M, g)$, $p = 2^*$ and $u > 0$, (1.1) is the celebrated Yamabe equation. We let $H^1(M)$ be the completion of $C^\infty(M)$ for $u \mapsto \|u\|_{H^1} := \|u\|_2 + \|\nabla u\|_2$. It has been known since the seminal work of Struwe [36] that families of solutions to (1.1) that are uniformly bounded in $H^1(M)$ may develop concentration phenomena in the form of “bubbles” (see also Druet-Hebey-Robert [12] for the case of a compact Riemannian manifold). These “bubbles” correspond to a loss of compactness of the solutions, and understanding the conditions in which they may appear and their mutual interactions has been the subject of many works in the last decades.

Despite the abundance of contributions for positive solutions of (1.1), the literature for sign-changing solutions is less developed. In the present work we take a step in this direction and study sign-changing solutions of (1.1) that blow-up as a single sign-changing bubble. In the sequel, we fix $r_g \in (0, i_g(M))$ where $i_g(M) > 0$ is the injectivity radius of $(M, g)$. Throughout this paper $\chi \in C_c^\infty(\mathbb{R}^n)$ will denote a cutoff function such that $\chi(t) = 1$ if $|t| \le r_g/2$ and $\chi(t) = 0$ if $|t| \ge r_g$. We say that a family $(B_\epsilon)_{\epsilon>0}$ of functions in $H^1(M)$ is a family of bubbles if there exists

<footnote>Date: April 24th, 2024.</footnote>

<footnote>The first author was supported by the FNRS CdR grant J.0135.19, the Fonds Thélam and an ARC Avancé 2020 grant.</footnote>
# ONE-BUBBLE NODAL BLOW-UP FOR ASYMPTOTICALLY CRITICAL STATIONARY SCHRÖDINGER-TYPE EQUATIONS

BRUNO PREMOSELLI AND FRÉDÉRIC ROBERT

ABSTRACT. We investigate in this work families $(u_\varepsilon)_{\varepsilon>0}$ of sign-changing blowing-up solutions of asymptotically critical stationary nonlinear Schrödinger equations of the following type:

$$ \Delta_g u_\varepsilon + h_\varepsilon u_\varepsilon = |u_\varepsilon|^{p_\varepsilon - 2} u_\varepsilon $$

in a closed Riemannian manifold $(M, g)$, where $h_\varepsilon$ converges to $h$ in $C^1(M)$ and $p_\varepsilon \to 2^*$ := $\frac{2n}{n-2}$ as $\varepsilon \to 0$. Assuming that $(u_\varepsilon)_{\varepsilon>0}$ blows-up as a single sign-changing bubble, we obtain necessary conditions for blow-up that constrain the localisation of blow-up points and exhibit a strong interaction between $h$, the geometry of $(M, g)$ and the bubble itself. These conditions are new and are a consequence of the sign-changing nature of $u_\varepsilon$.

## 1. INTRODUCTION

1.1. **Statement of the main results.** Let $(M^n, g)$, $n \ge 3$ be a smooth, connected and closed Riemannian manifold, where closed means compact without boundary. We study in this paper sign-changing solutions $u \in C^2(M)$ of the equation

$$ (1.1) \qquad \Delta_g u + hu = |u|^{p-2}u \quad \text{in } M $$

where $\Delta_g = -\text{div}_g(\nabla\cdot)$ is the Laplace-Beltrami operator, $h \in C^1(M)$ and $2 < p \le 2^*$ with $2^* = \frac{2n}{n-2}$. When $h \equiv \frac{n-2}{4(n-1)}S_g$, where $S_g$ is the scalar curvature of $(M, g)$, $p = 2^*$ and $u > 0$, (1.1) is the celebrated Yamabe equation. We let $H^1(M)$ be the completion of $C^\infty(M)$ for $u \mapsto \|u\|_{H^1} := \|u\|_2 + \|\nabla u\|_2$. It has been known since the seminal work of Struwe [36] that families of solutions to (1.1) that are uniformly bounded in $H^1(M)$ may develop concentration phenomena in the form of “bubbles” (see also Druet-Hebey-Robert [12] for the case of a compact Riemannian manifold). These “bubbles” correspond to a loss of compactness of the solutions, and understanding the conditions in which they may appear and their mutual interactions has been the subject of many works in the last decades.

Despite the abundance of contributions for positive solutions of (1.1), the literature for sign-changing solutions is less developed. In the present work we take a step in this direction and study sign-changing solutions of (1.1) that blow-up as a single sign-changing bubble. In the sequel, we fix $r_g \in (0, i_g(M))$ where $i_g(M) > 0$ is the injectivity radius of $(M, g)$. Throughout this paper $\chi \in C_c^\infty(\mathbb{R}^n)$ will denote a cutoff function such that $\chi(t) = 1$ if $|t| \le r_g/2$ and $\chi(t) = 0$ if $|t| \ge r_g$. We say that a family $(B_\varepsilon)_{\varepsilon>0}$ of functions in $H^1(M)$ is a family of bubbles if there exists

<footnote>Date: December 5th, 2024.</footnote>

<footnote>The first author was supported by the FNRS CdR grant J.0135.19, the Fonds Thélam and an ARC Avancé 2020 grant.</footnote>
a non-zero solution $V$ to

$$ (1.2) \qquad \Delta_{\xi} V = |V|^{2^* - 2} V \text{ in } \mathbb{R}^n, \quad V \in D^{1,2}(\mathbb{R}^n), $$

where $D^{1,2}(\mathbb{R}^n)$ is the completion of $C_c^{\infty}(\mathbb{R}^n)$ for the norm $u \mapsto \|\nabla u\|_2$, and families of points $(x_\varepsilon)_{\varepsilon>0}$ in $M$ and of positive real numbers $(\mu_\varepsilon)_{\varepsilon>0}$ such that

$$ B_{\varepsilon}(x) - \chi(d_g(\cdot, x_{\varepsilon}))\mu_{\varepsilon}^{-\frac{n-2}{2}}V\left(\frac{\exp^{-1}_{x_{\varepsilon}}(\cdot)}{\mu_{\varepsilon}}\right) \to 0 \text{ in } H^1(M) $$

as $\varepsilon \to 0$. Here $\exp_{x_\varepsilon}$ is the exponential chart at $x_\varepsilon$ and $\xi$ is the Euclidean metric, so that $\Delta_\xi = -\sum_{i=1}^n \partial_i^2$. We will say that such a family $(B_\varepsilon)_{\varepsilon>0}$ is centered at $x_\varepsilon$ and of radius $\mu_\varepsilon$ and modeled on $V$, and $V$ itself will sometimes be referred to as the bubble. A more precise definition is given in Definition 3.1 and properties of such families are investigated in Section 3 below. All the solutions considered in this paper are allowed to change sign unless mentioned otherwise.

Our main results establish a necessary condition for the existence of blowing-up solutions of (1.1) that blow-up with a single sign-changing bubble. We introduce one last notation as follows: if $B = (B_\varepsilon)_{\varepsilon>0}$ is a family of bubbles centered at $x_\varepsilon$, of radius $\mu_\varepsilon$ and modeled on $V$, we define

$$ \begin{aligned} \text{Weyl}_g \otimes B := & \frac{4n}{3(n-2)^2} \left( \int_{\mathbb{R}^n} |V|^{2^*} dx \right)^{-1} \times \\ & \text{Weyl}_g(x_0)_{i\alpha j\beta} \int_{\mathbb{R}^n} x^\alpha x^\beta \partial_{ij}^2 V \left( \frac{n-2}{2} V + x^l \partial_l V \right) dx \end{aligned} $$

where $x_0 = \lim_{\varepsilon \to 0} x_\varepsilon$. We first state a result when $n \ge 5$.

**Theorem 1.1.** Let $(M^n, g), n \ge 5$ be a smooth, connected and closed manifold, let $h, (h_\varepsilon)_{0<\varepsilon\le 1} \in C^1(M)$ such that $\lim_{\varepsilon\to 0} h_\varepsilon = h$ in $C^1(M)$ and $\Delta_g + h$ is coercive. Let $(p_\varepsilon)_{0<\varepsilon\le 1} \in (2, 2^*]$ be such that $\lim_{\varepsilon\to 0} p_\varepsilon = 2^*$. Let $(u_\varepsilon)_{0<\varepsilon\le 1}$ be a family of solutions of

$$ (1.3) \qquad \Delta_g u_\varepsilon + h_\varepsilon u_\varepsilon = |u_\varepsilon|^{p_\varepsilon - 2} u_\varepsilon \quad \text{in } M. $$

We assume that $(u_\varepsilon)_{0<\varepsilon\le 1}$ satisfies

$$ (1.4) \qquad u_\varepsilon = cB_\varepsilon + o(1) \quad \text{in } H^1(M) $$

for some $c > 0$, where $B = (B_\varepsilon)_{\varepsilon>0}$ is a bubble centered at $(x_\varepsilon)_\varepsilon \in M$ and with radius $(\mu_\varepsilon)_\varepsilon \in (0, +\infty)$ and modeled on $V \in D^{1,2}(\mathbb{R}^n)$ as in Definition 3.1 below.

Then $c = 1$, $\lim_{\varepsilon\to 0} \mu_\varepsilon^{2^*-p_\varepsilon} = 1$, the following limit exists:

$$ \Lambda := \lim_{\varepsilon \to 0} \frac{2^* - p_\varepsilon}{\mu_\varepsilon^2} \ge 0, $$

and $h$ satisfies:

$$ (1.5) \qquad h(x_0) - \frac{n-2}{4(n-1)} S_g(x_0) = \frac{(n-2)^2}{4n} \frac{\int_{\mathbb{R}^n} |V|^{2^*} dx}{\int_{\mathbb{R}^n} V^2 dx} (\Lambda - \text{Weyl}_g \otimes B), $$

where we have let $x_0 := \lim_{\varepsilon \to 0} x_\varepsilon$.

The new term $\text{Weyl}_g \otimes B$ appearing in (1.5) is well-defined when $n \ge 5$ and we prove in Section 3 that it does not depend on the representative chosen for the family $(B_\varepsilon)_{\varepsilon>0}$ appearing in (1.4) (see Proposition 3.1 below).

In dimensions 3 and 4 we prove the following analogue of Theorem 1.1:
a non-zero solution $V$ to

$$ (1.2) \qquad \Delta_{\xi} V = |V|^{2^* - 2} V \quad \text{in } \mathbb{R}^n, \quad V \in D^{1,2}(\mathbb{R}^n), $$

where $D^{1,2}(\mathbb{R}^n)$ is the completion of $C_c^\infty(\mathbb{R}^n)$ for the norm $u \mapsto \|\nabla u\|_2$, and families of points $(x_\varepsilon)_{\varepsilon>0}$ in $M$ and of positive real numbers $(\mu_\varepsilon)_{\varepsilon>0}$ such that

$$ B_{\varepsilon}(x) - \chi(d_g(\cdot, x_{\varepsilon}))\mu_{\varepsilon}^{-\frac{n-2}{2}}V\left(\frac{\exp^{-1}_{x_{\varepsilon}}(\cdot)}{\mu_{\varepsilon}}\right) \to 0 \text{ in } H^1(M) $$

as $\varepsilon \to 0$. Here $\exp_{x_\varepsilon}$ is the exponential chart at $x_\varepsilon$ and $\xi$ is the Euclidean metric, so that $\Delta_\xi = -\sum_{i=1}^n \partial_i^2$. We will say that such a family $(B_\varepsilon)_{\varepsilon>0}$ is centered at $x_\varepsilon$ and of radius $\mu_\varepsilon$ and modeled on $V$, and $V$ itself will sometimes be referred to as the bubble. A more precise definition is given in Definition 3.1 and properties of such families are investigated in Section 3 below. All the solutions considered in this paper are allowed to change sign unless mentioned otherwise.

Our main results establish a necessary condition for the existence of blowing-up solutions of (1.1) that blow-up with a single sign-changing bubble. We introduce one last notation as follows: if $B = (B_\varepsilon)_{\varepsilon>0}$ is a family of bubbles centered at $x_\varepsilon$, of radius $\mu_\varepsilon$ and modeled on $V$, we define

$$ \begin{aligned} \text{Weyl}_g \otimes B := & \frac{4n}{3(n-2)^2} \left( \int_{\mathbb{R}^n} |V|^{2^*} dx \right)^{-1} \times \\ & \text{Weyl}_g(x_0)_{i\alpha j\beta} \int_{\mathbb{R}^n} x^\alpha x^\beta \partial_{ij}^2 V \left( \frac{n-2}{2} V + x^l \partial_l V \right) dx \end{aligned} $$

where $x_0 = \lim_{\varepsilon \to 0} x_\varepsilon$. We first state a result when $n \ge 5$.

**Theorem 1.1.** Let $(M^n, g), n \ge 5$ be a smooth, connected and closed Riemannian manifold, let $h, (h_\varepsilon)_{0<\varepsilon\le 1} \in C^1(M)$ such that $\lim_{\varepsilon\to 0} h_\varepsilon = h$ in $C^1(M)$ and $\Delta_g + h$ is coercive. Let $(p_\varepsilon)_{0<\varepsilon\le 1} \in (2, 2^*]$ be such that $\lim_{\varepsilon\to 0} p_\varepsilon = 2^*$. Let $(u_\varepsilon)_{0<\varepsilon\le 1}$ be a family of solutions of

$$ (1.3) \qquad \Delta_g u_\varepsilon + h_\varepsilon u_\varepsilon = |u_\varepsilon|^{p_\varepsilon-2} u_\varepsilon \quad \text{in } M. $$

We assume that $(u_\varepsilon)_{0<\varepsilon\le 1}$ satisfies

$$ (1.4) \qquad u_\varepsilon = cB_\varepsilon + o(1) \quad \text{in } H^1(M) $$

for some $c > 0$, where $B = (B_\varepsilon)_{\varepsilon>0}$ is a bubble centered at $(x_\varepsilon)_\varepsilon \in M$ and with radius $(\mu_\varepsilon)_\varepsilon \in (0, +\infty)$ and modeled on $V \in D^{1,2}(\mathbb{R}^n)$ as in Definition 3.1 below.

Then $c = 1$, $\lim_{\varepsilon\to 0} \mu_\varepsilon^{2^*-p_\varepsilon} = 1$, the following limit exists:

$$ \Lambda := \lim_{\varepsilon \to 0} \frac{2^* - p_\varepsilon}{\mu_\varepsilon^2} \ge 0, $$

and $h$ satisfies:

$$ (1.5) \qquad h(x_0) - \frac{n-2}{4(n-1)} S_g(x_0) = \frac{(n-2)^2}{4n} \frac{\int_{\mathbb{R}^n} |V|^{2^*} dx}{\int_{\mathbb{R}^n} V^2 dx} (\Lambda - \text{Weyl}_g \otimes B), $$

where we have let $x_0 := \lim_{\varepsilon \to 0} x_\varepsilon$.

The new term $\text{Weyl}_g \otimes B$ appearing in (1.5) is well-defined when $n \ge 5$ and we prove in Section 3 that it does not depend on the representative chosen for the family $(B_\varepsilon)_{\varepsilon>0}$ appearing in (1.4) (see Proposition 3.1 below).

In dimensions 3 and 4 we prove the following analogue of Theorem 1.1:
**Theorem 1.2.** Let $(M^n, g), n = 3, 4$ be a smooth, connected and closed manifold, let $h, (h_\epsilon)_{0<\epsilon\le 1}, (p_\epsilon)_{0<\epsilon\le 1} \in (2, 2^*)$ be as in the statement of Theorem 1.1. Let $(u_\epsilon)_{0<\epsilon\le 1}$ be a family of solutions of (1.3) satisfying again (1.4) for some $c > 0$, where $B = (B_\epsilon)_{\epsilon>0}$ is a bubble centered at $(x_\epsilon)_\epsilon \in M$ and with radius $(\mu_\epsilon)_\epsilon \in (0, +\infty)$ and modeled on $V \in D^{1,2}(\mathbb{R}^n)$ as in Definition 3.1 below.

Then $c = 1$, $\lim_{\epsilon \to 0} \mu_\epsilon^{2^*-p_\epsilon} = 1$, the following limits exist:

$$ \Lambda := \lim_{\epsilon \to 0} \left\{ \begin{array}{ll} \frac{6 - p_\epsilon}{\mu_\epsilon} & n = 3 \\ \frac{4 - p_\epsilon}{\mu_\epsilon^2 \ln \frac{1}{\mu_\epsilon}} & n = 4 \end{array} \right\} \ge 0, $$

and $h$ satisfies

$$ (1.6) \qquad \frac{\left(\int_{\mathbb{R}^4} |V|^2 V \,dx\right)^2}{\omega_3 \int_{\mathbb{R}^4} |V|^4 \,dx} \left(h(x_0) - \frac{1}{6} S_g(x_0)\right) = \Lambda \quad \text{when } n = 4 $$

$$ (1.7) \qquad -6 \frac{\left(\int_{\mathbb{R}^3} |V|^4 V \,dx\right)^2}{\int_{\mathbb{R}^3} |V|^6 \,dx} m_h(x_0) = \Lambda \quad \text{when } n = 3, $$

where as before we have let $x_0 := \lim_{\epsilon \to 0} x_\epsilon$. In (1.6) $\omega_3$ is the area of the 3-sphere and in (1.7) $m_h(x_0) \in \mathbb{R}$ is the mass of the Green's function of $\Delta_g + h$ at $x_0 \in M$ (see Definition 3.3 below).

A few comments on Theorems 1.1 and 1.2 are in order. Assumption (1.4) implies that $(u_\epsilon)_{\epsilon>0}$ blows-up as a single, possibly sign-changing bubble. In other words, (1.4) is a Struwe decomposition with a single bubbling profile. No assumption is made on the energy of $V$ (except that it is finite) and $V$ is not assumed to be non-degenerate in the sense of Duyckaerts-Kenig-Merle [13]. Assumption (1.4) should be regarded, insofar as sign-changing solutions are considered, as the analogue of the classical notion of isolated and simple blow-up point for positive solutions introduced in Schoen [35].

Theorems 1.1 and 1.2 are not new for positive solutions. It has been known since Obata [26] and Struwe [36] that when $(u_\epsilon)_{\epsilon>0}$ is a family of positive solutions of (1.3) the bubble $V$ is positive and radial and given by (2.2) below. In this case, and when $n \ge 5$, the additional term $Weyl_g \otimes B$ in (1.5) vanishes by the symmetries of $Weyl_g$ (see Proposition 3.3 below), and condition (1.5) becomes

$$ h(x_0) - \frac{n-2}{4(n-1)} S_g(x_0) = C(n)\Lambda \ge 0, \quad \text{with } \Lambda := \lim_{\epsilon \to 0} \frac{2^* - p_\epsilon}{\mu_\epsilon^2} $$

for some dimensional constant $C(n) > 0$. When $n = 3, 4$, (1.6) and (1.7) show that

$$ h(x_0) - \frac{1}{6} S_g(x_0) = C(4)\Lambda \quad \text{when } n = 4 \quad \text{and} \quad m_h(x_0) = -C(3)\Lambda \quad \text{when } n = 3, $$

where $\Lambda$ is as in the statement of the Theorem. When $p_\epsilon = 2^*$ for all $\epsilon$, and for positive finite-energy solutions, these conditions have been known since the work of Li-Zhu [20] and Druet [10, 11].

If $(u_\epsilon)_{\epsilon>0}$ is a blowing-up family of sign-changing solutions of (1.3) with a single bubbling profile as in (1.4), however, the bubble $V$ is sign-changing in general and (1.5), (1.6) and (1.7) display new phenomena. When $n \ge 5$ the additional term $Weyl_g \otimes B$ cannot be expected to vanish a priori. We construct indeed in Section 3
**Theorem 1.2.** Let $(M^n, g), n = 3, 4$ be a smooth, connected and closed Riemannian manifold, let $h, (h_\epsilon)_{0<\epsilon\le 1}, (p_\epsilon)_{0<\epsilon\le 1} \in (2, 2^*)$ be as in the statement of Theorem 1.1. Let $(u_\epsilon)_{0<\epsilon\le 1}$ be a family of solutions of (1.3) satisfying again (1.4) for some $c > 0$, where $B = (B_\epsilon)_{\epsilon>0}$ is a bubble centered at $(x_\epsilon)_\epsilon \in M$ and with radius $(\mu_\epsilon)_\epsilon \in (0, +\infty)$ and modeled on $V \in D^{1,2}(\mathbb{R}^n)$ as in Definition 3.1 below.

Then $c = 1$, $\lim_{\epsilon \to 0} \mu_{\epsilon}^{2^*-p_{\epsilon}} = 1$, the following limits exist:

$$ \Lambda := \lim_{\epsilon \to 0} \left\{ \begin{array}{ll} \frac{6 - p_{\epsilon}}{\mu_{\epsilon}} & n = 3 \\ \frac{4 - p_{\epsilon}}{\mu_{\epsilon}^2 \ln \frac{1}{\mu_{\epsilon}}} & n = 4 \end{array} \right\} \geq 0, $$

and $h$ satisfies

$$ (1.6) \qquad \frac{\left(\int_{\mathbb{R}^4} |V|^2 V \,dx\right)^2}{\omega_3 \int_{\mathbb{R}^4} |V|^4 \,dx} \left(h(x_0) - \frac{1}{6} S_g(x_0)\right) = \Lambda \quad \text{when } n = 4 $$

$$ (1.7) \qquad -6 \frac{\left(\int_{\mathbb{R}^3} |V|^4 V \,dx\right)^2}{\int_{\mathbb{R}^3} |V|^6 \,dx} m_h(x_0) = \Lambda \quad \text{when } n = 3, $$

where as before we have let $x_0 := \lim_{\epsilon \to 0} x_\epsilon$. In (1.6) $\omega_3$ is the area of the 3-sphere and in (1.7) $m_h(x_0) \in \mathbb{R}$ is the mass of the Green's function of $\Delta_g + h$ at $x_0 \in M$ (see Definition 3.3 below).

A few comments on Theorems 1.1 and 1.2 are in order. Assumption (1.4) implies that $(u_\epsilon)_{\epsilon>0}$ blows-up as a single, possibly sign-changing bubble. In other words, (1.4) is a Struwe decomposition with a single bubbling profile. No assumption is made on the energy of $V$ (except that it is finite) and $V$ is not assumed to be non-degenerate in the sense of Duyckaerts-Kenig-Merle [13]. Assumption (1.4) should be regarded, insofar as sign-changing solutions are considered, as the analogue of the classical notion of isolated and simple blow-up point for positive solutions introduced in Schoen [35].

Theorems 1.1 and 1.2 are not new for positive solutions. It has been known since Obata [26] and Struwe [36] that when $(u_\epsilon)_{\epsilon>0}$ is a family of positive solutions of (1.3) the bubble $V$ is positive and radial and given by (2.2) below. In this case, and when $n \ge 5$, the additional term $Weyl_g \otimes B$ in (1.5) vanishes by the symmetries of $Weyl_g$ (see Proposition 3.3 below), and condition (1.5) becomes

$$ h(x_0) - \frac{n-2}{4(n-1)} S_g(x_0) = C(n)\Lambda \ge 0, \quad \text{with } \Lambda := \lim_{\epsilon \to 0} \frac{2^* - p_{\epsilon}}{\mu_{\epsilon}^2} $$

for some dimensional constant $C(n) > 0$. When $n = 3, 4$, (1.6) and (1.7) show that

$$ h(x_0) - \frac{1}{6} S_g(x_0) = C(4)\Lambda \quad \text{when } n = 4 \quad \text{and} \quad m_h(x_0) = -C(3)\Lambda \quad \text{when } n = 3, $$

where $\Lambda$ is as in the statement of the Theorem. When $p_\epsilon = 2^*$ for all $\epsilon$, and for positive finite-energy solutions, these conditions have been known since the work of Li-Zhu [20] and Druet [10, 11].

If $(u_\epsilon)_{\epsilon>0}$ is a blowing-up family of sign-changing solutions of (1.3) with a single bubbling profile as in (1.4), however, the bubble $V$ is sign-changing in general and (1.5), (1.6) and (1.7) display new phenomena. When $n \ge 5$ the additional term $Weyl_g \otimes B$ cannot be expected to vanish a priori. We construct indeed in Section
examples of manifolds $(M, g)$ and of solutions $V$ of (1.2) such that $Weyl_g \otimes B \neq 0$ (see Proposition 3.4 below). The term $Weyl_g \otimes B$ is thus a striking new feature of sign-changing blow-up that does not appear for positive solutions. When $n \geq 5$, condition (1.5) is entirely new and highlights a strong interaction between $h$, the geometry of $(M, g)$ and the limiting bubble $V$ itself. When $n = 3, 4, (1.6)$ and $(1.7)$ involve $\int_{\mathbb{R}^n} |V|^{2^*-2}Vdx$. But if $V$ is a sign-changing solution of (1.2) for $n \geq 3$, the latter integral may vanish if $V$ decays strongly at infinity (see (2.9) below). When this is the case conditions (1.6) and (1.7) indicate that, in dimensions 3 and 4, the interaction between $B_\varepsilon$ and the geometry of $(M, g)$ takes place at a higher-order than it did for positive solutions.

**1.2. Consequences of Theorem 1.1 on the blow-up picture of (1.3).** Theorems 1.1 and 1.2 shed new light on the stability of the set of sign-changing solutions of (1.3). We say that equation (1.1), with $p = 2^*$, is stable if for every family $(h_\varepsilon)_{\varepsilon>0}$ converging to $h$, every energy-bounded family of solutions $(u_\varepsilon)_{\varepsilon>0}$ of (1.3) converges in $C^2(M)$, up to a subsequence. Druet [10] proved that, for positive solutions, stability holds provided $h \neq \frac{n-2}{4(n-1)}S_g$ everywhere in $M$ (with a caveat in dimensions 3 and 6, see [10]). This result was generalised to the case of sign-changing solutions, when $n \geq 7$ and when $(M, g)$ is locally conformally flat, in Premoselli-Vétois [31]. When $(M, g)$ is not locally conformally flat, and at least when $n \geq 5$, however, (1.5) shows that necessary conditions for a one bubble blow-up involve $h$, the geometry of $(M, g)$ and the possible limiting bubble $V$ itself. The set of sign-changing solutions of (1.2) is still poorly understood and no classification result is known, but examples of large-energy sequences of solutions have been constructed by Del Pino-Musso-Pacard-Pistoia [7], Ding [9], Medina-Musso [22], Medina-Musso-Wei [23]. The wealth of solutions of (1.2) is the main obstacle to understanding the term $Weyl_g \otimes B$ in (1.5). It seems therefore highly unlikely, even in the single-bubble case, that stability for sign-changing solutions of (1.3) can be enforced solely by a global assumption on $h$, as was the case for positive solutions.

Theorems 1.1 and 1.2, however, allow us to rule out limiting bubbles $V$ that may appear in the blow-up in some cases. We start with an example in dimensions 3 and 4. If $V$ is a solution of (1.2) it satisfies

$$V(x) = \frac{\lambda(V)}{|x|^{n-2}} + O(|x|^{1-n}) \quad \text{as } |x| \to +\infty,$$

where $(n-2)\omega_{n-1}\lambda(V) = \int_{\mathbb{R}^n} |V|^{2^*-2}Vdx$ (see (2.4) and (2.9) below). Bubbles decaying at infinity as $|x|^{2-n}$ satisfy $\lambda(V) > 0$: this is the case for the positive bubbles given by (2.2) and for the bubbles of [7]. An immediate corollary of Theorem 1.2 is as follows:

**Corollary 1.3.** Let $(M^n, g)$, $n = 3, 4$, be a smooth, closed, connected Riemannian manifold and $h$, $(h_\varepsilon)_{0<\varepsilon\leq 1} \in C^1(M)$ be such that $\lim_{\varepsilon\to 0} h_\varepsilon = h$ in $C^1(M)$ and $\Delta_g + h$ is coercive. Let $(p_\varepsilon)_{0<\varepsilon\leq 1} \in (2, 2^*)$ be such that $\lim_{\varepsilon\to 0} p_\varepsilon = 2^*$. Let $(u_\varepsilon)_{0<\varepsilon\leq 1}$ be a family of solutions of (1.3) satisfying (1.4), where $B_\varepsilon$ is a family of bubbles centered at $x_\varepsilon$, of radius $\mu_\varepsilon$ and modeled on $V$. Assume that

* the mass of $\Delta_g + h$ is positive at every point of $M$ if $n = 3$

* $h < \frac{1}{6}S_g$ everywhere in $M$ if $n = 4$.

Then $\int_{\mathbb{R}^n} |V|^{2^*-2}Vdx = 0$. In particular, $V(x) = O(|x|^{1-n})$ as $|x| \to +\infty$ and $V$ is neither positive nor one of the bubbles of [7].
3 examples of Riemannian manifolds $(M, g)$ and of solutions $V$ of (1.2) such that Weyl$_g \otimes B \neq 0$ (see Proposition 3.4 below). The term Weyl$_g \otimes B$ is thus a striking new feature of sign-changing blow-up that does not appear for positive solutions. When $n \ge 5$, condition (1.5) is entirely new and highlights a strong interaction between $h$, the geometry of $(M, g)$ and the limiting bubble $V$ itself. When $n = 3, 4$, (1.6) and (1.7) involve $\int_{\mathbb{R}^n} |V|^{2*-2}Vdx$. But if $V$ is a sign-changing solution of (1.2) for $n \ge 3$, the latter integral may vanish if $V$ decays strongly at infinity (see (2.9) below). When this is the case conditions (1.6) and (1.7) indicate that, in dimensions 3 and 4, the interaction between $B_\varepsilon$ and the geometry of $(M, g)$ takes place at a higher-order than it did for positive solutions.

**1.2. Consequences of Theorem 1.1 on the blow-up picture of (1.3).** Theorems 1.1 and 1.2 shed new light on the stability of the set of sign-changing solutions of (1.3). We say that equation (1.1), with $p = 2^*$, is stable if for every family $(h_\varepsilon)_{\varepsilon>0}$ converging to $h$, every energy-bounded family of solutions $(u_\varepsilon)_{\varepsilon>0}$ of (1.3) converges in $C^2(M)$, up to a subsequence. Druet [10] proved that, for positive solutions, stability holds provided $h \neq \frac{n-2}{4(n-1)}S_g$ everywhere in $M$ (with a caveat in dimensions 3 and 6, see [10]). This result was generalised to the case of sign-changing solutions, when $n \ge 7$ and when $(M, g)$ is locally conformally flat, in Premoselli-Vétois [31]. When $(M, g)$ is not locally conformally flat, and at least when $n \ge 5$, however, (1.5) shows that necessary conditions for a one bubble blow-up involve $h$, the geometry of $(M, g)$ and the possible limiting bubble $V$ itself. The set of sign-changing solutions of (1.2) is still poorly understood and no classification result is known, but examples of large-energy sequences of solutions have been constructed by Del Pino-Musso-Pacard-Pistoia [7], Ding [9], Medina-Musso [22], Medina-Musso-Wei [23]. The wealth of solutions of (1.2) is the main obstacle to understanding the term Weyl$_g \otimes B$ in (1.5). It seems therefore highly unlikely, even in the single-bubble case, that stability for sign-changing solutions of (1.3) can be enforced solely by a global assumption on $h$, as was the case for positive solutions.

Theorems 1.1 and 1.2, however, allow us to rule out limiting bubbles $V$ that may appear in the blow-up in some cases. We start with an example in dimensions 3 and 4. If $V$ is a solution of (1.2) it satisfies

$$V(x) = \frac{\lambda(V)}{|x|^{n-2}} + O(|x|^{1-n}) \quad \text{as } |x| \to +\infty,$$

where $(n-2)\omega_{n-1}\lambda(V) = \int_{\mathbb{R}^n} |V|^{2*-2}Vdx$ (see (2.4) and (2.9) below). Bubbles decaying at infinity as $|x|^{2-n}$ satisfy $\lambda(V) > 0$: this is the case for the positive bubbles given by (2.2) and for the bubbles of [7]. An immediate corollary of Theorem 1.2 is as follows:

**Corollary 1.3.** Let $(M^n, g)$, $n = 3, 4$, be a smooth, closed, connected Riemannian manifold and $h$, $(h_\varepsilon)_{0<\varepsilon\le 1} \in C^1(M)$ be such that $\lim_{\varepsilon\to 0} h_\varepsilon = h$ in $C^1(M)$ and $\Delta_g + h$ is coercive. Let $(p_\varepsilon)_{0<\varepsilon\le 1} \in (2, 2^*)$ be such that $\lim_{\varepsilon\to 0} p_\varepsilon = 2^*$. Let $(u_\varepsilon)_{0<\varepsilon\le 1}$ be a family of solutions of (1.3) satisfying (1.4), where $B_\varepsilon$ is a family of bubbles centered at $x_\varepsilon$, of radius $\mu_\varepsilon$ and modeled on $V$. Assume that

* the mass of $\Delta_g + h$ is positive at every point of $M$ if $n = 3$

* $h < \frac{1}{6}S_g$ everywhere in $M$ if $n = 4$.

Then $\int_{\mathbb{R}^n} |V|^{2*-2}Vdx = 0$. In particular, $V(x) = O(|x|^{1-n})$ as $|x| \to +\infty$ and $V$ is neither positive nor one of the bubbles of [7].
In dimensions $n \ge 5$ (1.5) does not provide any information on the decay of $V$ at infinity, but we can still rule out the existence of symmetric bubbles in some cases. We recall that Ding [9] constructed a family of solutions of (1.2), possessing $O(p) \times O(n-p)$ symmetry for $2 \le p \le n-1$. We show that Ding solutions cannot appear as bubbles in a single-bubble blow-up on symmetric manifolds:

**Corollary 1.4.** Let $(M^n, g) = (\mathbb{S}^p \times \mathbb{S}^q, g_p \times g_q)$, where $g_p$ and $g_q$ are the round metrics on $\mathbb{S}^p$ and $\mathbb{S}^q$, $p, q \ge 2$ and $n = p+q \ge 5$. Let $h, (h_\epsilon)_{0<\epsilon\le 1} \in C^1(M)$ be such that $\lim_{\epsilon\to 0} h_\epsilon = h$ in $C^1(M)$ and $\Delta_g + h$ is coercive. Let $(p_\epsilon)_{0<\epsilon\le 1} \in (2, 2^*)$ be such that $\lim_{\epsilon\to 0} p_\epsilon = 2^*$. Let $(u_\epsilon)_{0<\epsilon\le 1}$ be a family of solutions of (1.3) satisfying (1.4), where $B_\epsilon$ is a family of bubbles centered at $x_\epsilon$, of radius $\mu_\epsilon$, modeled on $V$. Assume that $h \le \frac{n-2}{4(n-1)}S_g$ everywhere. Then

$$Weyl_g \otimes B \ge 0.$$

As a consequence, $V$ cannot be a Ding solution.

*Proof.* Since $h \le \frac{n-2}{4(n-1)}S_g$ everywhere, (1.5) shows that $Weyl_g \otimes B \ge 0$. But as shown in Corollary 3.5 below, $Weyl_g \otimes B < 0$ if $V$ is a Ding solution. $\square$

Corollaries 1.3 (when $n=3$) and 1.4 both apply to one-bubble sign-changing blow-up for the sign-changing Yamabe equation:

$$\Delta_g u + \frac{n-2}{4(n-1)} S_g u = |u|^{2^*-2} u.$$

Solutions of this equation naturally appear as minimisers for the second conformal eigenvalue (see Ammann-Humbert [1]), and their stability properties shed new light on the set of possible minimisers (see Premoselli-Vétois [30]). Other stability results for sign-changing solutions are in Premoselli-Vétois [29, 31]. Examples of sign-changing families of solutions of (1.3) have been constructed in Bonheure-Casteras-Premoselli [5], Micheletti-Pistoia-Vétois [25], Pistoia-Vétois [27] and Robert-Vétois [33, 34]. Deng-Musso-Wei [8] have considered an innovative construction of a family of sign-changing solutions $(u_\epsilon)_{\epsilon>0}$ to the (1.3) with $h_\epsilon \equiv h$ and $p_\epsilon := 2^* - \epsilon$ that blows-up like a single bubble as in (1.4). We discuss similar constructions in Section 6.

The structure of the paper is as follows. In Section 2 we introduce a few properties of sign-changing solutions of (1.2) and highlight the main differences with respect to the positive case. Section 3 contains the formal definition and properties of the families of bubbles that we investigate here. We also investigate the term $Weyl_g \otimes B$ there and construct examples of manifolds and bubbles where $Weyl_g \otimes B \ne 0$. Section 4 is the core of the analysis of this paper and we prove there that relation (1.4) can be improved into optimal pointwise bounds on $u_\epsilon$. This is the content of Proposition 4.1 and is the main ingredient in the proof of Theorems 1.1 and 1.2. Finally, Theorems 1.1 and 1.2 are proven in Section 5.

**Acknowledgements:** the authors would like to thank P.-D. Thizy and J. Vétois for many fruitful discussions in the early version of this work.

## 2. FINITE-ENERGY NODAL SOLUTIONS OF THE YAMABE EQUATION IN $\mathbb{R}^n$

We let $\Sigma$ be the set of non-zero finite-energy solutions of (1.2):

$$(2.1) \qquad \Sigma = \{V \in D^{1,2}(\mathbb{R}^n) \setminus \{0\}, \quad \Delta_{\xi} V = |V|^{2^*-2} V\}.$$
In dimensions $n \ge 5$ (1.5) does not provide any information on the decay of $V$ at infinity, but we can still rule out the existence of symmetric bubbles in some cases. We recall that Ding [9] constructed a family of solutions of (1.2), possessing $O(p) \times O(n-p)$ symmetry for $2 \le p \le n-1$. We show that Ding solutions cannot appear as bubbles in a single-bubble blow-up on symmetric manifolds:

**Corollary 1.4.** Let $(M^n, g) = (\mathbb{S}^p \times \mathbb{S}^q, g_p \times g_q)$, where $g_p$ and $g_q$ are the round metrics on $\mathbb{S}^p$ and $\mathbb{S}^q$, $p, q \ge 2$ and $n = p+q \ge 5$. Let $h, (h_\varepsilon)_{0<\varepsilon\le 1} \in C^1(M)$ be such that $\lim_{\varepsilon\to 0} h_\varepsilon = h$ in $C^1(M)$ and $\Delta_g + h$ is coercive. Let $(p_\varepsilon)_{0<\varepsilon\le 1} \in (2, 2^*)$ be such that $\lim_{\varepsilon\to 0} p_\varepsilon = 2^*$. Let $(u_\varepsilon)_{0<\varepsilon\le 1}$ be a family of solutions of (1.3) satisfying (1.4), where $B_\varepsilon$ is a family of bubbles centered at $x_\varepsilon$, of radius $\mu_\varepsilon$, modeled on $V$. Assume that $h \le \frac{n-2}{4(n-1)}S_g$ everywhere. Then

$$Weyl_g \otimes B \ge 0.$$

As a consequence, $V$ cannot be a Ding solution.

*Proof.* Since $h \le \frac{n-2}{4(n-1)}S_g$ everywhere, (1.5) shows that $Weyl_g \otimes B \ge 0$. But as shown in Corollary 3.5 below, $Weyl_g \otimes B < 0$ if $V$ is a Ding solution. $\square$

Corollaries 1.3 (when $n=3$) and 1.4 both apply to one-bubble sign-changing blow-up for the sign-changing Yamabe equation:

$$\Delta_g u + \frac{n-2}{4(n-1)} S_g u = |u|^{2^*-2} u.$$

Solutions of this equation naturally appear as minimisers for the second conformal eigenvalue (see Ammann-Humbert [1]), and their stability properties shed new light on the set of possible minimisers (see Premoselli-Vétois [30]). Other stability results for sign-changing solutions are in Premoselli-Vétois [29, 31]. Examples of sign-changing families of solutions of (1.3) have been constructed in Bonheure-Casteras-Premoselli [5], Micheletti-Pistoia-Vétois [25], Pistoia-Vétois [27] and Robert-Vétois [33, 34]. Deng-Musso-Wei [8] have considered an innovative construction of a family of sign-changing solutions $(u_\varepsilon)_{\varepsilon>0}$ to the (1.3) with $h_\varepsilon \equiv h$ and $p_\varepsilon := 2^* - \varepsilon$ that blows-up like a single bubble as in (1.4). We discuss similar constructions in Section 6.

The structure of the paper is as follows. In Section 2 we introduce a few properties of sign-changing solutions of (1.2) and highlight the main differences with respect to the positive case. Section 3 contains the formal definition and properties of the families of bubbles that we investigate here. We also investigate the term $Weyl_g \otimes B$ there and construct examples of manifolds and bubbles where $Weyl_g \otimes B \ne 0$. Section 4 is the core of the analysis of this paper and we prove there that relation (1.4) can be improved into optimal pointwise bounds on $u_\varepsilon$. This is the content of Proposition 4.1 and is the main ingredient in the proof of Theorems 1.1 and 1.2. Finally, Theorems 1.1 and 1.2 are proven in Section 5.

**Acknowledgements:** the authors would like to thank P.-D. Thizy and J. Vétois for many fruitful discussions in the early version of this work.

## 2. FINITE-ENERGY NODAL SOLUTIONS OF THE YAMABE EQUATION IN $\mathbb{R}^n$

We let $\Sigma$ be the set of non-zero finite-energy solutions of (1.2):

$$(2.1) \qquad \Sigma = \{V \in D^{1,2}(\mathbb{R}^n) \setminus \{0\}, \quad \Delta_{\xi} V = |V|^{2^*-2} V\}.$$
Regularity theory for critical equations (see eg Trudinger [37]) together with standard elliptic theory, shows that there exists $\alpha_n > 0$ such that $V \in C_{\text{loc}}^{3, \alpha_n}(\mathbb{R}^n)$ for all $V \in \Sigma$. If $V \in \Sigma$ is positive it is equal by Obata [26], up to translations and rescalings, to

$$ (2.2) \qquad B_0^+(x) = \left(1 + \frac{|x|^2}{n(n-2)}\right)^{-\frac{n-2}{2}}, \quad x \in \mathbb{R}^n. $$

By Caffarelli-Gidas-Spruck [6] this result remains true for all positive solutions of (1.2), without the assumption that they belong to $D^{1,2}(\mathbb{R}^n)$. Note also that, up to dilation, $B_0^+$ is the only non-zero radial solution of (1.2) which belongs to $D^{1,2}(\mathbb{R}^n)$ by a simple application of Pohozaev's identity.

If $V \in D^{1,2}(\mathbb{R}^n)$ we define its Kelvin transform as

$$ (2.3) \qquad V^*(x) = \frac{1}{|x|^{n-2}} V\left(\frac{x}{|x|^2}\right) \quad \text{for a.e. } x \in \mathbb{R}^n \setminus \{0\}. $$

It is well-known (see for instance Duyckaerts-Kenig-Merle [13, Proposition 3.1]) that $V \mapsto V^*$ defines an isometry of $D^{1,2}(\mathbb{R}^n)$ and of $L^{2^*}(\mathbb{R}^n)$ and that

$$ \Delta_{\xi} V^*(x) = \frac{1}{|x|^{n+2}} \Delta_{\xi} V \left( \frac{x}{|x|^2} \right) $$

in $\mathbb{R}^n \setminus \{0\}$, provided $V$ is of class $C^2$. As a consequence we have $V^* \in \Sigma$ whenever $V \in \Sigma$ (see again [13]) and the regularity theory for (1.2) shows that $V^* \in C^3(\mathbb{R}^n)$ when $V \in \Sigma$. The following result provides a precise description of the behavior at infinity of the elements in $\Sigma$:

**Lemma 2.1.** Let $V \in \Sigma$. There exist $\lambda(V) \in \mathbb{R}$ and $\alpha(V) \in \mathbb{R}^n$ such that the following asymptotic expansion holds as $|x| \to +\infty$ and can be differentiated:

$$ (2.4) \qquad V(x) = \frac{\lambda(V)}{|x|^{n-2}} + \frac{\langle \alpha(V), x \rangle}{|x|^n} + O\left(\frac{1}{|x|^n}\right). $$

When we say that expansion (2.4) can be differentiated we mean that the following holds as $|x| \to +\infty$: for any $1 \le i \le n$,

$$ (2.5) \qquad \begin{aligned} \partial_i V(x) = -(n-2)\lambda(V) \frac{x_i}{|x|^n} &+ \frac{\alpha(V)_i |x|^2 - n \langle \alpha(V), x \rangle x_i}{|x|^{n+2}} \\ &+ O\left(\frac{1}{|x|^{n+1}}\right). \end{aligned} $$

*Proof.* Let $V^*$ be, as before, the Kelvin transform of $V$. It is of class $C^3$ in $\mathbb{R}^n$ and we can consider its Taylor expansion at $0$: there exist $\lambda(V) \in \mathbb{R}$ and $\alpha(V) \in \mathbb{R}^n$ such that, as $y \to 0$, and for $1 \le i \le n$,

$$ (2.6) \qquad \begin{aligned} V^*(y) &= \lambda(V) + \langle \alpha(V), y \rangle + O(|y|^2) \\ \partial_i V^*(y) &= \alpha(V)_i + O(|y|) \end{aligned} $$

hold. Since for any $x \ne 0$ we have $V(x) = \frac{1}{|x|^{n-2}} V^*\left(\frac{x}{|x|^2}\right)$, (2.4) follows from (2.6). To prove (2.5) we write that by the chain rule, for $x \ne 0$,

$$ (2.7) \qquad \partial_i V(x) = -(n-2) \frac{x_i}{|x|^n} V^* \left( \frac{x}{|x|^2} \right) + \sum_{j=1}^{n} \frac{1}{|x|^{n-2}} \left( \frac{\delta_{ij}}{|x|^2} - 2 \frac{x_i x_j}{|x|^4} \right) \partial_j V^* \left( \frac{x}{|x|^2} \right) $$
Regularity theory for critical equations (see eg Trudinger [37]) together with standard elliptic theory, shows that there exists $\alpha_n > 0$ such that $V \in C_{\text{loc}}^{3,\alpha_n}(\mathbb{R}^n)$ for all $V \in \Sigma$. If $V \in \Sigma$ is positive it is equal by Obata [26], up to translations and rescalings, to

$$ (2.2) \qquad B_0^+(x) = \left(1 + \frac{|x|^2}{n(n-2)}\right)^{-\frac{n-2}{2}}, \quad x \in \mathbb{R}^n. $$

By Caffarelli-Gidas-Spruck [6] this result remains true for all positive solutions of (1.2), without the assumption that they belong to $D^{1,2}(\mathbb{R}^n)$. Note also that, up to dilation, $B_0^+$ is the only non-zero radial solution of (1.2) which belongs to $D^{1,2}(\mathbb{R}^n)$ by a simple application of Pohozaev's identity.

If $V \in D^{1,2}(\mathbb{R}^n)$ we define its Kelvin transform as

$$ (2.3) \qquad V^*(x) = \frac{1}{|x|^{n-2}} V\left(\frac{x}{|x|^2}\right) \quad \text{for a.e. } x \in \mathbb{R}^n \setminus \{0\}. $$

It is well-known (see for instance Duyckaerts-Kenig-Merle [13, Proposition 3.1]) that $V \mapsto V^*$ defines an isometry of $D^{1,2}(\mathbb{R}^n)$ and of $L^{2^*}(\mathbb{R}^n)$ and that

$$ \Delta_{\xi} V^*(x) = \frac{1}{|x|^{n+2}} \Delta_{\xi} V \left( \frac{x}{|x|^2} \right) $$

in $\mathbb{R}^n \setminus \{0\}$, provided $V$ is of class $C^2$. As a consequence we have $V^* \in \Sigma$ whenever $V \in \Sigma$ (see again [13]) and the regularity theory for (1.2) shows that $V^* \in C^3(\mathbb{R}^n)$ when $V \in \Sigma$. The following result provides a precise description of the behavior at infinity of the elements in $\Sigma$:

**Lemma 2.1.** Let $V \in \Sigma$. There exist $\lambda(V) \in \mathbb{R}$ and $\alpha(V) \in \mathbb{R}^n$ such that the following asymptotic expansion holds as $|x| \to +\infty$ and can be differentiated:

$$ (2.4) \qquad V(x) = \frac{\lambda(V)}{|x|^{n-2}} + \frac{\langle \alpha(V), x \rangle}{|x|^n} + O\left(\frac{1}{|x|^n}\right). $$

When we say that expansion (2.4) can be differentiated we mean that the following holds as $|x| \to +\infty$: for any $1 \le i \le n$,

$$ (2.5) \qquad \begin{aligned} \partial_i V(x) = &-(n-2)\lambda(V) \frac{x_i}{|x|^n} + \frac{\alpha(V)_i |x|^2 - n \langle \alpha(V), x \rangle x_i}{|x|^{n+2}} \\ &+ O\left(\frac{1}{|x|^{n+1}}\right). \end{aligned} $$

*Proof.* Let $V^*$ be, as before, the Kelvin transform of $V$. It is of class $C^3$ in $\mathbb{R}^n$ and we can consider its Taylor expansion at $0$: there exist $\lambda(V) \in \mathbb{R}$ and $\alpha(V) \in \mathbb{R}^n$ such that, as $y \to 0$, and for $1 \le i \le n$,

$$ (2.6) \qquad \begin{aligned} V^*(y) &= \lambda(V) + \langle \alpha(V), y \rangle + O(|y|^2) \\ \partial_i V^*(y) &= \alpha(V)_i + O(|y|) \end{aligned} $$

hold. Since for any $x \ne 0$ we have $V(x) = \frac{1}{|x|^{n-2}} V^*\left(\frac{x}{|x|^2}\right)$, (2.4) follows from (2.6). To prove (2.5) we write that by the chain rule, for $x \ne 0$,

$$ (2.7) \qquad \partial_i V(x) = -(n-2) \frac{x_i}{|x|^n} V^* \left( \frac{x}{|x|^2} \right) + \sum_{j=1}^{n} \frac{1}{|x|^{n-2}} \left( \frac{\delta_{ij}}{|x|^2} - 2 \frac{x_i x_j}{|x|^4} \right) \partial_j V^* \left( \frac{x}{|x|^2} \right) $$
holds, and we conclude using again (2.6).

<figure><img src="image_3.png" /></figure>

As a consequence of (2.4) and (2.5), there exists $C = C(n, V)$ such that

$$ (2.8) \quad |V(x)| + (1 + |x|)|\nabla V(x)| + (1 + |x|)^2 |\nabla^2 V(x)| \le C(1 + |x|)^{2-n} $$

for all $x \in \mathbb{R}^n$. This can also be shown by the conformal invariance of (1.2), see e.g. Premoselli [28, Lemma 2.2]. As shown in Premoselli [28, Lemma 2.1], we have

$$ (2.9) \quad \lambda(V) = \frac{1}{(n-2)\omega_{n-1}} \int_{\mathbb{R}^n} |V|^{2*-2}Vdx \quad \text{for all } V \in \Sigma. $$

If $B_0^+$ is as in (2.2) we have $\lambda(B_0^+) = (n(n-2))^{\frac{n-2}{2}}$. The solutions constructed in [7] also satisfy $\lambda(V) \ne 0$. As proven in Lemma 2.1, $\lambda(V) = V^*(0)$. As a consequence, examples where $\lambda(V) = 0$ are easily obtained by Kelvin-transforming sign-changing solutions of $\Sigma$ at a point where they vanish.

*Remark 2.2.* An accurate inspection of the proof of Lemma 2.1 shows that we can state a slightly more precise result. Let $V \in \Sigma$.

• Assume that $n = 3, 4$. Then there exists $k \in \mathbb{N}$ and a non-zero homogeneous polynomial $P$ of degree $k$ such that

$$ V(x) = \frac{P(x)}{|x|^{n-2+2k}} + O\left(\frac{1}{|x|^{n-1+k}}\right) \quad \text{as } |x| \to +\infty. $$

• Assume that $n = 5$. Then there exists $k \in \{0, \dots, 4\}$ and a homogeneous polynomial $P$ of degree $k$ such that

$$ V(x) = \frac{P(x)}{|x|^{3+2k}} + \begin{cases} O\left(\frac{1}{|x|^{4+k}}\right) & \text{if } k \le 3 \\ O\left(\frac{1}{|x|^{7+\frac{1}{3}}}\right) & \text{if } k = 4 \end{cases} \quad \text{as } |x| \to +\infty. $$

• Assume finally that $n \ge 6$. Let $0 < \alpha < 1$ if $n = 6$ and $0 < \alpha \le \frac{4}{n-2}$ if $n \ge 7$. Then there exists $k \in \{0, \dots, 3\}$ and a homogeneous polynomial $P$ of degree $k$ such that $V$ satisfies

$$ V(x) = \frac{P(x)}{|x|^{n-2+2k}} + O\left(\frac{1}{|x|^{n-2+k+\alpha}}\right) \quad \text{as } |x| \to +\infty. $$

All these expansions can be differentiated. In dimensions $n = 3, 4$, $P$ is the first non-zero homogeneous polynomial in the Taylor expansion of $V^*$ at 0, which always exists by classical finite continuation results (see e.g. Aronszajn [2]) since $V \in D^{1,2}(\mathbb{R}^n)\setminus\{0\}$. In dimensions $n \ge 5$, $P$ may vanish. Examples of solutions in any dimension $n \ge 3$ for which the degree of $P$ is larger or equal than 2 are still unknown and the question of their existence was raised in Duyckaerts-Kenig-Merle[13].

### 3. RIEMANNIAN PRELIMINARIES

3.1. **Sign-changing families of bubbles in $M$.** We define the notion of Riemannian bubble that we will consider in this work:

**Definition 3.1.** Let $(x_\varepsilon)_{\varepsilon>0} \in M$ and $(\mu_\varepsilon)_{\varepsilon>0} \in \mathbb{R}_{>0}$ be families of points in $M$ and positive real numbers, with $\mu_\varepsilon \to 0$ as $\varepsilon \to 0$. We let $x_0 := \lim_{\varepsilon\to 0} x_\varepsilon$. Let $V \in \Sigma$ be fixed and let $B = (B_\varepsilon)_{\varepsilon>0}$ be a family of functions in $H^1(M)$. We say that $B$ is a family of bubbles centered at $x_\varepsilon$ with radius $\mu_\varepsilon$ and modeled on $V$ if
holds, and we conclude using again (2.6).

<figure><img src="image_2.png" /></figure>

As a consequence of (2.4) and (2.5), there exists $C = C(n, V)$ such that

$$ (2.8) \qquad |V(x)| + (1 + |x|)|\nabla V(x)| + (1 + |x|)^2 |\nabla^2 V(x)| \le C(1 + |x|)^{2-n} $$

for all $x \in \mathbb{R}^n$. This can also be shown by the conformal invariance of (1.2), see e.g. Premoselli [28, Lemma 2.2]. As shown in Premoselli [28, Lemma 2.1], we have

$$ (2.9) \qquad \lambda(V) = \frac{1}{(n-2)\omega_{n-1}} \int_{\mathbb{R}^n} |V|^{2*-2}Vdx \quad \text{for all } V \in \Sigma. $$

If $B_0^+$ is as in (2.2) we have $\lambda(B_0^+) = (n(n-2))^{\frac{n-2}{2}}$. The solutions constructed in [7] also satisfy $\lambda(V) \ne 0$. As proven in Lemma 2.1, $\lambda(V) = V^*(0)$. As a consequence, examples where $\lambda(V) = 0$ are easily obtained by Kelvin-transforming sign-changing solutions of $\Sigma$ at a point where they vanish.

*Remark 2.2.* An accurate inspection of the proof of Lemma 2.1 shows that we can state a slightly more precise result. Let $V \in \Sigma$.

• Assume that $n = 3, 4$. Then there exists $k \in \mathbb{N}$ and a non-zero homogeneous polynomial $P$ of degree $k$ such that

$$ V(x) = \frac{P(x)}{|x|^{n-2+2k}} + O\left(\frac{1}{|x|^{n-1+k}}\right) \quad \text{as } |x| \to +\infty. $$

• Assume that $n = 5$. Then there exists $k \in \{0, \dots, 4\}$ and a homogeneous polynomial $P$ of degree $k$ such that

$$ V(x) = \frac{P(x)}{|x|^{3+2k}} + \begin{cases} O\left(\frac{1}{|x|^{4+k}}\right) & \text{if } k \le 3 \\ O\left(\frac{1}{|x|^{7+\frac{1}{3}}}\right) & \text{if } k = 4 \end{cases} \quad \text{as } |x| \to +\infty. $$

• Assume finally that $n \ge 6$. Let $0 < \alpha < 1$ if $n = 6$ and $0 < \alpha \le \frac{4}{n-2}$ if $n \ge 7$. Then there exists $k \in \{0, \dots, 3\}$ and a homogeneous polynomial $P$ of degree $k$ such that $V$ satisfies

$$ V(x) = \frac{P(x)}{|x|^{n-2+2k}} + O\left(\frac{1}{|x|^{n-2+k+\alpha}}\right) \quad \text{as } |x| \to +\infty. $$

All these expansions can be differentiated. In dimensions $n = 3, 4$, $P$ is the first non-zero homogeneous polynomial in the Taylor expansion of $V^*$ at 0, which always exists by classical finite continuation results (see e.g. Aronszajn [2]) since $V \in D^{1,2}(\mathbb{R}^n)\setminus\{0\}$. In dimensions $n \ge 5$, $P$ may vanish. Examples of solutions in any dimension $n \ge 3$ for which the degree of $P$ is larger or equal than 2 are still unknown and the question of their existence was raised in Duyckaerts-Kenig-Merle[13].

### 3. RIEMANNIAN PRELIMINARIES

3.1. **Sign-changing families of bubbles in $M$.** We define the notion of Riemannian bubble that we will consider in this work:

**Definition 3.1.** Let $(x_\varepsilon)_{\varepsilon>0} \in M$ and $(\mu_\varepsilon)_{\varepsilon>0} \in \mathbb{R}_{>0}$ be families of points in $M$ and positive real numbers, with $\mu_\varepsilon \to 0$ as $\varepsilon \to 0$. We let $x_0 := \lim_{\varepsilon\to 0} x_\varepsilon$. Let $V \in \Sigma$ be fixed and let $B = (B_\varepsilon)_{\varepsilon>0}$ be a family of functions in $H^1(M)$. We say that $B$ is a family of bubbles centered at $x_\varepsilon$ with radius $\mu_\varepsilon$ and modeled on $V$ if
there exists a family of charts $(\varphi_{\varepsilon})_{\varepsilon \ge 0}$ such that $\varphi_{\varepsilon} : B_g(x_0, r_g) \to \mathbb{R}^n$ is normal at $x_{\varepsilon}$ for all $\varepsilon > 0$, $\lim_{\varepsilon \to 0} \varphi_{\varepsilon} = \varphi_0$ in $C^k$ for all $k$ and such that

$$ (3.1) \qquad B_{\varepsilon} - \chi(d_g(\cdot, x_{\varepsilon}))\mu_{\varepsilon}^{-\frac{n-2}{2}} V \left( \frac{\varphi_{\varepsilon}(\cdot) - \varphi_{\varepsilon}(x_{\varepsilon})}{\mu_{\varepsilon}} \right) \to 0 \quad \text{in } H^1(M) $$
as $\varepsilon \to 0$.

Here $\Sigma$ is as in (2.1). As the following simple result shows, Definition 3.1 is independent of the choice of the local chart up to an isometry. In the sequel, we let $O(n)$ be the group of isometries of the Euclidean space $\mathbb{R}^n$:

**Proposition 3.1.** Let $(x_{\varepsilon})_{\varepsilon>0} \in M$ and $(\mu_{\varepsilon})_{\varepsilon>0} \in \mathbb{R}_{>0}$ be families of points in $M$ and positive real numbers, with $\mu_{\varepsilon} \to 0$ as $\varepsilon \to 0$, and let $B = (B_{\varepsilon})_{\varepsilon}$ and $\bar{B} = (\bar{B}_{\varepsilon})_{\varepsilon}$ be families of bubbles centered at $x_{\varepsilon}$ with radius $\mu_{\varepsilon}$ and respectively modeled on $V, \bar{V} \in \Sigma$ as in Definition (3.1). Then

$$ B_{\varepsilon} - \bar{B}_{\varepsilon} \to 0 \text{ in } H^1(M) \text{ as } \varepsilon \to 0 \iff \bar{V} = V \circ O $$
for some orthogonal mapping $O \in O(n)$.

Families of bubbles $(B_{\varepsilon})_{\varepsilon>0}$ and $(\bar{B}_{\varepsilon})_{\varepsilon>0}$ centered at $x_{\varepsilon}$ with radius $\mu_{\varepsilon}$ and satisfying $B_{\varepsilon} - \bar{B}_{\varepsilon} \to 0$ in $H^1(M)$ as $\varepsilon \to 0$ will be called equivalent families.

*Proof.* Let $\varphi_{\varepsilon}, \bar{\varphi}_{\varepsilon} : B_g(x_0, r_g) \to \mathbb{R}^n$ be normal coordinate charts for $g$ at $x_{\varepsilon}$ associated to $B$ and $\bar{B}$ as in definition 3.1. If $x \in \mathbb{R}^n$ we have

$$ (3.2) \qquad \varphi_{\varepsilon} \circ \bar{\varphi}_{\varepsilon}^{-1}(\bar{\varphi}_{\varepsilon}(x_{\varepsilon}) + x) = \varphi_{\varepsilon}(x_{\varepsilon}) + O_{\varepsilon}x + O(|x|^2) \quad \text{as } x \to 0, $$

where we have let $O_{\varepsilon} = D\varphi_{\varepsilon}(x_{\varepsilon}) \circ (D\bar{\varphi}_{\varepsilon}(x_{\varepsilon}))^{-1}$. Since $\varphi_{\varepsilon}$ and $\bar{\varphi}_{\varepsilon}$ are normal charts at $x_{\varepsilon}$, $O_{\varepsilon} \in O(n)$. Straightforward computations using (2.8) and (3.2) show that

$$ \chi(d_g(\cdot, x_{\varepsilon})) \left[ \mu_{\varepsilon}^{-\frac{n-2}{2}} V \left( \frac{\varphi_{\varepsilon}(\cdot) - \varphi_{\varepsilon}(x_{\varepsilon})}{\mu_{\varepsilon}} \right) - \mu_{\varepsilon}^{-\frac{n-2}{2}} V \left( \frac{O_{\varepsilon}(\bar{\varphi}_{\varepsilon}(\cdot) - \bar{\varphi}_{\varepsilon}(x_{\varepsilon}))}{\mu_{\varepsilon}} \right) \right] \to 0 $$
in $H^1(M)$ as $\varepsilon \to 0$. Since $B_{\varepsilon} - \bar{B}_{\varepsilon} \to 0$ in $H^1(M)$ by assumption we obtain

$$ \chi(d_g(\cdot, x_{\varepsilon})) \left[ \mu_{\varepsilon}^{-\frac{n-2}{2}} V \left( \frac{O_{\varepsilon}(\bar{\varphi}_{\varepsilon}(\cdot) - \bar{\varphi}_{\varepsilon}(x_{\varepsilon}))}{\mu_{\varepsilon}} \right) - \mu_{\varepsilon}^{-\frac{n-2}{2}} \bar{V} \left( \frac{\bar{\varphi}_{\varepsilon}(\cdot) - \bar{\varphi}_{\varepsilon}(x_{\varepsilon})}{\mu_{\varepsilon}} \right) - \right] \to 0 $$
in $H^1(M)$ as $\varepsilon \to 0$. A simple scaling argument then shows that $V \circ O = \bar{V}$ where $O := \lim_{\varepsilon \to 0} O_{\varepsilon}$. $\square$

Equations (3.2) and (3.3) show that, up to composing $V$ with an euclidean isometry, we can always assume that $\varphi_{\varepsilon}$ is an exponential chart at $x_{\varepsilon}$ in Definition 3.1. Choose, for any $p \in B_g(x_0, r_g)$, an identification of $T_pM$ with $\mathbb{R}^n$. This defines an exponential map at $p$, that we will denote by $\exp_p : \mathbb{R}^n \to M$ in all of this paper. Mimicking the arguments in (3.2) and (3.3) we similarly have, for $V \in \Sigma$,

$$ \chi(d_g(\cdot, x_{\varepsilon})) \left[ \mu_{\varepsilon}^{-\frac{n-2}{2}} V \left( \frac{\varphi_{\varepsilon}(\cdot) - \varphi_{\varepsilon}(x_{\varepsilon})}{\mu_{\varepsilon}} \right) - \mu_{\varepsilon}^{-\frac{n-2}{2}} (V \circ O) \left( \frac{\exp_{x_{\varepsilon}}^{-1}(\cdot)}{\mu_{\varepsilon}} \right) \right] \to 0 $$
for some $O \in O(n)$. We have thus shown that $B = (B_{\varepsilon})_{\varepsilon>0}$ is a family of bubbles centered at $x_{\varepsilon}$ with radius $\mu_{\varepsilon}$ and modeled on $V \in \Sigma$ if and only if

$$ (3.4) \qquad B_{\varepsilon} - \chi(d_g(\cdot, x_{\varepsilon}))\mu_{\varepsilon}^{-\frac{n-2}{2}}(V \circ O)\left(\frac{\exp_{x_{\varepsilon}}^{-1}(\cdot)}{\mu_{\varepsilon}}\right) \to 0 \quad \text{in } H^1(M) $$
there exists a family of charts $(\varphi_{\varepsilon})_{\varepsilon \ge 0}$ such that $\varphi_{\varepsilon} : B_g(x_0, r_g) \to \mathbb{R}^n$ is normal at $x_{\varepsilon}$ for all $\varepsilon > 0$, $\lim_{\varepsilon \to 0} \varphi_{\varepsilon} = \varphi_0$ in $C^k$ for all $k$ and such that

$$ (3.1) \qquad B_{\varepsilon} - \chi(d_g(\cdot, x_{\varepsilon}))\mu_{\varepsilon}^{-\frac{n-2}{2}} V \left( \frac{\varphi_{\varepsilon}(\cdot) - \varphi_{\varepsilon}(x_{\varepsilon})}{\mu_{\varepsilon}} \right) \to 0 \quad \text{in } H^1(M) $$
as $\varepsilon \to 0$.

Here $\Sigma$ is as in (2.1). As the following simple result shows, Definition 3.1 is independent of the choice of the local chart up to an isometry. In the sequel, we let $O(n)$ be the group of isometries of the Euclidean space $\mathbb{R}^n$:

**Proposition 3.1.** Let $(x_{\varepsilon})_{\varepsilon>0} \in M$ and $(\mu_{\varepsilon})_{\varepsilon>0} \in \mathbb{R}_{>0}$ be families of points in $M$ and positive real numbers, with $\mu_{\varepsilon} \to 0$ as $\varepsilon \to 0$, and let $B = (B_{\varepsilon})_{\varepsilon}$ and $\bar{B} = (\bar{B}_{\varepsilon})_{\varepsilon}$ be families of bubbles centered at $x_{\varepsilon}$ with radius $\mu_{\varepsilon}$ and respectively modeled on $V, \bar{V} \in \Sigma$ as in Definition 3.1. Then, up to extraction,

$$ B_{\varepsilon} - \bar{B}_{\varepsilon} \to 0 \text{ in } H^1(M) \text{ as } \varepsilon \to 0 \iff \bar{V} = V \circ O $$

for some orthogonal mapping $O \in O(n)$.

Families of bubbles $(B_{\varepsilon})_{\varepsilon>0}$ and $(\bar{B}_{\varepsilon})_{\varepsilon>0}$ centered at $x_{\varepsilon}$ with radius $\mu_{\varepsilon}$ and satisfying $B_{\varepsilon} - \bar{B}_{\varepsilon} \to 0$ in $H^1(M)$ as $\varepsilon \to 0$ will be called equivalent families.

*Proof.* Let $\varphi_{\varepsilon}, \bar{\varphi}_{\varepsilon} : B_g(x_0, r_g) \to \mathbb{R}^n$ be normal coordinate charts for $g$ at $x_{\varepsilon}$ associated to $B$ and $\bar{B}$ as in definition 3.1. If $x \in \mathbb{R}^n$ we have

$$ (3.2) \qquad \varphi_{\varepsilon} \circ \bar{\varphi}_{\varepsilon}^{-1}(\bar{\varphi}_{\varepsilon}(x_{\varepsilon}) + x) = \varphi_{\varepsilon}(x_{\varepsilon}) + O_{\varepsilon}x + O(|x|^2) \quad \text{as } x \to 0, $$

where we have let $O_{\varepsilon} = D\varphi_{\varepsilon}(x_{\varepsilon}) \circ (D\bar{\varphi}_{\varepsilon}(x_{\varepsilon}))^{-1}$. Since $\varphi_{\varepsilon}$ and $\bar{\varphi}_{\varepsilon}$ are normal charts at $x_{\varepsilon}$, $O_{\varepsilon} \in O(n)$. Straightforward computations using (2.8) and (3.2) show that

$$ (3.3) \qquad \chi(d_g(\cdot, x_{\varepsilon})) \left[ \mu_{\varepsilon}^{-\frac{n-2}{2}} V \left( \frac{\varphi_{\varepsilon}(\cdot) - \varphi_{\varepsilon}(x_{\varepsilon})}{\mu_{\varepsilon}} \right) - \mu_{\varepsilon}^{-\frac{n-2}{2}} V \left( \frac{O_{\varepsilon}(\bar{\varphi}_{\varepsilon}(\cdot) - \bar{\varphi}_{\varepsilon}(x_{\varepsilon}))}{\mu_{\varepsilon}} \right) \right] \to 0 $$
in $H^1(M)$ as $\varepsilon \to 0$. Since $B_{\varepsilon} - \bar{B}_{\varepsilon} \to 0$ in $H^1(M)$ by assumption we obtain

$$ \chi(d_g(\cdot, x_{\varepsilon})) \left[ \mu_{\varepsilon}^{-\frac{n-2}{2}} V \left( \frac{O_{\varepsilon}(\bar{\varphi}_{\varepsilon}(\cdot) - \bar{\varphi}_{\varepsilon}(x_{\varepsilon}))}{\mu_{\varepsilon}} \right) - \mu_{\varepsilon}^{-\frac{n-2}{2}} \bar{V} \left( \frac{\bar{\varphi}_{\varepsilon}(\cdot) - \bar{\varphi}_{\varepsilon}(x_{\varepsilon})}{\mu_{\varepsilon}} \right) - \right] \to 0 $$
in $H^1(M)$ as $\varepsilon \to 0$. A simple scaling argument then shows that $V \circ O = \bar{V}$ where $O := \lim_{\varepsilon \to 0} O_{\varepsilon}$ up to extraction due to the compactness of $O(n)$. $\square$

Equations (3.2) and (3.3) show that, up to composing $V$ with an euclidean isometry, we can always assume that $\varphi_{\varepsilon}$ is an exponential chart at $x_{\varepsilon}$ in Definition 3.1. Choose, for any $p \in B_g(x_0, r_g)$, an identification of $T_pM$ with $\mathbb{R}^n$. This defines an exponential map at $p$, that we will denote by $\exp_p : \mathbb{R}^n \to M$ in all of this paper. Mimicking the arguments in (3.2) and (3.3) we similarly have, for $V \in \Sigma$,

$$ \chi(d_g(\cdot, x_{\varepsilon})) \left[ \mu_{\varepsilon}^{-\frac{n-2}{2}} V \left( \frac{\varphi_{\varepsilon}(\cdot) - \varphi_{\varepsilon}(x_{\varepsilon})}{\mu_{\varepsilon}} \right) - \mu_{\varepsilon}^{-\frac{n-2}{2}} (V \circ O) \left( \frac{\exp_{x_{\varepsilon}}^{-1}(\cdot)}{\mu_{\varepsilon}} \right) \right] \to 0 $$

for some $O \in O(n)$. We have thus shown that $B = (B_{\varepsilon})_{\varepsilon>0}$ is a family of bubbles centered at $x_{\varepsilon}$ with radius $\mu_{\varepsilon}$ and modeled on $V \in \Sigma$ if and only if

$$ (3.4) \qquad B_{\varepsilon} - \chi(d_g(\cdot, x_{\varepsilon}))\mu_{\varepsilon}^{-\frac{n-2}{2}}(V \circ O)\left(\frac{\exp_{x_{\varepsilon}}^{-1}(\cdot)}{\mu_{\varepsilon}}\right) \to 0 \quad \text{in } H^1(M) $$
as $\epsilon \to 0$, for some $O \in O(n)$. Any two exponential charts at $x_\epsilon$ are related by an isometry of $\mathbb{R}^n$. A consequence of (3.4) is that Definition 3.1 is independent of the choice of the exponential chart at $x_\epsilon$ up to composing $V$ with an isometry. The classical example of a family of bubbles is the so-called family of standard bubbles:

$$B_{\epsilon}(x) := \chi(d_{g}(x, x_{\epsilon}))B_{\epsilon}^{+}(x), \text{ where} \\ (3.5) \qquad B_{\epsilon}^{+}(x) := \left( \frac{\mu_{\epsilon}}{\mu_{\epsilon}^{2} + \frac{d_{g}(x, x_{\epsilon})^{2}}{n(n-2)}} \right)^{\frac{n-2}{2}} = \mu_{\epsilon}^{-\frac{n-2}{2}} B_{0}^{+} \left( \frac{\exp_{x_{\epsilon}}^{-1}(x)}{\mu_{\epsilon}} \right),$$

where $B_0^+$ is given by (2.2). This example is peculiar: since $B_0^+$ is radially symmetrical (as a consequence of positivity), any exponential chart at $x_\epsilon$ yields the same family of bubbles, and the function $B_\epsilon^+$ thus obtained is well-defined on the whole of $M$. Naively composing with an exponential chart does not, however, allow to globally pull-back a sign-changing element $V \in \Sigma$ in $M$ without losing regularity at the cut-locus, hence the need for Definition 3.1.

The following proposition is straightforward:

**Proposition 3.2.** Let $(B_\epsilon)_{\epsilon>0}$ be a family of bubble as in Definition 3.1. Then

* There exists $C > 0$ such that $\|B_\epsilon\|_{H^1} \le C$ for all $\epsilon > 0$

* $\lim_{\epsilon \to 0} \|B_\epsilon\|_{L^2} = 0$

* $B_\epsilon \to 0$ weakly in $H^1(M)$ as $\epsilon \to 0$

* $\lim_{\epsilon \to 0} \|\nabla B_\epsilon\|_2^2 = \int_{\mathbb{R}^n} |\nabla V|^2 dx$.

3.2. The product $Weyl_g \otimes B$. We now investigate the properties of the term $Weyl_g \otimes B$ appearing in the Statement of Theorem 1.1.

**Definition 3.2.** Assume that $n \ge 5$. Let $B = (B_\epsilon)_{\epsilon>0}$ be a family of bubbles centered at $x_\epsilon$ with radius $\mu_\epsilon$ and modeled on $V \in \Sigma$ as in Definition 3.1. We define

$$ (3.6) \qquad \begin{aligned} Weyl_g \otimes B := & \frac{4n}{3(n-2)^2} \left( \int_{\mathbb{R}^n} |V|^{2^*} dx \right)^{-1} \times \\ & Weyl_g(x_0)_{i\alpha j\beta} \int_{\mathbb{R}^n} x^\alpha x^\beta \partial_{ij}^2 V \left( \frac{n-2}{2}V + x^l \partial_l V \right) dx \end{aligned} $$

where $x_0 = \lim_{\epsilon \to 0} x_\epsilon$ and where the coordinates of the Weyl tensor are taken with respect to the chart $\exp_{x_0}$.

It is intended in (3.6) that repeated indices are summed over. That the second integral in (3.6) is finite follows from (2.8) since $n \ge 5$. Note that the term $Weyl_g \otimes B$ depends on $x_0$ since the family $B = (B_\epsilon)_\epsilon > 0$ does. We recall the symmetries of the Weyl tensor of $g$ at any point $x \in M$:

$$ Weyl_g(x)_{ijk\ell} = -Weyl_g(x)_{jik\ell} = -Weyl_g(x)_{ij\ell k} = Weyl_g(x)_{k\ell ij} $$

$$ Weyl_g(x)_{iji\ell} = 0 \text{ and } Weyl_g(x)_{ijk\ell} + Weyl_g(x)_{ji\ell k} + Weyl_g(x)_{kij\ell} = 0 $$
as $\epsilon \to 0$, for some $O \in O(n)$. Any two exponential charts at $x_\epsilon$ are related by an isometry of $\mathbb{R}^n$. A consequence of (3.4) is that Definition 3.1 is independent of the choice of the exponential chart at $x_\epsilon$ up to composing $V$ with an isometry. The classical example of a family of bubbles is the so-called family of standard bubbles:

$$B_{\epsilon}(x) := \chi(d_{g}(x, x_{\epsilon}))B_{\epsilon}^{+}(x), \text{ where} \\ (3.5) \qquad B_{\epsilon}^{+}(x) := \left( \frac{\mu_{\epsilon}}{\mu_{\epsilon}^{2} + \frac{d_{g}(x, x_{\epsilon})^{2}}{n(n-2)}} \right)^{\frac{n-2}{2}} = \mu_{\epsilon}^{-\frac{n-2}{2}} B_{0}^{+} \left( \frac{\exp_{x_{\epsilon}}^{-1}(x)}{\mu_{\epsilon}} \right),$$

where $B_0^+$ is given by (2.2). This example is peculiar: since $B_0^+$ is radially symmetrical (as a consequence of positivity), any exponential chart at $x_\epsilon$ yields the same family of bubbles, and the function $B_\epsilon^+$ thus obtained is well-defined on the whole of $M$. Naively composing with an exponential chart does not, however, allow to globally pull-back a sign-changing element $V \in \Sigma$ in $M$ without losing regularity at the cut-locus, hence the need for Definition 3.1.

The following proposition is straightforward:

**Proposition 3.2.** Let $(B_\epsilon)_{\epsilon>0}$ be a family of bubble as in Definition 3.1. Then

* There exists $C > 0$ such that $\|B_\epsilon\|_{H^1} \le C$ for all $\epsilon > 0$

* $\lim_{\epsilon \to 0} \|B_\epsilon\|_{L^2} = 0$

* $B_\epsilon \to 0$ weakly in $H^1(M)$ as $\epsilon \to 0$

* $\lim_{\epsilon \to 0} \|\nabla B_\epsilon\|_2^2 = \int_{\mathbb{R}^n} |\nabla V|^2 dx$.

3.2. The product $Weyl_g \otimes B$. We now investigate the properties of the term $Weyl_g \otimes B$ appearing in the Statement of Theorem 1.1.

**Definition 3.2.** Assume that $n \ge 5$. Let $B = (B_\epsilon)_{\epsilon>0}$ be a family of bubbles centered at $x_\epsilon$ with radius $\mu_\epsilon$ and modeled on $V \in \Sigma$ as in Definition 3.1. We define

$$ (3.6) \qquad \begin{aligned} Weyl_g \otimes B := & \frac{4n}{3(n-2)^2} \left( \int_{\mathbb{R}^n} |V|^{2^*} dx \right)^{-1} \times \\ & Weyl_g(x_0)_{i\alpha j\beta} \int_{\mathbb{R}^n} x^\alpha x^\beta \partial_{ij}^2 V \left( \frac{n-2}{2}V + x^l \partial_l V \right) dx \end{aligned} $$

where $x_0 = \lim_{\epsilon \to 0} x_\epsilon$ and where the coordinates of the Weyl tensor are taken with respect to the chart $\exp_{x_0}$.

It is intended in (3.6) that repeated indices are summed over. That the second integral in (3.6) is finite follows from (2.8) since $n \ge 5$. Note that the term $Weyl_g \otimes B$ depends on $x_0$ since the family $B = (B_\epsilon)_\epsilon > 0$ does. We recall the symmetries of the Weyl tensor of $g$ at any point $x \in M$:

$$ Weyl_g(x)_{ijk\ell} = -Weyl_g(x)_{jik\ell} = -Weyl_g(x)_{ij\ell k} = Weyl_g(x)_{k\ell ij} $$

$$ Weyl_g(x)_{iji\ell} = 0 \text{ and } Weyl_g(x)_{ijk\ell} + Weyl_g(x)_{ji\ell k} + Weyl_g(x)_{kij\ell} = 0 $$
for any $i, j, k, \ell \in \{1, \dots, n\}$, where repeated indices are summed over. As a consequence of these symmetries, simple integration by parts (see for instance Hebey-Vaugon [17] or Lemma 3.2 in Mesmar-Robert [24]) show that

$$ (3.7) \qquad \begin{aligned} \mathrm{Weyl}_g \otimes B = & \frac{4n}{3(n-2)^2} \left( \int_{\mathbb{R}^n} |V|^{2^*} dx \right)^{-1} \times \\ & \mathrm{Weyl}_g(x_0)_{i\alpha j\beta} \int_{\mathbb{R}^n} x^\alpha x^\beta \partial_i V \partial_j V dx. \end{aligned} $$

Again the integral on the right is finite when $n \ge 5$ by (2.8). A consequence of (3.7) and of the tensorial nature of $\mathrm{Weyl}_g(x_0)$ is that $\mathrm{Weyl}_g \otimes B$ does not change if $V$ is replaced by $V \circ O$ for some $O \in O(n)$, provided the coordinates of $\mathrm{Weyl}_g(x_0)$ are now computed in the chart $O \circ \exp_{x_0}^{-1}$. Hence $\mathrm{Weyl}_g \otimes B$ does not depend on the choice of the exponential chart at $x_0$ and only depends on the equivalence class of families of bubbles at $x_0$ in the sense of Proposition 3.1.

We describe two simple cases where $\mathrm{Weyl}_g \otimes B$ vanishes globally:

**Proposition 3.3.** Assume that $n \ge 5$ and let $B$ be a family of bubbles centered at $(x_\varepsilon)_\varepsilon \in M$, with radius $(\mu_\varepsilon)_\varepsilon \in \mathbb{R}_{>0}$ and modeled on $V \in \Sigma$, and $x_0 = \lim_{\varepsilon \to 0} x_\varepsilon$. Then $\mathrm{Weyl}_g \otimes B = 0$ in the following (non-exhaustive) situations:

* $\mathrm{Weyl}_g(x_0) = 0$

* $V$ is radial (and hence $V = B_0^+$ given by (2.2)).

*Proof.* By (3.6) the first point is trivial. If $V = B_0^+$ is radial, the second point follows from the radiality of $B_0^+$, from the equivalent expression (3.7) and from the total antisymmetry of the Weyl tensor. $\square$

The next two results shows that $\mathrm{Weyl}_g \otimes B$ cannot be expected to vanish in general:

**Proposition 3.4.** Let $(\mathbb{S}^p, g_p)$ and $(\mathbb{S}^q, g_q)$ be the unit spheres of dimensions $p, q \ge 2$ endowed with their round metrics. Assume that $n = p + q \ge 5$ and that $V \in D^{1,2}(\mathbb{R}^n) \setminus \{0\}$ solves (1.2) and is invariant under the action of $O(p) \times O(q)$ but not under the action of $O(n)$. Let $M = \mathbb{S}^p \times \mathbb{S}^q$ endowed with the metric $g = g_p \otimes g_q$. Let $B = (B_\varepsilon)_{\varepsilon>0}$ be a family of bubbles in $(M, g)$ modeled on $V$ as in Definition 3.1, centered at $(x_\varepsilon)_{\varepsilon>0} \in M$ and with radius $(\mu_\varepsilon)_{\varepsilon>0} \in \mathbb{R}_{>0}$ and $x_0 = \lim_{\varepsilon \to 0} x_\varepsilon$. Then $\mathrm{Weyl}_g \otimes B < 0$.

*Proof.* We let $M := \mathbb{S}^p \times \mathbb{S}^q$ where $p, q \ge 2$ and $p + q = n \ge 5$ and $g := g_p \otimes g_q$ be the product metric of the round metrics $g_p$ and $g_q$ respectively on $\mathbb{S}^p$ and $\mathbb{S}^q$. We follow the geometric notations of Besse [4] and Hebey [15]. The Riemann tensors for each metric are

$$ \mathrm{Rm}_{g_p} = \frac{1}{2} g_p \odot g_p \quad \text{and} \quad \mathrm{Rm}_{g_q} = \frac{1}{2} g_q \odot g_q, $$

where $\odot$ is the Kulkarni-Nomizu product, given in coordinates by $(T \odot S)_{ijkl} = T_{ik}S_{jl} - T_{il}S_{jk} + T_{jl}S_{ik} - T_{jk}S_{il}$. Standard results on product manifolds give that, after the necessary projections, $\mathrm{Rm}_g = \frac{1}{2}g_p \odot g_p + \frac{1}{2}g_q \odot g_q$, $\mathrm{Ric}_g = (p-1)g_p + (q-1)g_q$ and $S_g = p(p-1) + q(q-1)$. The Weyl tensor of $g$ is defined as

$$ \mathrm{Weyl}_g = \mathrm{Rm}_g - \frac{1}{n-2} \mathrm{Ric}_g \odot g + \frac{S_g}{2(n-1)(n-2)} g \odot g. $$
for any $i, j, k, \ell \in \{1, \dots, n\}$, where repeated indices are summed over. As a consequence of these symmetries, simple integration by parts (see for instance Hebey-Vaugon [17] or Lemma 3.2 in Mesmar-Robert [24]) show that

$$ (3.7) \qquad \begin{aligned} \mathrm{Weyl}_g \otimes B &= \frac{4n}{3(n-2)^2} \left( \int_{\mathbb{R}^n} |V|^{2^*} dx \right)^{-1} \times \\ & \qquad \mathrm{Weyl}_g(x_0)_{i\alpha j\beta} \int_{\mathbb{R}^n} x^\alpha x^\beta \partial_i V \partial_j V dx. \end{aligned} $$

Again the integral on the right is finite when $n \ge 5$ by (2.8). A consequence of (3.7) and of the tensorial nature of $\mathrm{Weyl}_g(x_0)$ is that $\mathrm{Weyl}_g \otimes B$ does not change if $V$ is replaced by $V \circ O$ for some $O \in O(n)$, provided the coordinates of $\mathrm{Weyl}_g(x_0)$ are now computed in the chart $O \circ \exp_{x_0}^{-1}$. Hence $\mathrm{Weyl}_g \otimes B$ does not depend on the choice of the exponential chart at $x_0$ and only depends on the equivalence class of families of bubbles at $x_0$ in the sense of Proposition 3.1.

We describe two simple cases where $\mathrm{Weyl}_g \otimes B$ vanishes globally:

**Proposition 3.3.** Assume that $n \ge 5$ and let $B$ be a family of bubbles centered at $(x_\varepsilon)_\varepsilon \in M$, with radius $(\mu_\varepsilon)_\varepsilon \in \mathbb{R}_{>0}$ and modeled on $V \in \Sigma$, and $x_0 = \lim_{\varepsilon \to 0} x_\varepsilon$. Then $\mathrm{Weyl}_g \otimes B = 0$ in the following (non-exhaustive) situations:

* $\mathrm{Weyl}_g(x_0) = 0$

* $V$ is radial (and hence $V = B_0^+$ given by (2.2)).

*Proof.* By (3.6) the first point is trivial. If $V = B_0^+$ is radial, the second point follows from the radiality of $B_0^+$, from the equivalent expression (3.7) and from the total antisymmetry of the Weyl tensor. $\square$

The next two results shows that $\mathrm{Weyl}_g \otimes B$ cannot be expected to vanish in general:

**Proposition 3.4.** Let $(\mathbb{S}^p, g_p)$ and $(\mathbb{S}^q, g_q)$ be the unit spheres of dimensions $p, q \ge 2$ endowed with their round metrics. Assume that $n = p + q \ge 5$ and that $V \in D^{1,2}(\mathbb{R}^n) \setminus \{0\}$ solves (1.2) and is invariant under the action of $O(p) \times O(q)$ but not under the action of $O(n)$. Let $M = \mathbb{S}^p \times \mathbb{S}^q$ endowed with the metric $g = g_p \otimes g_q$. Let $B = (B_\varepsilon)_{\varepsilon>0}$ be a family of bubbles in $(M, g)$ modeled on $V$ as in Definition 3.1, centered at $(x_\varepsilon)_{\varepsilon>0} \in M$ and with radius $(\mu_\varepsilon)_{\varepsilon>0} \in \mathbb{R}_{>0}$ and $x_0 = \lim_{\varepsilon \to 0} x_\varepsilon$. Then $\mathrm{Weyl}_g \otimes B < 0$.

*Proof.* We let $M := \mathbb{S}^p \times \mathbb{S}^q$ where $p, q \ge 2$ and $p + q = n \ge 5$ and $g := g_p \otimes g_q$ be the product metric of the round metrics $g_p$ and $g_q$ respectively on $\mathbb{S}^p$ and $\mathbb{S}^q$. We follow the geometric notations of Besse [4] and Hebey [15]. The Riemann tensors for each metric are

$$ \mathrm{Rm}_{g_p} = \frac{1}{2} g_p \odot g_p \quad \text{and} \quad \mathrm{Rm}_{g_q} = \frac{1}{2} g_q \odot g_q, $$

where $\odot$ is the Kulkarni-Nomizu product, given in coordinates by $(T \odot S)_{ijkl} = T_{ik}S_{jl} - T_{il}S_{jk} + T_{jl}S_{ik} - T_{jk}S_{il}$. Standard results on product manifolds give that, after the necessary projections, $\mathrm{Rm}_g = \frac{1}{2}g_p \odot g_p + \frac{1}{2}g_q \odot g_q$, $\mathrm{Ric}_g = (p-1)g_p + (q-1)g_q$ and $S_g = p(p-1) + q(q-1)$. The Weyl tensor of $g$ is defined as

$$ \mathrm{Weyl}_g = \mathrm{Rm}_g - \frac{1}{n-2} \mathrm{Ric}_g \odot g + \frac{S_g}{2(n-1)(n-2)} g \odot g. $$
If $x_0 = (x_0^p, x_0^q) \in \mathbb{S}^p \times \mathbb{S}^q$ we compute the coordinates of $\text{Weyl}_g(x_0)$ in the chart obtained from the exponential chart on $\mathbb{S}^p$ at $x_0^p$ and the exponential chart on $\mathbb{S}^q$ at $x_0^q$; this yields a normal chart on $\mathbb{S}^p \times \mathbb{S}^q$. The Latin letters $i, j, k...$ refer to coordinates in $\mathbb{S}^p$, the Greek letter $\alpha, \beta, ...$ refer to $\mathbb{S}^q$. The previous observations show that, at any point around $x_0$ we have

$$ (\text{Weyl}_g)_{ijkl} = C_1(g_{ik}g_{jl} - g_{il}g_{jk}) $$

$$ (\text{Weyl}_g)_{\alpha\beta\gamma\delta} = C_2(g_{\alpha\gamma}g_{\beta\delta} - g_{\alpha\delta}g_{\beta\gamma}) $$

$$ (\text{Weyl}_g)_{i\alpha j\beta} = -C_3 g_{ij} g_{\alpha\beta} $$

$$ (\text{Weyl}_g)_{ij\alpha\beta} = 0, $$

with $C_1 = \frac{2q(q-1)}{(n-1)(n-2)}$, $C_2 = \frac{2p(p-1)}{(n-1)(n-2)}$ and $C_3 = \frac{2(p-1)(q-1)}{(n-1)(n-2)}$. The mixed terms $(\text{Weyl}_g)_{ijk\alpha}$ and $(\text{Weyl}_g)_{i\alpha\beta\gamma}$ are null, and the other terms are obtained via the symmetries of the Weyl tensor. We denote the coordinates in $\mathbb{R}^n$ by $(X, Y)$, where $X \in \mathbb{R}^p$ and $Y \in \mathbb{R}^q$. Let $V \in \Sigma$ be $O(p) \times O(q)$ invariant and let $C_{n,V} = \frac{4n}{3(n-2)^2} (\int_{\mathbb{R}^n} |V|^{2*} dx)^{-1}$. Using (3.7) and the explicit expression of $\text{Weyl}_g(x_0)$ we have

$$ C_{n,V}^{-1} \text{Weyl}_g \otimes B = C_1 \int_{\mathbb{R}^n} [|X|^2 |\nabla_1 V|^2 - (X, \nabla_1 V)^2] dX dY $$

$$ + C_2 \int_{\mathbb{R}^n} [|Y|^2 |\nabla_2 V|^2 - (Y, \nabla_2 V)^2] dX dY $$

$$ - C_3 \int_{\mathbb{R}^n} [|X|^2 |\nabla_2 V|^2 + |Y|^2 |\nabla_1 V|^2 - 2(X, \nabla_1 V)(Y, \nabla_2 V)] dX dY $$

where $\nabla_1 V$ and $\nabla_2 V$ denote respectively the gradient of $V$ with respect to $X \in \mathbb{R}^p$ and $Y \in \mathbb{R}^q$. The $O(p) \times O(q)$ invariance of $V$ implies that $V(X, Y) = V(|X|, |Y|)$, and as a consequence we have $\nabla_1 V(X, Y) = \partial_1 V(X, Y)_{\frac{X}{|X|}}$ and $\nabla_2 V(X, Y) = \partial_2 V(X, Y)_{\frac{Y}{|Y|}}$ where $\partial_1 V$ and $\partial_2 V$ are the partial derivatives of $(r_1, r_2) \mapsto V(r_1, r_2)$. We obtain in the end

$$ C_{n,V}^{-1} \text{Weyl}_g \otimes B = -C_3 \int_{\mathbb{R}^n} [|X|^2 (\partial_2 V)^2 + |Y|^2 (\partial_1 V)^2 - 2|X||Y|\partial_1 V \partial_2 V] dX dY $$

$$ = -C_3 \int_{\mathbb{R}^n} (|X|\partial_2 V - |Y|\partial_1 V)^2 dX dY \le 0. $$

If $\text{Weyl}_g \otimes B = 0$, then $r_2\partial_1 V - r_1\partial_2 V$ for all $r_1 = |X| > 0$ and $r_2 = |Y| > 0$. This shows that $V = V(r)$, for $r^2 = r_1^2 + r_2^2$ when $r_1, r_2 > 0$, and this equality extends continuously to $\mathbb{R}^n \setminus \{0\}$. We then get that $V$ is radially symmetrical, contradicting our assumption. Therefore $\text{Weyl}_g \otimes B < 0$. $\square$

**Corollary 3.5.** There are manifolds $(M, g)$ of dimension $n \ge 5$ and families of bubbles $B = (B_\varepsilon)_{\varepsilon>0}$ as in Definition 3.1 such that $\text{Weyl}_g \otimes B < 0$.

*Proof.* We choose $p, q \ge 2$ such that $n := p + q \ge 5$. As before we denote by $g_n$ the round metric in $\mathbb{S}^n$. It follows from the result of Ding [9] that for all $L > 0$, there exists $\tilde{V} \in C^2(\mathbb{S}^n)$ that solves $\Delta_{g_n} \tilde{V} + \frac{n(n-2)}{4} \tilde{V} = |\tilde{V}|^{2*-2} \tilde{V}$, satisfies $\int_{\mathbb{S}^n} |\tilde{V}|^{2*} dv_{g_n} \ge L$ and is $O(p) \times O(q+1)$-invariant, that is $\tilde{V}(x, y) = \tilde{V}(|x|, |y|)$ for $(x, y) \in \mathbb{R}^p \times \mathbb{R}^{q+1} = \mathbb{R}^{n+1}$. Let $\pi_N : \mathbb{S}^n \setminus \{N\} \to \mathbb{R}^n$ be the stereographic
If $x_0 = (x_0^p, x_0^q) \in \mathbb{S}^p \times \mathbb{S}^q$ we compute the coordinates of $\text{Weyl}_g(x_0)$ in the chart obtained from the exponential chart on $\mathbb{S}^p$ at $x_0^p$ and the exponential chart on $\mathbb{S}^q$ at $x_0^q$: this yields a normal chart on $\mathbb{S}^p \times \mathbb{S}^q$. The Latin letters $i, j, k...$ refer to coordinates in $\mathbb{S}^p$, the Greek letter $\alpha, \beta, ...$ refer to $\mathbb{S}^q$. The previous observations show that, at any point around $x_0$ we have

$$ (\text{Weyl}_g)_{ijkl} = C_1(g_{ik}g_{jl} - g_{il}g_{jk}) $$

$$ (\text{Weyl}_g)_{\alpha\beta\gamma\delta} = C_2(g_{\alpha\gamma}g_{\beta\delta} - g_{\alpha\delta}g_{\beta\gamma}) $$

$$ (\text{Weyl}_g)_{i\alpha j\beta} = -C_3 g_{ij} g_{\alpha\beta} $$

$$ (\text{Weyl}_g)_{ij\alpha\beta} = 0, $$

with $C_1 = \frac{2q(q-1)}{(n-1)(n-2)}$, $C_2 = \frac{2p(p-1)}{(n-1)(n-2)}$ and $C_3 = \frac{2(p-1)(q-1)}{(n-1)(n-2)}$. The mixed terms $(\text{Weyl}_g)_{ijk\alpha}$ and $(\text{Weyl}_g)_{i\alpha\beta\gamma}$ are null, and the other terms are obtained via the symmetries of the Weyl tensor. We denote the coordinates in $\mathbb{R}^n$ by $(X, Y)$, where $X \in \mathbb{R}^p$ and $Y \in \mathbb{R}^q$. Let $V \in \Sigma$ be $O(p) \times O(q)$ invariant and let $C_{n,V} = \frac{4n}{3(n-2)^2} (\int_{\mathbb{R}^n} |V|^{2^*} dx)^{-1}$. Using (3.7) and the explicit expression of $\text{Weyl}_g(x_0)$ we have

$$ C_{n,V}^{-1} \text{Weyl}_g \otimes B = C_1 \int_{\mathbb{R}^n} [|X|^2 |\nabla_1 V|^2 - (X, \nabla_1 V)^2] dX dY $$

$$ + C_2 \int_{\mathbb{R}^n} [|Y|^2 |\nabla_2 V|^2 - (Y, \nabla_2 V)^2] dX dY $$

$$ - C_3 \int_{\mathbb{R}^n} [|X|^2 |\nabla_2 V|^2 + |Y|^2 |\nabla_1 V|^2 - 2(X, \nabla_1 V)(Y, \nabla_2 V)] dX dY $$

where $\nabla_1 V$ and $\nabla_2 V$ denote respectively the gradient of $V$ with respect to $X \in \mathbb{R}^p$ and $Y \in \mathbb{R}^q$. The $O(p) \times O(q)$ invariance of $V$ implies that $V(X, Y) = V(|X|, |Y|)$, and as a consequence we have $\nabla_1 V(X, Y) = \partial_1 V(X, Y)_{\frac{X}{|X|}}$ and $\nabla_2 V(X, Y) = \partial_2 V(X, Y)_{\frac{Y}{|Y|}}$ where $\partial_1 V$ and $\partial_2 V$ are the partial derivatives of $(r_1, r_2) \mapsto V(r_1, r_2)$. We obtain in the end

$$ C_{n,V}^{-1} \text{Weyl}_g \otimes B = -C_3 \int_{\mathbb{R}^n} [|X|^2 (\partial_2 V)^2 + |Y|^2 (\partial_1 V)^2 - 2|X||Y|\partial_1 V \partial_2 V] dX dY $$

$$ = -C_3 \int_{\mathbb{R}^n} (|X|\partial_2 V - |Y|\partial_1 V)^2 dX dY \le 0. $$

If $\text{Weyl}_g \otimes B = 0$, then $r_2\partial_1 V - r_1\partial_2 V$ for all $r_1 = |X| > 0$ and $r_2 = |Y| > 0$. This shows that $V = V(r)$, for $r^2 = r_1^2 + r_2^2$ when $r_1, r_2 > 0$, and this equality extends continuously to $\mathbb{R}^n \setminus \{0\}$. We then get that $V$ is radially symmetrical, contradicting our assumption. Therefore $\text{Weyl}_g \otimes B < 0$. $\square$

**Corollary 3.5.** There are Riemannian manifolds $(M, g)$ of dimension $n \ge 5$ and families of bubbles $B = (B_\varepsilon)_{\varepsilon>0}$ as in Definition 3.1 such that $\text{Weyl}_g \otimes B < 0$.

*Proof.* We choose $p, q \ge 2$ such that $n := p + q \ge 5$. As before we denote by $g_n$ the round metric in $\mathbb{S}^n$. It follows from the result of Ding [9] that for all $L > 0$, there exists $\tilde{V} \in C^2(\mathbb{S}^n)$ that solves $\Delta_{g_n} \tilde{V} + \frac{n(n-2)}{4} \tilde{V} = |\tilde{V}|^{2^*-2} \tilde{V}$, satisfies $\int_{\mathbb{S}^n} |\tilde{V}|^{2^*} dv_{g_n} \ge L$ and is $O(p) \times O(q+1)$-invariant, that is $\tilde{V}(x, y) = \tilde{V}(|x|, |y|)$ for $(x, y) \in \mathbb{R}^p \times \mathbb{R}^{q+1} = \mathbb{R}^{n+1}$. Let $\pi_N : \mathbb{S}^n \setminus \{N\} \to \mathbb{R}^n$ be the stereographic
projection from the North pole. Let, for $x \in \mathbb{R}^n$,

$$V(x) := \left( \frac{2}{1 + |x|^2} \right)^{\frac{n-2}{2}} \tilde{V}(\pi_N^{-1}(x)).$$

Then $V \in C^2(\mathbb{R}^n) \cap D^{1,2}(\mathbb{R}^n)$ solves $\Delta_{\xi}V = |V|^{2^*-2}V$ in $\mathbb{R}^n$, satisfies $\int_{\mathbb{R}^n} |V|^{2^*} dx \ge L$ and using the coordinate expression of the stereographic projection it is easily seen that it is $O(p) \times O(q)$-invariant, that is $V(x', x'') = V(|x'|, |x''|)$ for $(x', x'') \in \mathbb{R}^p \times \mathbb{R}^q = \mathbb{R}^n$. Up to taking $L$ large enough, we get that $V$ is not $O(n)$-invariant since by Obata [26] finite-energy radial solutions of (1.2) are translations and rescalings of $B_0^+$ given by (2.2) and all have fixed minimal Dirichlet energy. Let now $M := S^p \times S^q$ and $g := g_p \otimes g_q$ be the product metric of the round metrics on $S^p$ and $S^q$. Let $B = (B_{\varepsilon})_{\varepsilon>0}$ be any family of bubbles in $(M, g)$ modeled on $V$ as in Definition 3.1, centered at $(x_{\varepsilon})_{\varepsilon} \in M$ and with radius $(\mu_{\varepsilon})_{\varepsilon} \in \mathbb{R}_{>0}$ and $x_0 = \lim_{\varepsilon \to 0} x_{\varepsilon}$, where $x_0 \in M$ is any point. For $L$ large enough Proposition 3.4 applies and shows that Weyl$_g \otimes B < 0$. $\square$

3.3. **The mass of the Green's function in 3D.** We conclude this section by recalling the definition of the mass when $n=3$, following Li-Zhu [20]:

**Definition 3.3.** Let $(M, g)$ be a smooth, closed, connected 3-dimensional Riemannian manifold. Let $h \in C^1(M)$ be such that $\Delta_g + h$ is coercive. Let $G_h$ be the Green's function for $\Delta_g + h$. Then for all $x_0 \in M$, there exists $m_h(x_0) \in \mathbb{R}$, denoted as the mass, such that

$$G_h(x, x_0) = \frac{1}{4\pi d_g(x, x_0)} + m_h(x_0) + o(1) \text{ as } x \to x_0.$$

When $n=3$, $\Delta_g + h$ is coercive and $x_0 \in M$, it follows from the construction of the Green's function (see Aubin [3], Robert [32]) that there exists $\beta_{x_0} \in C^0(M) \cap H_2^p(M)$ for all $p < 3$ such that

$$(3.8) \qquad G_h(x, x_0) = \frac{\chi(x)}{4\pi d_g(x, x_0)} + \beta_{x_0}(x) \text{ for all } x \in M \setminus \{x_0\},$$

where $\chi \in C^\infty(M)$ is a cutoff function with $\chi(x) = 1$ in $B_g(x_0, \frac{r_g}{2})$ and $\chi(x) = 0$ in $M \setminus B_g(x_0, r_g)$. We have $(\Delta_g + h)\beta_{x_0} = -(\Delta_g + h)((4\pi)^{-1}\chi d_g(\cdot, x_0)^{-1}) = O(d_g(\cdot, x_0)^{-1})$ and standard elliptic theory then shows that

$$|\beta_{x_0}(x)| \le C, \quad |\nabla \beta_{x_0}(x)| \le C(1 + |\ln d_g(x_0, x)|) \quad \text{for all } x \in B_g\left(x_0, \frac{r_g}{2}\right) \setminus \{0\}.$$

It now follows from Definition 3.3 that

$$m_h(x_0) := \beta_{x_0}(x_0).$$

## 4. A PRIORI POINTWISE ESTIMATES FOR FAMILIES OF BUBBLES

4.1. **Statement of the results.** In this section we will assume, as in the introduction, that $(u_{\varepsilon})_{0<\varepsilon\le 1}$ is a family of possibly sign-changing solutions of (1.3) satisfying (1.4):

$$u_{\varepsilon} = cB_{\varepsilon} + o(1) \quad \text{in } H^1(M),$$

for some $c > 0$, where $B_{\varepsilon}$ is a family of bubbles as in Definition 3.1, centered at $x_{\varepsilon}$ with radius $\mu_{\varepsilon}$. Here $(x_{\varepsilon})_{0<\varepsilon\le 1}$ and $(\mu_{\varepsilon})_{0<\varepsilon\le 1}$ are families of points in $M$ and of positive numbers with $\mu_{\varepsilon} \to 0$ as $\varepsilon \to 0$. We let $V \in \Sigma$, where $\Sigma$ is as in (2.1), be
projection from the North pole. Let, for $x \in \mathbb{R}^n$,

$$V(x) := \left( \frac{2}{1 + |x|^2} \right)^{\frac{n-2}{2}} \tilde{V}(\pi_N^{-1}(x)).$$

Then $V \in C^2(\mathbb{R}^n) \cap D^{1,2}(\mathbb{R}^n)$ solves $\Delta_{\xi}V = |V|^{2^*-2}V$ in $\mathbb{R}^n$, satisfies $\int_{\mathbb{R}^n} |V|^{2^*} dx \ge L$ and using the coordinate expression of the stereographic projection it is easily seen that it is $O(p) \times O(q)$-invariant, that is $V(x', x'') = V(|x'|, |x''|)$ for $(x', x'') \in \mathbb{R}^p \times \mathbb{R}^q = \mathbb{R}^n$. Up to taking $L$ large enough, we get that $V$ is not $O(n)$-invariant since by Obata [26] finite-energy radial solutions of (1.2) are translations and rescalings of $B_0^+$ given by (2.2) and all have fixed minimal Dirichlet energy. Let now $M := S^p \times S^q$ and $g := g_p \otimes g_q$ be the product metric of the round metrics on $S^p$ and $S^q$. Let $B = (B_{\varepsilon})_{\varepsilon>0}$ be any family of bubbles in $(M, g)$ modeled on $V$ as in Definition 3.1, centered at $(x_{\varepsilon})_{\varepsilon} \in M$ and with radius $(\mu_{\varepsilon})_{\varepsilon} \in \mathbb{R}_{>0}$ and $x_0 = \lim_{\varepsilon \to 0} x_{\varepsilon}$, where $x_0 \in M$ is any point. For $L$ large enough Proposition 3.4 applies and shows that Weyl$_g \otimes B < 0$. $\square$

### 3.3. The mass of the Green's function in 3D.
We conclude this section by recalling the definition of the mass when $n=3$, following Li-Zhu [20]:

**Definition 3.3.** Let $(M, g)$ be a smooth, closed, connected 3-dimensional Riemannian manifold. Let $h \in C^1(M)$ be such that $\Delta_g + h$ is coercive. Let $G_h$ be the Green's function for $\Delta_g + h$. Then for all $x_0 \in M$, there exists $m_h(x_0) \in \mathbb{R}$, denoted as the mass, such that

$$G_h(x, x_0) = \frac{1}{4\pi d_g(x, x_0)} + m_h(x_0) + o(1) \text{ as } x \to x_0.$$

When $n=3$, $\Delta_g + h$ is coercive and $x_0 \in M$, it follows from the construction of the Green's function (see Aubin [3], Robert [32]) that there exists $\beta_{x_0} \in C^0(M) \cap H_2^p(M)$ for all $p < 3$ such that

$$(3.8) \qquad G_h(x, x_0) = \frac{\chi(x)}{4\pi d_g(x, x_0)} + \beta_{x_0}(x) \text{ for all } x \in M \setminus \{x_0\},$$

where $\chi \in C^\infty(M)$ is a cutoff function with $\chi(x) = 1$ in $B_g(x_0, \frac{r_g}{2})$ and $\chi(x) = 0$ in $M \setminus B_g(x_0, r_g)$. We have $(\Delta_g + h)\beta_{x_0} = -(\Delta_g + h)((4\pi)^{-1}\chi d_g(\cdot, x_0)^{-1}) = O(d_g(\cdot, x_0)^{-1})$ and standard elliptic theory then shows that

$$|\beta_{x_0}(x)| \le C, \quad |\nabla \beta_{x_0}(x)| \le C(1 + |\ln d_g(x_0, x)|) \quad \text{for all } x \in B_g\left(x_0, \frac{r_g}{2}\right) \setminus \{0\}.$$

It now follows from Definition 3.3 that

$$m_h(x_0) := \beta_{x_0}(x_0).$$

## 4. A PRIORI POINTWISE ESTIMATES FOR FAMILIES OF BUBBLES

### 4.1. Statement of the results.
In this section we will assume, as in the introduction, that $(u_{\varepsilon})_{0<\varepsilon\le 1}$ is a family of possibly sign-changing solutions of (1.3) satisfying (1.4):

$$u_{\varepsilon} = cB_{\varepsilon} + o(1) \quad \text{in } H^1(M),$$

for some $c > 0$, where $B_{\varepsilon}$ is a family of bubbles as in Definition 3.1, centered at $x_{\varepsilon}$ with radius $\mu_{\varepsilon}$. Here $(x_{\varepsilon})_{0<\varepsilon\le 1}$ and $(\mu_{\varepsilon})_{0<\varepsilon\le 1}$ are families of points in $M$ and of positive numbers with $\mu_{\varepsilon} \to 0$ as $\varepsilon \to 0$. We let $V \in \Sigma$, where $\Sigma$ is as in (2.1), be
the function that appears in (3.4). Under this assumption $(u_\epsilon)_{0<\epsilon\le 1}$ blows-up with
only one possibly sign-changing bubble and a zero weak limit. Throughout this
paper we will adopt the following conventions: for nonnegative families $(a_\epsilon)_{0<\epsilon\le 1}$
and $(b_\epsilon)_{0<\epsilon\le 1}$ of functions in $M$ (or real numbers), we will write $a_\epsilon \lesssim b_\epsilon$, or
equivalently $a_\epsilon = O(b_\epsilon)$, when there exists a positive constant $C$ such that

$$
a_{\epsilon}(x) \le C b_{\epsilon}(x) \quad \text{for all } 0 < \epsilon \le 1 \text{ and any } x \in M.
$$

Since $\Delta_g + h_\epsilon$ is coercive we let, for any $0 < \epsilon \le 1$, $G_{h_\epsilon}$ be the positive Green's function of $\Delta_g + h_\epsilon$ (see Robert [32]). Since $h_\epsilon \to h$ in $C^1(M)$ as $\epsilon \to 0$ and $\Delta_g + h$ is coercive there exists $C > 0$ such that

$$
(4.1) \quad \left| d_g(x, y)^{n-2} G_{h_\epsilon}(x, y) - \frac{1}{(n-2)\omega_{n-1}} \right| \le C d_g(x, y)
$$

holds true for any $(x, y) \in M \times M \setminus D(M)$ and for any $0 < \epsilon \le 1$, where $D(M) = \{(x, x), x \in M\}$ and where $\omega_{n-1}$ is the area of the standard sphere $\mathbb{S}^{n-1}$. Define, for $x \in M$,

$$
F_{\epsilon}(x_{\epsilon}, x) = (n-2)\omega_{n-1}d_{g}(x_{\epsilon}, x)^{n-2}G_{h_{\epsilon}}(x_{\epsilon}, x).
$$

We define, for any $x \in M$ and any $0 < \epsilon \le 1$:

$$
(4.2) \qquad
\begin{aligned}
B^{\mu_{\epsilon}, x_{\epsilon}}(x) ={}& \chi(d_g(x_{\epsilon}, x)) F_{\epsilon}(x_{\epsilon}, x) \mu_{\epsilon}^{-\frac{n-2}{2}} V\left(\frac{1}{\mu_{\epsilon}} \exp_{x_{\epsilon}}^{-1}(x)\right) \\
& + (1 - \chi(d_g(x_{\epsilon}, x))) (n-2) \omega_{n-1} \lambda(V) \mu_{\epsilon}^{\frac{n-2}{2}} G_{h_{\epsilon}}(x_{\epsilon}, x),
\end{aligned}
$$

where $\lambda(V)$ is as in (2.9) and for $\chi \in C_c^\infty(\mathbb{R})$ with $\chi \equiv 1$ in $[0, \frac{i_g(M)}{4}]$ and $\chi \equiv 0$ in $[\frac{i_g(M)}{2}, +\infty)$, and where $i_g(M)$ is the injectivity radius of $(M, g)$. By (4.1) and properties of the Green's function (see [32]) we have

$$
F_{\epsilon}(x_{\epsilon}, x) = 1 + O(d_{g}(x_{\epsilon}, x)) \quad \text{and} \quad |\nabla_{x} F_{\epsilon}(x_{\epsilon}, x)| = O(1) \quad \text{for } x \in M.
$$

It is therefore easily seen with (3.4) that

$$
B^{\mu_{\epsilon}, x_{\epsilon}} = B_{\epsilon} + o(1) \quad \text{in } H^1(M),
$$

so that $(B^{\mu_{\epsilon}, x_{\epsilon}})_{0<\epsilon\le 1}$ is a family of bubbles which is equivalent to $(B_{\epsilon})_{0<\epsilon\le 1}$ in the
sense of Proposition 3.1. The original assumption (1.4) implies that $u_{\epsilon}$ still satisfies

$$
(4.3) \qquad u_{\epsilon} = c B^{\mu_{\epsilon}, x_{\epsilon}} + o(1) \quad \text{in } H^1(M)
$$

for some $c > 0$. We will thus, in the subsequent analysis, work with the family $(B^{\mu_{\epsilon}, x_{\epsilon}})_{0<\epsilon\le 1}$. If $(y_{\epsilon})_{\epsilon}$ is a family of points in $M$ satisfying $\mu_{\epsilon} \ll d_g(x_{\epsilon}, y_{\epsilon}) \le \frac{i_g(M)}{4}$ we also have, by (2.4), (4.1) and (4.2),

$$
\begin{align}
B^{\mu_{\epsilon}, x_{\epsilon}}(y_{\epsilon}) &= F_{\epsilon}(x_{\epsilon}, y_{\epsilon}) \mu_{\epsilon}^{-\frac{n-2}{2}} V \left( \frac{1}{\mu_{\epsilon}} \exp_{x_{\epsilon}}^{-1}(y_{\epsilon}) \right) \nonumber \\
&= (n-2) \omega_{n-1} \mu_{\epsilon}^{\frac{n-2}{2}} G_{h_{\epsilon}}(x_{\epsilon}, y_{\epsilon}) \left( \lambda(V) + O \left( \frac{\mu_{\epsilon}}{d_g(x_{\epsilon}, y_{\epsilon})} \right) \right). \tag{4.4}
\end{align}
$$

The function $B^{\mu_{\epsilon}, x_{\epsilon}}$ thus interpolates the (rescaled) pull-back of $V$ in $M$ to its
first-order expansion at infinity.

In this section we prove that (4.3) still holds true globally pointwise in $M$, up
to an error term that is of the order of the positive standard bubble $B_{\epsilon}^{+}$ given by
(3.5). The main result of this section is the following:
the function that appears in (3.4). Under this assumption $(u_\epsilon)_{0<\epsilon\le 1}$ blows-up with
only one possibly sign-changing bubble and a zero weak limit. Throughout this
paper we will adopt the following conventions: for nonnegative families $(a_\epsilon)_{0<\epsilon\le 1}$
and $(b_\epsilon)_{0<\epsilon\le 1}$ of functions in $M$ (or real numbers), we will write $a_\epsilon \lesssim b_\epsilon$, or
equivalently $a_\epsilon = O(b_\epsilon)$, when there exists a positive constant $C$ such that

$$
a_{\epsilon}(x) \le C b_{\epsilon}(x) \quad \text{for all } 0 < \epsilon \le 1 \text{ and any } x \in M.
$$

Since $\Delta_g + h_\epsilon$ is coercive we let, for any $0 < \epsilon \le 1$, $G_{h_\epsilon}$ be the positive Green's function of $\Delta_g + h_\epsilon$ (see Robert [32]). Since $h_\epsilon \to h$ in $C^1(M)$ as $\epsilon \to 0$ and $\Delta_g + h$ is coercive there exists $C > 0$ such that

$$
(4.1) \quad \left| d_g(x, y)^{n-2} G_{h_\epsilon}(x, y) - \frac{1}{(n-2)\omega_{n-1}} \right| \le C d_g(x, y)
$$

holds true for any $(x, y) \in M \times M \setminus D(M)$ and for any $0 < \epsilon \le 1$, where $D(M) = \{(x, x), x \in M\}$ and where $\omega_{n-1}$ is the area of the standard sphere $S^{n-1}$. Define, for $x \in M$,

$$
F_{\epsilon}(x_{\epsilon}, x) = (n-2)\omega_{n-1}d_{g}(x_{\epsilon}, x)^{n-2}G_{h_{\epsilon}}(x_{\epsilon}, x).
$$

We define, for any $x \in M$ and any $0 < \epsilon \le 1$:

$$
(4.2) \qquad
\begin{aligned}
B^{\mu_{\epsilon}, x_{\epsilon}}(x) ={}& \chi(d_g(x_{\epsilon}, x)) F_{\epsilon}(x_{\epsilon}, x) \mu_{\epsilon}^{-\frac{n-2}{2}} V\left(\frac{1}{\mu_{\epsilon}} \exp_{x_{\epsilon}}^{-1}(x)\right) \\
& + (1 - \chi(d_g(x_{\epsilon}, x))) (n-2) \omega_{n-1} \lambda(V) \mu_{\epsilon}^{\frac{n-2}{2}} G_{h_{\epsilon}}(x_{\epsilon}, x),
\end{aligned}
$$

where $\lambda(V)$ is as in (2.9) and for $\chi \in C_c^\infty(\mathbb{R})$ with $\chi \equiv 1$ in $[0, \frac{i_g(M)}{4}]$ and $\chi \equiv 0$ in $[\frac{i_g(M)}{2}, +\infty)$, and where $i_g(M)$ is the injectivity radius of $(M, g)$. By (4.1) and properties of the Green's function (see [32]) we have

$$
F_{\epsilon}(x_{\epsilon}, x) = 1 + O(d_{g}(x_{\epsilon}, x)) \quad \text{and} \quad |\nabla_{x} F_{\epsilon}(x_{\epsilon}, x)| = O(1) \quad \text{for } x \in M.
$$

It is therefore easily seen with (3.4) that

$$
B^{\mu_{\epsilon}, x_{\epsilon}} = B_{\epsilon} + o(1) \quad \text{in } H^1(M),
$$

so that $(B^{\mu_{\epsilon}, x_{\epsilon}})_{0<\epsilon\le 1}$ is a family of bubbles which is equivalent to $(B_{\epsilon})_{0<\epsilon\le 1}$ in the
sense of Proposition 3.1. The original assumption (1.4) implies that $u_{\epsilon}$ still satisfies

$$
(4.3) \qquad u_{\epsilon} = c B^{\mu_{\epsilon}, x_{\epsilon}} + o(1) \quad \text{in } H^1(M)
$$

for some $c > 0$. We will thus, in the subsequent analysis, work with the family $(B^{\mu_{\epsilon}, x_{\epsilon}})_{0<\epsilon\le 1}$. If $(y_{\epsilon})_{\epsilon}$ is a family of points in $M$ satisfying $\mu_{\epsilon} \ll d_g(x_{\epsilon}, y_{\epsilon}) \le \frac{i_g(M)}{4}$ we also have, by (2.4), (4.1) and (4.2),

$$
\begin{align}
B^{\mu_{\epsilon}, x_{\epsilon}}(y_{\epsilon}) &= F_{\epsilon}(x_{\epsilon}, y_{\epsilon}) \mu_{\epsilon}^{-\frac{n-2}{2}} V \left( \frac{1}{\mu_{\epsilon}} \exp_{x_{\epsilon}}^{-1}(y_{\epsilon}) \right) \nonumber \\
&= (n-2) \omega_{n-1} \mu_{\epsilon}^{\frac{n-2}{2}} G_{h_{\epsilon}}(x_{\epsilon}, y_{\epsilon}) \left( \lambda(V) + O \left( \frac{\mu_{\epsilon}}{d_g(x_{\epsilon}, y_{\epsilon})} \right) \right). \tag{4.4}
\end{align}
$$

The function $B^{\mu_{\epsilon}, x_{\epsilon}}$ thus interpolates the (rescaled) pull-back of $V$ in $M$ to its
first-order expansion at infinity.

In this section we prove that (4.3) still holds true globally pointwise in $M$, up
to an error term that is of the order of the positive standard bubble $B_{\epsilon}^{+}$ given by
(3.5). The main result of this section is the following:
**Proposition 4.1.** Let $(M^n, g), n \ge 3$ be a smooth, connected and closed manifold, let $(h_\epsilon)_{0<\epsilon\le 1}$ be a family of $C^1$ functions in $M$ converging in $C^1(M)$ towards $h$ and let $(p_\epsilon)_{0<\epsilon\le 1}$ be a family of real numbers satisfying $2<p_\epsilon \le 2^*$ for all $0<\epsilon \le 1$ and such that $\lim_{\epsilon\to 0} p_\epsilon = 2^*$. Assume that $\Delta_g + h$ is coercive. Let $(u_\epsilon)_{0<\epsilon\le 1}$ be a family of solutions of (1.3) satisfying (1.4) for some $c>0$.

Then $c=1$ in (1.4) and there exists a family $(\sigma_\epsilon)_{0<\epsilon\le 1}$ of positive real numbers with $\lim_{\epsilon\to 0} \sigma_\epsilon = 0$ such that,

$$ (4.5) \qquad \left\| \frac{u_\epsilon - B^{\mu_\epsilon, x_\epsilon}}{B_\epsilon^+} \right\|_{L^\infty(M)} \le \sigma_\epsilon $$

for any $0 < \epsilon \le 1$, where $B^{\mu_\epsilon, x_\epsilon}$ is as in (4.2) and $B_\epsilon^+$ is as in (3.5).

The function $V \in \Sigma$ that appears in the definition of $B^{\mu_\epsilon, x_\epsilon}$ is the one that defines the family $(B_\epsilon)_{0<\epsilon\le 1}$ and is given by (3.4). In particular, for any $x \in M$ and any $0 < \epsilon \le 1$ Proposition 4.1 shows that

$$ |u_\epsilon(x) - B^{\mu_\epsilon, x_\epsilon}(x)| \le \sigma_\epsilon B_\epsilon^+(x). $$

As a consequence, if $(y_\epsilon)_{0<\epsilon\le 1}$ is any family of points in $M$ we have

$$ u_\epsilon(y_\epsilon) = B^{\mu_\epsilon, x_\epsilon}(y_\epsilon) + o(B_\epsilon^+(y_\epsilon)) $$

as $\epsilon \to 0$.

Proposition 4.1 is sharp if $\lambda(V) \ne 0$ which, by (2.4), corresponds to a bubble that decays to infinity in $\mathbb{R}^n$ as $|x|^{2-n}$. This is the case for the positive bubble $B_0^+$ given by (2.2) and for the sign-changing solutions of (1.2) constructed in DelPino-Musso-Pacard-Pistoia [7]. Proposition 4.1 for positive solutions has been known since the work of Li-Zhu [20] and Druet [10] and Druet-Hebey-Robert [12]] (see also [16]). In this section we prove Proposition 4.1 by following the approach in Ghoussoub-Mazumdar-Robert [14], that draws inspiration from the techniques in [12, 16]. Since we consider sign-changing solutions and bubbles we have to adapt the existing techniques that have been developed for positive solutions. In the exactly critical case $p_\epsilon \equiv 2^*$, Proposition 4.1 has already been proven in Premoselli [28] in the more general case where multiple sign-changing bubbles may appear.

4.2. **Proof of Proposition 4.1.** As a preliminary observation we claim that

$$ (4.6) \qquad c = \lim_{\epsilon \to 0} \mu_\epsilon^{-(\frac{n-2}{2})(2^*-p_\epsilon)} \in [1, +\infty). $$

We set $\tilde{U}_\epsilon(x) := \mu_\epsilon^{\frac{n-2}{2}} u_\epsilon(\exp_{x_\epsilon}(\mu_\epsilon x))$ for $x \in B(0, \frac{i_g(M)}{2\mu_\epsilon})$. It follows from (1.4) and (3.4) that $\tilde{U}_\epsilon \to cV$ in $H_{\text{loc}}^1(\mathbb{R}^n)$. With a change of variable, equation (1.3) writes

$$ \Delta_{g_\epsilon} \tilde{U}_\epsilon + \mu_\epsilon^2 h_\epsilon (\exp_{x_\epsilon}(\mu_\epsilon x)) \tilde{U}_\epsilon = \mu_\epsilon^{\frac{n-2}{2}(2^*-p_\epsilon)} |\tilde{U}_\epsilon|^{2^*-2} \tilde{U}_\epsilon \text{ in } B(0, \frac{i_g(M)}{2\mu_\epsilon}), $$

where we have let $g_\epsilon(x) = \exp_{x_\epsilon}^* g(\mu_\epsilon x)$. Passing to the limit $\epsilon \to 0$ yields

$$ \Delta_{\xi}(cV) = \left( \lim_{\epsilon \to 0} \mu_{\epsilon}^{\frac{n-2}{2}(2^*-p_{\epsilon})} \right) c^{2^*-1} |V|^{2^*-2V} \text{ weakly in } \mathbb{R}^n. $$

Since $\mu_\epsilon \to 0$ as $\epsilon \to 0$ and $p_\epsilon \le 2^*$ we have $\mu_\epsilon^{2^*-p_\epsilon} \le 1$ for $\epsilon$ small enough. Since $V \in D^{1,2}(\mathbb{R}^n)\setminus\{0\}$ it cannot be harmonic, which shows that $\lim_{\epsilon\to 0} \mu_\epsilon^{\frac{n-2}{2}(2^*-p_\epsilon)} \in (0, 1]$. Finally, $V$ solves (1.2), which proves the claim.
**Proposition 4.1.** Let $(M^n, g)$, $n \ge 3$ be a smooth, connected and closed Riemannian manifold, let $(h_\epsilon)_{0<\epsilon\le 1}$ be a family of $C^1$ functions in $M$ converging in $C^1(M)$ towards $h$ and let $(p_\epsilon)_{0<\epsilon\le 1}$ be a family of real numbers satisfying $2<p_\epsilon\le 2^*$ for all $0<\epsilon\le 1$ and such that $\lim_{\epsilon\to 0} p_\epsilon = 2^*$. Assume that $\Delta_g + h$ is coercive. Let $(u_\epsilon)_{0<\epsilon\le 1}$ be a family of solutions of (1.3) satisfying (1.4) for some $c>0$.

Then $c=1$ in (1.4) and there exists a family $(\sigma_\epsilon)_{0<\epsilon\le 1}$ of positive real numbers with $\lim_{\epsilon\to 0} \sigma_\epsilon = 0$ such that,

$$ (4.5) \qquad \left\| \frac{u_{\epsilon} - B^{\mu_{\epsilon}, x_{\epsilon}}}{B_{\epsilon}^{+}} \right\|_{L^{\infty}(M)} \le \sigma_{\epsilon} $$

for any $0 < \epsilon \le 1$, where $B^{\mu_\epsilon, x_\epsilon}$ is as in (4.2) and $B_\epsilon^+$ is as in (3.5).

The function $V \in \Sigma$ that appears in the definition of $B^{\mu_\epsilon, x_\epsilon}$ is the one that defines the family $(B_\epsilon)_{0<\epsilon\le 1}$ and is given by (3.4). In particular, for any $x \in M$ and any $0 < \epsilon \le 1$ Proposition 4.1 shows that

$$ |u_{\epsilon}(x) - B^{\mu_{\epsilon}, x_{\epsilon}}(x)| \le \sigma_{\epsilon} B_{\epsilon}^{+}(x). $$

As a consequence, if $(y_\epsilon)_{0<\epsilon\le 1}$ is any family of points in $M$ we have

$$ u_{\epsilon}(y_{\epsilon}) = B^{\mu_{\epsilon}, x_{\epsilon}}(y_{\epsilon}) + o(B_{\epsilon}^{+}(y_{\epsilon})) $$

as $\epsilon \to 0$.

Proposition 4.1 is sharp if $\lambda(V) \ne 0$ which, by (2.4), corresponds to a bubble that decays to infinity in $\mathbb{R}^n$ as $|x|^{2-n}$. This is the case for the positive bubble $B_0^+$ given by (2.2) and for the sign-changing solutions of (1.2) constructed in DelPino-Musso-Pacard-Pistoia [7]. Proposition 4.1 for positive solutions has been known since the work of Li-Zhu [20] and Druet [10] and Druet-Hebey-Robert [12]] (see also [16]). In this section we prove Proposition 4.1 by following the approach in Ghoussoub-Mazumdar-Robert [14], that draws inspiration from the techniques in [12, 16]. Since we consider sign-changing solutions and bubbles we have to adapt the existing techniques that have been developed for positive solutions. In the exactly critical case $p_\epsilon \equiv 2^*$, Proposition 4.1 has already been proven in Premoselli [28] in the more general case where multiple sign-changing bubbles may appear.

**4.2. Proof of Proposition 4.1.** As a preliminary observation we claim that

$$ (4.6) \qquad c = \lim_{\epsilon \to 0} \mu_{\epsilon}^{-\frac{(n-2)^2}{8}(2^*-p_{\epsilon})} \in [1, +\infty). $$

We set $\tilde{U}_\epsilon(x) := \mu_\epsilon^{\frac{n-2}{2}} u_\epsilon(\exp_{x_\epsilon}(\mu_\epsilon x))$ for $x \in B(0, \frac{i_g(M)}{2\mu_\epsilon})$. It follows from (1.4) and (3.4) that $\tilde{U}_\epsilon \to cV$ in $H_{\text{loc}}^1(\mathbb{R}^n)$. With a change of variable, equation (1.3) writes

$$ \Delta_{g_\epsilon} \tilde{U}_\epsilon + \mu_\epsilon^2 h_\epsilon(\exp_{x_\epsilon}(\mu_\epsilon x)) \tilde{U}_\epsilon = \mu_\epsilon^{\frac{n-2}{2}(2^*-p_\epsilon)} |\tilde{U}_\epsilon|^{2^*-2} \tilde{U}_\epsilon \text{ in } B(0, \frac{i_g(M)}{2\mu_\epsilon}), $$

where we have let $g_\epsilon(x) = \exp_{x_\epsilon}^* g(\mu_\epsilon x)$ and where $\exp_p^* g$ denotes the pull-back of the metric $g$ via the map $\exp_p$. Passing to the limit $\epsilon \to 0$ yields

$$ \Delta_{\xi}(cV) = \left( \lim_{\epsilon \to 0} \mu_{\epsilon}^{\frac{n-2}{2}(2^*-p_{\epsilon})} \right) c^{2^*-1} |V|^{2^*-2} V \quad \text{weakly in } \mathbb{R}^n. $$

Since $\mu_\epsilon \to 0$ as $\epsilon \to 0$ and $p_\epsilon \le 2^*$ we have $\mu_\epsilon^{2^*-p_\epsilon} \le 1$ for $\epsilon$ small enough. Since $V \in D^{1,2}(\mathbb{R}^n)\setminus\{0\}$ it cannot be harmonic, which shows that $\lim_{\epsilon\to 0} \mu_\epsilon^{\frac{n-2}{2}(2^*-p_\epsilon)} \in (0, 1]$. Finally, $V$ solves (1.2), which proves the claim.
As already remarked, we can assume that (4.3) holds. This implies that the family $(u_\varepsilon)_{0<\varepsilon\le 1}$ blows-up, that is

$$ \lim_{\varepsilon \to 0} \|u_\varepsilon\|_{L^\infty(M)} = +\infty $$

By (4.6) it is also easily seen that

$$ (4.7) \qquad \mu_{\varepsilon}^{\frac{n-2}{2} - \frac{2}{p_{\varepsilon}-2}} \to c \quad \text{as } \varepsilon \to 0. $$

Decomposition (4.3) with (4.2) shows in particular that $u_\varepsilon \to 0$ in $H^1(M\setminus B_g(x_\varepsilon, \delta))$ as $\varepsilon \to 0$, for any $\delta > 0$ fixed. The regularity theory of Trudinger [37], together with standard elliptic theory, then shows that

$$ (4.8) \qquad u_\varepsilon \to 0 \quad \text{in} \quad C_{\text{loc}}^2(M\setminus\{x_0\}) $$

as $\varepsilon \to 0$, where we have let $x_0 = \lim_{\varepsilon \to 0} x_\varepsilon$. We prove Proposition 4.1 in several steps.

**Step 1:** Let, for any $x \in B(0, \frac{i_g(M)}{2\mu_\varepsilon})$,

$$ U_\varepsilon(x) = \mu_\varepsilon^{\frac{2}{p_\varepsilon-2}} u_\varepsilon (\exp_{x_\varepsilon}(\mu_\varepsilon x)). $$

We claim that

$$ (4.9) \qquad U_\varepsilon \to V \quad \text{in} \quad C_{\text{loc}}^2(\mathbb{R}^n) $$

as $\varepsilon \to 0$.

*Proof of (4.9).* Let, for $x \in B(0, \frac{i_g(M)}{2\mu_\varepsilon})$, $g_\varepsilon(x) = \exp_{x_\varepsilon}^* g(\mu_\varepsilon x)$. By (1.3) $U_\varepsilon$ satisfies

$$ (4.10) \qquad \Delta_{g_\varepsilon} U_\varepsilon + \mu_\varepsilon^2 \tilde{h}_\varepsilon U_\varepsilon = |U_\varepsilon|^{p_\varepsilon-2} U_\varepsilon \quad \text{in} \quad B\left(0, \frac{i_g(M)}{2\mu_\varepsilon}\right), $$

where we have let $\tilde{h}_\varepsilon(x) = h_\varepsilon(\exp_{x_\varepsilon}(\mu_\varepsilon x))$. Equation (4.3) together with (4.7) shows that, for any compact set $K \subset \mathbb{R}^n$, we have

$$ \begin{align*} \int_K |U_\varepsilon - V|^{p_\varepsilon} dv_{g_\varepsilon} &= \mu_\varepsilon^{\frac{2p_\varepsilon}{p_\varepsilon-2}-n} \int_{\exp_{x_\varepsilon}(\mu_\varepsilon K)} \left| u_\varepsilon - \mu_\varepsilon^{\frac{n-2}{2} - \frac{2}{p_\varepsilon-2}} B^{\mu_\varepsilon, x_\varepsilon} \right|^{p_\varepsilon} dv_g \\ &= \mu_\varepsilon^{\frac{2p_\varepsilon}{p_\varepsilon-2}-n} \int_{\exp_{x_\varepsilon}(\mu_\varepsilon K)} \left| u_\varepsilon - (c+o(1)) B^{\mu_\varepsilon, x_\varepsilon} \right|^{p_\varepsilon} dv_g \\ &\lesssim \int_{\exp_{x_\varepsilon}(\mu_\varepsilon K)} \left| u_\varepsilon - (c+o(1)) B^{\mu_\varepsilon, x_\varepsilon} \right|^{p_\varepsilon} dv_g \\ &= o(1) \end{align*} $$

as $\varepsilon \to 0$, where the second line follows from (4.7) and in the third line we used that $\mu_\varepsilon^{\frac{2p_\varepsilon}{p_\varepsilon-2}-n} \lesssim 1$ since $\frac{2p_\varepsilon}{p_\varepsilon-2} - n \ge 0$. The local regularity theory of Trudinger [37] and standard elliptic theory using (4.10) then show that the convergence of $U_\varepsilon$ towards $V$ takes place in $C_{\text{loc}}^2(\mathbb{R}^n)$, which concludes the proof of (4.9). $\square$

By (4.3), (4.9) and Sobolev's inequality, we get

$$ (4.11) \qquad \lim_{R \to +\infty} \limsup_{\varepsilon \to 0} \|u_\varepsilon\|_{L^{2^*}(M \setminus B_g(x_\varepsilon, R\mu_\varepsilon))} = 0. $$

**Step 2:** Let $0 < \delta < i_g(M)$ be fixed. We claim that

$$ (4.12) \qquad \lim_{R \to +\infty} \limsup_{\varepsilon \to 0} \max_{B_g(x_\varepsilon, \delta) \setminus B_g(x_\varepsilon, R\mu_\varepsilon)} d_g(x_\varepsilon, x)^{\frac{2}{p_\varepsilon-2}} |u_\varepsilon(x)| = 0. $$
As already remarked, we can assume that (4.3) holds. This implies that the family $(u_\epsilon)_{0<\epsilon\le 1}$ blows-up, that is

$$ \lim_{\epsilon \to 0} \|u_\epsilon\|_{L^\infty(M)} = +\infty $$

By (4.6) it is also easily seen that

$$ (4.7) \qquad \mu_{\epsilon}^{\frac{n-2}{2} - \frac{2}{p_{\epsilon}-2}} \to c \quad \text{as } \epsilon \to 0. $$

Decomposition (4.3) with (4.2) shows in particular that $u_\epsilon \to 0$ in $H^1(M\setminus B_g(x_\epsilon, \delta))$ as $\epsilon \to 0$, for any $\delta > 0$ fixed. The regularity theory of Trudinger [37], together with standard elliptic theory, then shows that

$$ (4.8) \qquad u_{\epsilon} \to 0 \quad \text{in} \quad C_{\text{loc}}^{2}(M\setminus\{x_0\}) $$

as $\epsilon \to 0$, where we have let $x_0 = \lim_{\epsilon \to 0} x_\epsilon$. We prove Proposition 4.1 in several steps.

**Step 1:** Let, for any $x \in B(0, \frac{i_g(M)}{2\mu_\epsilon})$,

$$ U_{\epsilon}(x) = \mu_{\epsilon}^{\frac{2}{p_{\epsilon}-2}} u_{\epsilon} (\exp_{x_{\epsilon}}(\mu_{\epsilon}x)). $$

We claim that

$$ (4.9) \qquad U_{\epsilon} \to V \quad \text{in} \quad C_{\text{loc}}^{2}(\mathbb{R}^{n}) $$

as $\epsilon \to 0$.

*Proof of (4.9).* Let, for $x \in B(0, \frac{i_g(M)}{2\mu_\epsilon})$, $g_\epsilon(x) = \exp_{x_\epsilon}^* g(\mu_\epsilon x)$. By (1.3) $U_\epsilon$ satisfies

$$ (4.10) \qquad \Delta_{g_\epsilon} U_\epsilon + \mu_\epsilon^2 \tilde{h}_\epsilon U_\epsilon = |U_\epsilon|^{p_\epsilon-2} U_\epsilon \quad \text{in} \quad B\left(0, \frac{i_g(M)}{2\mu_\epsilon}\right), $$

where we have let $\tilde{h}_\epsilon(x) = h_\epsilon(\exp_{x_\epsilon}(\mu_\epsilon x))$. Equation (4.3) together with (4.7) shows that, for any compact set $K \subset \mathbb{R}^n$, we have

$$ \begin{align*} \int_K |U_\epsilon - V|^{p_\epsilon} dv_{g_\epsilon} &= \mu_\epsilon^{\frac{2p_\epsilon}{p_\epsilon-2}-n} \int_{\exp_{x_\epsilon}(\mu_\epsilon K)} \left| u_\epsilon - \mu_\epsilon^{\frac{n-2}{2} - \frac{2}{p_\epsilon-2}} B^{\mu_\epsilon, x_\epsilon} \right|^{p_\epsilon} dv_g \\ &= \mu_\epsilon^{\frac{2p_\epsilon}{p_\epsilon-2}-n} \int_{\exp_{x_\epsilon}(\mu_\epsilon K)} \left| u_\epsilon - (c+o(1)) B^{\mu_\epsilon, x_\epsilon} \right|^{p_\epsilon} dv_g \\ &\lesssim \int_{\exp_{x_\epsilon}(\mu_\epsilon K)} \left| u_\epsilon - (c+o(1)) B^{\mu_\epsilon, x_\epsilon} \right|^{p_\epsilon} dv_g \\ &= o(1) \end{align*} $$

as $\epsilon \to 0$, where the second line follows from (4.7) and in the third line we used that $\mu_\epsilon^{\frac{2p_\epsilon}{p_\epsilon-2}-n} \lesssim 1$ since $\frac{2p_\epsilon}{p_\epsilon-2} - n \ge 0$. The local regularity theory of Trudinger [37] and standard elliptic theory using (4.10) then show that the convergence of $U_\epsilon$ towards $V$ takes place in $C_{\text{loc}}^2(\mathbb{R}^n)$, which concludes the proof of (4.9). $\square$

By (4.3), (4.9) and Sobolev's inequality, we get

$$ (4.11) \qquad \lim_{R \to +\infty} \limsup_{\epsilon \to 0} \|u_\epsilon\|_{L^{2^*}(M \setminus B_g(x_\epsilon, R\mu_\epsilon))} = 0. $$

**Step 2:** Let $0 < \delta < i_g(M)$ be fixed. We claim that

$$ (4.12) \qquad \lim_{R \to +\infty} \limsup_{\epsilon \to 0} \max_{B_g(x_\epsilon, \delta) \setminus B_g(x_\epsilon, R\mu_\epsilon)} d_g(x_\epsilon, x)^{\frac{2}{p_\epsilon-2}} |u_\epsilon(x)| = 0. $$
*Proof.* To prove (4.12) we will prove more generally that

$$ (4.13) \quad \limsup_{\varepsilon \to 0} \max_{x \in M} (\mu_{\varepsilon} + d_g(x_{\varepsilon}, x))^{\frac{2}{p_{\varepsilon}-2}} |u_{\varepsilon} - cB^{\mu_{\varepsilon}, x_{\varepsilon}}|(x) = 0. $$

It is easily seen using (4.2) that (4.12) follows from (4.13). To prove (4.13) we proceed by contradiction and assume that, up to passing to a subsequence $(\varepsilon_k)_{k \ge 0}$ with $\lim_{k \to +\infty} \varepsilon_k = 0$, there exists $\eta_0 > 0$ and a sequence $(y_k)_{k \ge 0}$ of points of $M$ such that

$$ (4.14) \quad \begin{aligned} & (\mu_{\varepsilon_k} + d_g(x_{\varepsilon_k}, y_k))^{\frac{2}{p_{\varepsilon_k}-2}} |u_{\varepsilon_k} - cB^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}|(y_k) \\ & = \max_{x \in M} (\mu_{\varepsilon_k} + d_g(x_{\varepsilon_k}, x))^{\frac{2}{p_{\varepsilon_k}-2}} |u_{\varepsilon_k} - cB^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}|(x) \ge \eta_0. \end{aligned} $$

By (4.7) and (4.9) we get that

$$ (4.15) \quad \frac{d_g(x_{\varepsilon_k}, y_k)}{\mu_{\varepsilon_k}} \to +\infty $$

as $k \to +\infty$. By (4.2) and since $\frac{2}{p_{\varepsilon_k}-2} - \frac{n-2}{2} \ge 0$ we thus have

$$ (4.16) \quad \begin{aligned} & (\mu_{\varepsilon_k} + d_g(x_{\varepsilon_k}, y_k))^{\frac{2}{p_{\varepsilon_k}-2}} |B^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}|(y_k) \\ & \lesssim \left( \frac{\mu_{\varepsilon_k}}{\mu_{\varepsilon} + d_g(x_{\varepsilon_k}, y_k)} \right)^{\frac{n-2}{2}} (\mu_{\varepsilon} + d_g(x_{\varepsilon_k}, y_k))^{\frac{2}{p_{\varepsilon_k}-2} - \frac{n-2}{2}} \\ & = o(1) \end{aligned} $$

as $k \to +\infty$. Using (4.2) and (4.8), (4.14) also shows that $d_g(x_{\varepsilon_k}, y_k) \to 0$, and hence, with (4.16), that $|u_{\varepsilon_k}(y_k)| \to +\infty$ as $k \to +\infty$. We let $\nu_k = |u_{\varepsilon_k}(y_k)|^{-\frac{p_{\varepsilon_k}-2}{2}}$, so that $\nu_k \to 0$, and for $x \in B(0, \frac{i_g(M)}{2\nu_k})$ we let

$$ w_k(x) = \nu_k^{\frac{2}{p_{\varepsilon_k}-2}} u_{\varepsilon_k} \left( \exp_{y_k} (\nu_k x) \right). $$

By (1.3) $w_k$ satisfies

$$ (4.17) \quad \Delta_{g_k} w_k + \nu_k^2 \tilde{h}_k w_k = |w_k|^{p_{\varepsilon_k}-2} w_k \quad \text{in } B\left(0, \frac{i_g(M)}{2\nu_k}\right), $$

where we have let $\tilde{h}_k(x) = h_{\varepsilon_k}(\exp_{y_k}(\nu_k x))$ and $g_k = \exp_{y_k}^* g(\nu_k x)$. Independently, by (4.14), (4.15) and (4.16) we see that

$$ (4.18) \quad d_g(x_{\varepsilon_k}, y_k) \gtrsim \nu_k. $$

With the latter, (4.14) shows that there exists $r_0 > 0$ small enough such that

$$ \|w_k\|_{L^\infty(B(0,2r_0))} \le 2 $$

for any $k$ large enough. With (4.17) and standard elliptic theory, and since $|w_k(0)| = 1$, $w_k$ converges in $C^2(\bar{B}(0, r_0))$ as $k \to +\infty$, up to a subsequence, to a non-zero solution $w_0$ of

$$ \Delta_{\xi} w_0 = |w_0|^{2^*-2} w_0 \quad \text{in } B(0, r_0), $$

that thus satisfies

$$ (4.19) \quad \int_{B(0,r_0)} |w_0|^{2^*} dx > 0. $$
*Proof.* To prove (4.12) we will prove more generally that

$$ (4.13) \quad \limsup_{\varepsilon \to 0} \max_{x \in M} (\mu_{\varepsilon} + d_g(x_{\varepsilon}, x))^{\frac{2}{p_{\varepsilon}-2}} |u_{\varepsilon} - cB^{\mu_{\varepsilon}, x_{\varepsilon}}|(x) = 0. $$

It is easily seen using (4.2) that (4.12) follows from (4.13). To prove (4.13) we proceed by contradiction and assume that, up to passing to a subsequence $(\varepsilon_k)_{k \ge 0}$ with $\lim_{k \to +\infty} \varepsilon_k = 0$, there exists $\eta_0 > 0$ and a sequence $(y_k)_{k \ge 0}$ of points of $M$ such that

$$ (4.14) \quad \begin{aligned} & (\mu_{\varepsilon_k} + d_g(x_{\varepsilon_k}, y_k))^{\frac{2}{p_{\varepsilon_k}-2}} |u_{\varepsilon_k} - cB^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}|(y_k) \\ & = \max_{x \in M} (\mu_{\varepsilon_k} + d_g(x_{\varepsilon_k}, x))^{\frac{2}{p_{\varepsilon_k}-2}} |u_{\varepsilon_k} - cB^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}|(x) \ge \eta_0. \end{aligned} $$

By (4.7) and (4.9) we get that

$$ (4.15) \quad \frac{d_g(x_{\varepsilon_k}, y_k)}{\mu_{\varepsilon_k}} \to +\infty $$

as $k \to +\infty$. By (4.2) and since $\frac{2}{p_{\varepsilon_k}-2} - \frac{n-2}{2} \ge 0$ we thus have

$$ (4.16) \quad \begin{aligned} & (\mu_{\varepsilon_k} + d_g(x_{\varepsilon_k}, y_k))^{\frac{2}{p_{\varepsilon_k}-2}} |B^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}|(y_k) \\ & \lesssim \left( \frac{\mu_{\varepsilon_k}}{\mu_{\varepsilon} + d_g(x_{\varepsilon_k}, y_k)} \right)^{\frac{n-2}{2}} (\mu_{\varepsilon} + d_g(x_{\varepsilon_k}, y_k))^{\frac{2}{p_{\varepsilon_k}-2} - \frac{n-2}{2}} \\ & = o(1) \end{aligned} $$

as $k \to +\infty$. Using (4.2) and (4.8), (4.14) also shows that $d_g(x_{\varepsilon_k}, y_k) \to 0$, and hence, with (4.16), that $|u_{\varepsilon_k}(y_k)| \to +\infty$ as $k \to +\infty$. We let $\nu_k = |u_{\varepsilon_k}(y_k)|^{-\frac{p_{\varepsilon_k}-2}{2}}$, so that $\nu_k \to 0$, and for $x \in B(0, \frac{i_g(M)}{2\nu_k})$ we let

$$ w_k(x) = \nu_k^{\frac{2}{p_{\varepsilon_k}-2}} u_{\varepsilon_k} \left( \exp_{y_k} (\nu_k x) \right). $$

By (1.3) $w_k$ satisfies

$$ (4.17) \quad \Delta_{g_k} w_k + \nu_k^2 \tilde{h}_k w_k = |w_k|^{p_{\varepsilon_k}-2} w_k \quad \text{in } B\left(0, \frac{i_g(M)}{2\nu_k}\right), $$

where we have let $\tilde{h}_k(x) = h_{\varepsilon_k}(\exp_{y_k}(\nu_k x))$ and $g_k = \exp_{y_k}^* g(\nu_k x)$. Independently, by (4.14), (4.15) and (4.16) we see that

$$ (4.18) \quad d_g(x_{\varepsilon_k}, y_k) \gtrsim \nu_k. $$

With the latter, (4.14) shows that there exists $r_0 > 0$ small enough such that

$$ \|w_k\|_{L^\infty(B(0,2r_0))} \le 2 $$

for any $k$ large enough. With (4.17) and standard elliptic theory, and since $|w_k(0)| = 1$, $w_k$ converges in $C^2(\bar{B}(0, r_0))$ as $k \to +\infty$, up to a subsequence, to a non-zero solution $w_0$ of

$$ \Delta_{\xi} w_0 = |w_0|^{2^*-2} w_0 \quad \text{in } B(0, r_0), $$

that thus satisfies

$$ (4.19) \quad \int_{B(0,r_0)} |w_0|^{2^*} dx > 0. $$
Independently, (4.15) and (4.18) show that, up to assuming that $r_0$ is small enough, we have

$$ (4.20) \qquad \frac{1}{\mu_{\varepsilon_k}} d_g(x_{\varepsilon_k}, B_g(y_k, r_0 \nu_k)) \to +\infty $$

as $k \to +\infty$. We let $q_k = \frac{n}{2}p_{\varepsilon_k} - n$. As is easily seen, it satisfies $q_k \le 2^*$ (since $p_{\varepsilon_k} \le 2^*$), $q_k \to 2^*$ as $k \to +\infty$ and

$$ \frac{2q_k}{p_{\varepsilon_k} - 2} = n. $$

We then have, by the $C^2$ convergence of $w_k$ towards $w_0$ and by Hölder's inequality,

$$ \begin{align*} \int_{B(0,r_0)} |w_0|^{2^*} dx &= \limsup_{k \to +\infty} \int_{B(0,r_0)} |w_k|^{q_k} dv_{g_k} \\ &= \limsup_{k \to +\infty} \int_{B_g(y_k, r_0 \nu_k)} |u_{\varepsilon_k}|^{q_k} dv_g \\ &\lesssim \limsup_{k \to +\infty} \left( \int_{B_g(y_k, r_0 \nu_k)} |u_{\varepsilon_k}|^{2^*} dv_g \right)^{\frac{q_k}{2^*}} \\ &= 0, \end{align*} $$

where the last equality follows from (4.11) and (4.20). This is a contradiction with (4.19). Hence (4.14) cannot occur, which proves (4.13). $\square$

For $0 < \delta < i_g(M)$ we let

$$ (4.21) \qquad \eta_{\varepsilon}(\delta) = \|u_{\varepsilon}\|_{L^{\infty}(M \setminus B_g(x_{\varepsilon}, \delta))}. $$

By (4.8), $\eta_{\varepsilon}(\delta) \to 0$ as $\varepsilon \to 0$ for any $0 < \delta < i_g(M)$ fixed.

**Step 3:** We claim that for any $0 < \alpha < \frac{1}{2}$ there exists $0 < \delta_{\alpha} < i_g(M)$ and $C_{\alpha} > 0$ such that, for any $x \in M \setminus B_g(x_{\varepsilon}, \mu_{\varepsilon})$ and any $0 < \varepsilon \le 1$

$$ (4.22) \qquad |u_{\varepsilon}(x)| \le C_{\alpha} \left( \mu_{\varepsilon}^{\frac{2}{p_{\varepsilon}-2}(1-2\alpha)} d_g(x_{\varepsilon}, x)^{-\frac{4(1-\alpha)}{p_{\varepsilon}-2}} + \eta_{\varepsilon}(\delta_{\alpha}) d_g(x_{\varepsilon}, x)^{-\frac{4\alpha}{p_{\varepsilon}-2}} \right) $$

holds.

*Proof of (4.22).* Let $0 < \alpha < \frac{1}{2}$. For $x \in M \setminus B_g(x_{\varepsilon}, \mu_{\varepsilon})$ we let

$$ \Phi_{\varepsilon}^{\alpha}(x) = \mu_{\varepsilon}^{\frac{2}{p_{\varepsilon}-2}(1-2\alpha)} G_{h_{\varepsilon}}(x_{\varepsilon}, x)^{\frac{4(1-\alpha)}{(n-2)(p_{\varepsilon}-2)}} + \eta_{\varepsilon}(\delta) G_{h_{\varepsilon}}(x_{\varepsilon}, x)^{\frac{4\alpha}{(n-2)(p_{\varepsilon}-2)}}, $$

where $G_{h_{\varepsilon}}$ is the Green's function of $\Delta_g + h_{\varepsilon}$. It follows from properties of the Green's function that there exists $\delta'_{\alpha} > 0$ such that

$$ d_g(x, x_{\varepsilon})^2 \left( \frac{|\nabla_x G_{h_{\varepsilon}}(x_{\varepsilon}, x)|^2}{G_{h_{\varepsilon}}(x_{\varepsilon}, x)^2} - h_{\varepsilon}(x) \right) \ge \frac{(n-2)^2}{4} \text{ for } d_g(x, x_{\varepsilon}) < \delta'_{\alpha}. $$

Let $\delta_{\alpha} \in (0, \delta'_{\alpha})$ be such that

$$ 2 \|h\|_{L^{\infty}(M)} \delta_{\alpha}^{2} < (n-2)^{2} \alpha (1-\alpha) $$

and let $0 < \delta \le \delta_{\alpha}$ be fixed. We claim that there exists $C > 0$ such that, for any $0 < \varepsilon \le 1$,

$$ (4.23) \qquad \left\| \frac{u_{\varepsilon}}{\Phi_{\varepsilon}^{\alpha}} \right\|_{L^{\infty}(M \setminus B_g(x_{\varepsilon}, \mu_{\varepsilon}))} \le C $$
Independently, (4.15) and (4.18) show that, up to assuming that $r_0$ is small enough, we have

$$ (4.20) \qquad \frac{1}{\mu_{\varepsilon_k}} d_g(x_{\varepsilon_k}, B_g(y_k, r_0 \nu_k)) \to +\infty $$

as $k \to +\infty$. We let $q_k = \frac{n}{2}p_{\varepsilon_k} - n$. As is easily seen, it satisfies $q_k \le 2^*$ (since $p_{\varepsilon_k} \le 2^*$), $q_k \to 2^*$ as $k \to +\infty$ and

$$ \frac{2q_k}{p_{\varepsilon_k} - 2} = n. $$

We then have, by the $C^2$ convergence of $w_k$ towards $w_0$ and by Hölder's inequality,

$$ \begin{align*} \int_{B(0,r_0)} |w_0|^{2^*} dx &= \limsup_{k \to +\infty} \int_{B(0,r_0)} |w_k|^{q_k} dv_{g_k} \\ &= \limsup_{k \to +\infty} \int_{B_g(y_k, r_0 \nu_k)} |u_{\varepsilon_k}|^{q_k} dv_g \\ &\lesssim \limsup_{k \to +\infty} \left( \int_{B_g(y_k, r_0 \nu_k)} |u_{\varepsilon_k}|^{2^*} dv_g \right)^{\frac{q_k}{2^*}} \\ &= 0, \end{align*} $$

where the last equality follows from (4.11) and (4.20). This is a contradiction with (4.19). Hence (4.14) cannot occur, which proves (4.13). $\square$

For $0 < \delta < i_g(M)$ we let

$$ (4.21) \qquad \eta_{\varepsilon}(\delta) = \|u_{\varepsilon}\|_{L^{\infty}(M \setminus B_g(x_{\varepsilon}, \delta))}. $$

By (4.8), $\eta_{\varepsilon}(\delta) \to 0$ as $\varepsilon \to 0$ for any $0 < \delta < i_g(M)$ fixed.

**Step 3:** We claim that for any $0 < \alpha < \frac{1}{2}$ there exists $0 < \delta_{\alpha} < i_g(M)$ and $C_{\alpha} > 0$ such that, for any $x \in M \setminus B_g(x_{\varepsilon}, \mu_{\varepsilon})$ and any $0 < \varepsilon \le 1$

$$ (4.22) \qquad |u_{\varepsilon}(x)| \le C_{\alpha} \left( \mu_{\varepsilon}^{\frac{2}{p_{\varepsilon}-2}(1-2\alpha)} d_g(x_{\varepsilon}, x)^{-\frac{4(1-\alpha)}{p_{\varepsilon}-2}} + \eta_{\varepsilon}(\delta_{\alpha}) d_g(x_{\varepsilon}, x)^{-\frac{4\alpha}{p_{\varepsilon}-2}} \right) $$

holds.

*Proof of (4.22).* Let $0 < \alpha < \frac{1}{2}$. For $x \in M \setminus B_g(x_{\varepsilon}, \mu_{\varepsilon})$ we let

$$ \Phi_{\varepsilon}^{\alpha}(x) = \mu_{\varepsilon}^{\frac{2}{p_{\varepsilon}-2}(1-2\alpha)} G_{h_{\varepsilon}}(x_{\varepsilon}, x)^{\frac{4(1-\alpha)}{(n-2)(p_{\varepsilon}-2)}} + \eta_{\varepsilon}(\delta) G_{h_{\varepsilon}}(x_{\varepsilon}, x)^{\frac{4\alpha}{(n-2)(p_{\varepsilon}-2)}}, $$

where $G_{h_{\varepsilon}}$ is the Green's function of $\Delta_g + h_{\varepsilon}$. The existence follows from the uniform coercivity of $\Delta_g + h$ and the $C^1$-convergence of $(h_{\varepsilon})_{\varepsilon}$ to $h$. It follows from properties of the Green's function that there exists $\delta'_{\alpha} > 0$ such that

$$ d_g(x, x_{\varepsilon})^2 \left( \frac{|\nabla_x G_{h_{\varepsilon}}(x_{\varepsilon}, x)|^2}{G_{h_{\varepsilon}}(x_{\varepsilon}, x)^2} - h_{\varepsilon}(x) \right) \ge \frac{(n-2)^2}{4} \quad \text{for } d_g(x, x_{\varepsilon}) < \delta'_{\alpha}. $$

Let $\delta_{\alpha} \in (0, \delta'_{\alpha})$ be such that

$$ 2\|h\|_{L^{\infty}(M)} \delta_{\alpha}^2 < (n-2)^2 \alpha(1-\alpha) $$

and let $0 < \delta \le \delta_{\alpha}$ be fixed. We claim that there exists $C > 0$ such that, for any $0 < \varepsilon \le 1$,

$$ (4.23) \qquad \left\| \frac{u_{\varepsilon}}{\Phi_{\varepsilon}^{\alpha}} \right\|_{L^{\infty}(M \setminus B_g(x_{\varepsilon}, \mu_{\varepsilon}))} \le C $$
holds. By (4.1) it is easily seen that (4.22) follows from (4.23). We thus prove (4.23). We let, for any $0 < \varepsilon \le 1$, $y_\varepsilon \in M \setminus B_g(x_\varepsilon, \mu_\varepsilon)$ be such that

$$ \left| \frac{u_\varepsilon(y_\varepsilon)}{\Phi_\varepsilon^\alpha(y_\varepsilon)} \right| = \max_{x \in M \setminus B_g(x_\varepsilon, \mu_\varepsilon)} \left| \frac{u_\varepsilon(x)}{\Phi_\varepsilon^\alpha(x)} \right|. $$

We proceed by contradiction and assume that, up to passing to a subsequence $(\varepsilon_k)_{k \ge 0}$ with $\lim_{k \to +\infty} \varepsilon_k = 0$, and up to replacing $u_{\varepsilon_k}$ by $-u_{\varepsilon_k}$ (which still solves (1.3)), we have

$$ (4.24) \qquad \frac{u_{\varepsilon_k}(y_{\varepsilon_k})}{\Phi_{\varepsilon_k}^{\alpha}(y_{\varepsilon_k})} \to +\infty $$

as $k \to +\infty$. We first assume that, up to a subsequence,

$$ (4.25) \qquad \mu_{\varepsilon_k} = o(d_g(y_{\varepsilon_k}, x_{\varepsilon_k})) \quad \text{and} \quad d_g(x_{\varepsilon_k}, y_{\varepsilon_k}) \le \delta_\alpha $$

as $k \to +\infty$. By definition of $y_{\varepsilon_k}$ we have

$$ (4.26) \qquad \frac{\Delta_g u_{\varepsilon_k}(y_{\varepsilon_k})}{u_{\varepsilon_k}(y_{\varepsilon_k})} \ge \frac{\Delta_g \Phi_{\varepsilon_k}^{\alpha}(y_{\varepsilon_k})}{\Phi_{\varepsilon_k}^{\alpha}(y_{\varepsilon_k})}. $$

On the one side, straightforward computations show that

$$ d_g(x_{\varepsilon_k}, y_{\varepsilon_k})^2 \frac{\Delta_g \Phi_{\varepsilon_k}^{\alpha}(y_{\varepsilon_k})}{\Phi_{\varepsilon_k}^{\alpha}(y_{\varepsilon_k})} \ge \frac{(n-2)^2}{2} \alpha(1-\alpha) > 0 \text{ for } d_g(x, x_\varepsilon) < \delta_\alpha $$

as $k \to +\infty$. On the other side, (1.3) and (4.25) show that

$$ d_g(x_{\varepsilon_k}, y_{\varepsilon_k})^2 \frac{\Delta_g u_{\varepsilon_k}(y_{\varepsilon_k})}{u_{\varepsilon_k}(y_{\varepsilon_k})} \le d_g(x_{\varepsilon_k}, y_{\varepsilon_k})^2 |u_{\varepsilon_k}(y_{\varepsilon_k})|^{p_{\varepsilon_k}-2} + 2 \|h\|_{L^\infty(M)} \delta_\alpha^2 $$

as $k \to +\infty$. With (4.12), (4.26) and the choice of $\delta_\alpha$ this yields a contradiction for $k$ large enough, and shows that (4.25) cannot happen. Thus we either have, up to a subsequence, $\frac{d_g(x_{\varepsilon_k}, y_{\varepsilon_k})}{\mu_{\varepsilon_k}} \to D \in [1, +\infty)$ or $d_g(x_{\varepsilon_k}, y_{\varepsilon_k}) \ge \delta_\alpha$ as $k \to +\infty$. Using (4.9) in the first case and (4.21) in the second case we obtain a contradiction with (4.24). This proves (4.23) and thus (4.22). $\square$

**Step 4:** We claim that

$$ (4.27) \qquad |u_\varepsilon(x)| \lesssim B_\varepsilon^+(x) \quad \text{for any } x \in M. $$

where $B_\varepsilon^+$ is defined in (3.5).

*Proof of (4.27).* We first prove a weaker statement. Let $0 < \delta < i_g(M)$ be fixed. We claim that there exists $C_\delta > 0$ such that, for any $0 < \varepsilon \le 1$ and $x \in M$,

$$ (4.28) \qquad |u_\varepsilon(x)| \le C_\delta (B_\varepsilon^+(x) + \eta_\varepsilon(\delta)) $$

holds, where $B_\varepsilon^+$ is as in (3.5). Let $(y_\varepsilon)_{0<\varepsilon\le 1}$ be any family of points in $M$. First, if $d_g(y_\varepsilon, x_\varepsilon) \ge \delta$, (4.28) follows from (4.21). We can thus assume that $y_\varepsilon \in B_g(x_\varepsilon, \delta)$ for all $0 < \varepsilon \le 1$. A representation formula for $u_\varepsilon$ at $y_\varepsilon$ gives, with (1.3),

$$ (4.29) \qquad \begin{aligned} u_\varepsilon(y_\varepsilon) = & \int_{B_g(x_\varepsilon, \mu_\varepsilon)} G_{h_\varepsilon}(y_\varepsilon, y) |u_\varepsilon(y)|^{p_\varepsilon-2} u_\varepsilon(y) dy \\ & + \int_{M \setminus B_g(x_\varepsilon, \mu_\varepsilon)} G_{h_\varepsilon}(y_\varepsilon, y) |u_\varepsilon(y)|^{p_\varepsilon-2} u_\varepsilon(y) dy. \end{aligned} $$
holds. By (4.1) it is easily seen that (4.22) follows from (4.23). We thus prove (4.23). We let, for any $0 < \varepsilon \le 1$, $y_\varepsilon \in M \setminus B_g(x_\varepsilon, \mu_\varepsilon)$ be such that

$$ \left| \frac{u_\varepsilon(y_\varepsilon)}{\Phi_\varepsilon^\alpha(y_\varepsilon)} \right| = \max_{x \in M \setminus B_g(x_\varepsilon, \mu_\varepsilon)} \left| \frac{u_\varepsilon(x)}{\Phi_\varepsilon^\alpha(x)} \right|. $$

We proceed by contradiction and assume that, up to passing to a subsequence $(\varepsilon_k)_{k \ge 0}$ with $\lim_{k \to +\infty} \varepsilon_k = 0$, and up to replacing $u_{\varepsilon_k}$ by $-u_{\varepsilon_k}$ (which still solves (1.3)), we have

$$ (4.24) \qquad \frac{u_{\varepsilon_k}(y_{\varepsilon_k})}{\Phi_{\varepsilon_k}^{\alpha}(y_{\varepsilon_k})} \to +\infty $$

as $k \to +\infty$. We first assume that, up to a subsequence,

$$ (4.25) \qquad \mu_{\varepsilon_k} = o(d_g(y_{\varepsilon_k}, x_{\varepsilon_k})) \quad \text{and} \quad d_g(x_{\varepsilon_k}, y_{\varepsilon_k}) \le \delta_\alpha $$

as $k \to +\infty$. By definition of $y_{\varepsilon_k}$ we have

$$ (4.26) \qquad \frac{\Delta_g u_{\varepsilon_k}(y_{\varepsilon_k})}{u_{\varepsilon_k}(y_{\varepsilon_k})} \ge \frac{\Delta_g \Phi_{\varepsilon_k}^{\alpha}(y_{\varepsilon_k})}{\Phi_{\varepsilon_k}^{\alpha}(y_{\varepsilon_k})}. $$

On the one side, straightforward computations show that

$$ d_g(x_{\varepsilon_k}, y_{\varepsilon_k})^2 \frac{\Delta_g \Phi_{\varepsilon_k}^{\alpha}(y_{\varepsilon_k})}{\Phi_{\varepsilon_k}^{\alpha}(y_{\varepsilon_k})} \ge \frac{(n-2)^2}{2} \alpha(1-\alpha) > 0 \text{ for } d_g(x, x_\varepsilon) < \delta_\alpha $$

as $k \to +\infty$. On the other side, (1.3) and (4.25) show that

$$ d_g(x_{\varepsilon_k}, y_{\varepsilon_k})^2 \frac{\Delta_g u_{\varepsilon_k}(y_{\varepsilon_k})}{u_{\varepsilon_k}(y_{\varepsilon_k})} \le d_g(x_{\varepsilon_k}, y_{\varepsilon_k})^2 |u_{\varepsilon_k}(y_{\varepsilon_k})|^{p_{\varepsilon_k}-2} + 2 \|h\|_{L^\infty(M)} \delta_\alpha^2 $$

as $k \to +\infty$. With (4.12), (4.26) and the choice of $\delta_\alpha$ this yields a contradiction for $k$ large enough, and shows that (4.25) cannot happen. Thus we either have, up to a subsequence, $\frac{d_g(x_{\varepsilon_k}, y_{\varepsilon_k})}{\mu_{\varepsilon_k}} \to D \in [1, +\infty)$ or $d_g(x_{\varepsilon_k}, y_{\varepsilon_k}) \ge \delta_\alpha$ as $k \to +\infty$. Using (4.9) in the first case and (4.21) in the second case we obtain a contradiction with (4.24). This proves (4.23) and thus (4.22). □

**Step 4:** We claim that

$$ (4.27) \qquad |u_\varepsilon(x)| \lesssim B_\varepsilon^+(x) \quad \text{for any } x \in M. $$

where $B_\varepsilon^+$ is defined in (3.5).

*Proof of (4.27).* We first prove a weaker statement. Let $0 < \delta < i_g(M)$ be fixed. We claim that there exists $C_\delta > 0$ such that, for any $0 < \varepsilon \le 1$ and $x \in M$,

$$ (4.28) \qquad |u_\varepsilon(x)| \le C_\delta (B_\varepsilon^+(x) + \eta_\varepsilon(\delta)) $$

holds, where $B_\varepsilon^+$ is as in (3.5). Let $(y_\varepsilon)_{0<\varepsilon\le 1}$ be any family of points in $M$. First, if $d_g(y_\varepsilon, x_\varepsilon) \ge \delta$, (4.28) follows from (4.21). We can thus assume that $y_\varepsilon \in B_g(x_\varepsilon, \delta)$ for all $0 < \varepsilon \le 1$. A representation formula for $u_\varepsilon$ at $y_\varepsilon$ gives, with (1.3),

$$ (4.29) \qquad \begin{aligned} u_\varepsilon(y_\varepsilon) = & \int_{B_g(x_\varepsilon, \mu_\varepsilon)} G_{h_\varepsilon}(y_\varepsilon, y) |u_\varepsilon(y)|^{p_\varepsilon-2} u_\varepsilon(y) dy \\ & + \int_{M \setminus B_g(x_\varepsilon, \mu_\varepsilon)} G_{h_\varepsilon}(y_\varepsilon, y) |u_\varepsilon(y)|^{p_\varepsilon-2} u_\varepsilon(y) dy. \end{aligned} $$
Straightforward computations with (4.7), (4.9) and (4.1) show that

$$ (4.30) \qquad \left| \int_{B_g(x_\varepsilon, \mu_\varepsilon)} G_{h_\varepsilon}(y_\varepsilon, y) |u_\varepsilon(y)|^{p_\varepsilon-2} u_\varepsilon(y) dy \right| \lesssim B_\varepsilon^+(y_\varepsilon) $$

(see e.g. [16, Proposition 6.1]). Straightforward computations using (4.7), (4.22) and (4.1) show that

$$ (4.31) \qquad \left| \int_{M \setminus B_g(x_\varepsilon, \mu_\varepsilon)} G_{h_\varepsilon}(y_\varepsilon, y) |u_\varepsilon(y)|^{p_\varepsilon-2} u_\varepsilon(y) dy \right| \lesssim B_\varepsilon^+(y_\varepsilon) + \eta_\varepsilon(\delta)^{p_\varepsilon-1}. $$

Plugging (4.30) and (4.31) in (4.29) and since $\eta_\varepsilon(\delta) \to 0$ as $\varepsilon \to 0$ proves (4.28). We now claim that, for any $0 < \delta < i_g(M)$ fixed,

$$ (4.32) \qquad \eta_\varepsilon(\delta) \lesssim \mu_\varepsilon^{\frac{n-2}{2}} $$

holds for $\varepsilon$ small enough. With (3.5), (4.28) and (4.32) this will conclude the proof of (4.27). We prove (4.32) by contradiction: we fix $0 < \delta < i_g(M)$ and, up to passing to a subsequence $(\varepsilon_k)_{k \ge 0}$ with $\lim_{k \to +\infty} \varepsilon_k = 0$ we assume that

$$ (4.33) \qquad \frac{\eta_{\varepsilon_k}(\delta)}{\mu_{\varepsilon_k}^{\frac{n-2}{2}}} \to +\infty $$

as $k \to +\infty$. We let, for $k \ge 0$ and for any $x \in M$,

$$ \hat{u}_k(x) = \frac{u_{\varepsilon_k}(x)}{\eta_{\varepsilon_k}(\delta)}, $$

and let $y_k \in M \setminus B_g(x_{\varepsilon_k}, \delta)$ be such that $|\hat{u}_k(y_k)| = 1$. By (3.5), (4.28) and (4.33) the family $(\hat{u}_k)_{k \ge 0}$ is uniformly bounded in compact subsets of $M \setminus \{x_0\}$, where $x_0 = \lim_{\varepsilon \to 0} x_\varepsilon$. By (1.3), standard elliptic theory and since $\eta_{\varepsilon_k}(\delta) \to 0$ as $k \to +\infty$, $\hat{u}_k$ thus converges in $C_{\text{loc}}^2(M \setminus \{x_0\})$, up to a subsequence, to some function $\hat{u}_0 \in C^2(M \setminus \{x_0\})$ that satisfies $\Delta_g \hat{u}_0 + h \hat{u}_0 = 0$ in $M \setminus \{x_0\}$ and $|\hat{u}_0(y_0)| = 1$, where $y_0 = \lim_{k \to +\infty} y_k \in M \setminus B_g(x_0, \delta)$. Passing (4.28) to the limit pointwise shows in addition that we have

$$ |\hat{u}_0(x)| \le C_\delta \quad \text{for all } x \in M \setminus \{x_0\}. $$

Classical arguments thus show that the singularity of $\hat{u}_0$ at $x_0$ is removable and that $\hat{u}_0$ satisfies $\Delta_g \hat{u}_0 + h \hat{u}_0 = 0$ in $M$ in a strong sense. This implies $\hat{u}_0 \equiv 0$ since $\Delta_g + h$ is coercive, which is a contradiction with $|\hat{u}_0(y_0)| = 1$. Thus (4.32) holds true and (4.27) is proven. $\square$

A consequence of (4.27), (1.3) and standard elliptic theory is the following estimate on the first and second derivatives of $u_\varepsilon$: for $0 \le k \le 2$,

$$ (4.34) \qquad |\nabla^k u_\varepsilon(x)| \lesssim \frac{\mu_\varepsilon^{\frac{n-2}{2}}}{(\mu_\varepsilon + d_g(x_\varepsilon, x))^{n-2+k}} \quad \text{for all } x \in M \text{ and } 0 < \varepsilon \le 1. $$

Estimate (4.34) follows from a scaling argument and from standard elliptic theory.

**Step 5:** we now claim that

$$ (4.35) \qquad 2^* - p_\varepsilon = \begin{cases} O(\mu_\varepsilon) & \text{if } n = 3 \\ O\left(\mu_\varepsilon^2 \ln \frac{1}{\mu_\varepsilon}\right) & \text{if } n = 4. \\ O(\mu_\varepsilon^2) & \text{if } n \ge 5 \end{cases}. $$
Straightforward computations with (4.7), (4.9) and (4.1) show that

$$ (4.30) \qquad \left| \int_{B_g(x_\varepsilon, \mu_\varepsilon)} G_{h_\varepsilon}(y_\varepsilon, y) |u_\varepsilon(y)|^{p_\varepsilon-2} u_\varepsilon(y) dy \right| \lesssim B_\varepsilon^+(y_\varepsilon) $$

(see e.g. [16, Proposition 6.1]). Straightforward computations using (4.7), (4.22) and (4.1) show that

$$ (4.31) \qquad \left| \int_{M \setminus B_g(x_\varepsilon, \mu_\varepsilon)} G_{h_\varepsilon}(y_\varepsilon, y) |u_\varepsilon(y)|^{p_\varepsilon-2} u_\varepsilon(y) dy \right| \lesssim B_\varepsilon^+(y_\varepsilon) + \eta_\varepsilon(\delta)^{p_\varepsilon-1}. $$

Plugging (4.30) and (4.31) in (4.29) and since $\eta_\varepsilon(\delta) \to 0$ as $\varepsilon \to 0$ proves (4.28). We now claim that, for any $0 < \delta < i_g(M)$ fixed,

$$ (4.32) \qquad \eta_\varepsilon(\delta) \lesssim \mu_\varepsilon^{\frac{n-2}{2}} $$

holds for $\varepsilon$ small enough. With (3.5), (4.28) and (4.32) this will conclude the proof of (4.27). We prove (4.32) by contradiction: we fix $0 < \delta < i_g(M)$ and, up to passing to a subsequence $(\varepsilon_k)_{k \ge 0}$ with $\lim_{k \to +\infty} \varepsilon_k = 0$ we assume that

$$ (4.33) \qquad \frac{\eta_{\varepsilon_k}(\delta)}{\mu_{\varepsilon_k}^{\frac{n-2}{2}}} \to +\infty $$

as $k \to +\infty$. We let, for $k \ge 0$ and for any $x \in M$,

$$ \hat{u}_k(x) = \frac{u_{\varepsilon_k}(x)}{\eta_{\varepsilon_k}(\delta)}, $$

and let $y_k \in M \setminus B_g(x_{\varepsilon_k}, \delta)$ be such that $|\hat{u}_k(y_k)| = 1$. By (3.5), (4.28) and (4.33) the family $(\hat{u}_k)_{k \ge 0}$ is uniformly bounded in compact subsets of $M \setminus \{x_0\}$, where $x_0 = \lim_{\varepsilon \to 0} x_\varepsilon$. By (1.3), standard elliptic theory and since $\eta_{\varepsilon_k}(\delta) \to 0$ as $k \to +\infty$, $\hat{u}_k$ thus converges in $C_{\text{loc}}^2(M \setminus \{x_0\})$, up to a subsequence, to some function $\hat{u}_0 \in C^2(M \setminus \{x_0\})$ that satisfies $\Delta_g \hat{u}_0 + h \hat{u}_0 = 0$ in $M \setminus \{x_0\}$ and $|\hat{u}_0(y_0)| = 1$, where $y_0 = \lim_{k \to +\infty} y_k \in M \setminus B_g(x_0, \delta)$. Passing (4.28) to the limit pointwise shows in addition that we have

$$ |\hat{u}_0(x)| \le C_\delta \quad \text{for all } x \in M \setminus \{x_0\}. $$

Classical arguments thus show that the singularity of $\hat{u}_0$ at $x_0$ is removable and that $\hat{u}_0$ satisfies $\Delta_g \hat{u}_0 + h \hat{u}_0 = 0$ in $M$ in a strong sense. This implies $\hat{u}_0 \equiv 0$ since $\Delta_g + h$ is coercive, which is a contradiction with $|\hat{u}_0(y_0)| = 1$. Thus (4.32) holds true and (4.27) is proven. $\square$

A consequence of (4.27), (1.3) and standard elliptic theory is the following estimate on the first and second derivatives of $u_\varepsilon$: for $0 \le k \le 2$,

$$ (4.34) \qquad |\nabla^k u_\varepsilon(x)| \lesssim \frac{\mu_\varepsilon^{\frac{n-2}{2}}}{(\mu_\varepsilon + d_g(x_\varepsilon, x))^{n-2+k}} \quad \text{for all } x \in M \text{ and } 0 < \varepsilon \le 1. $$

Estimate (4.34) follows from a scaling argument and from standard elliptic theory.

**Step 5:** we now claim that

$$ (4.35) \qquad 2^* - p_\varepsilon = \begin{cases} O(\mu_\varepsilon) & \text{if } n=3 \\ O\left(\mu_\varepsilon^2 \ln \frac{1}{\mu_\varepsilon}\right) & \text{if } n=4. \\ O(\mu_\varepsilon^2) & \text{if } n \ge 5 \end{cases}. $$
as $\varepsilon \to 0$. In particular $\lim_{\varepsilon \to 0} \mu_{\varepsilon}^{2^* - p_{\varepsilon}} = 1$.

*Proof of (4.35).* Let $0 < \delta < i_g(M)$ be fixed and let

$$ \Omega_{\varepsilon} = B_g(x_{\varepsilon}, \delta). $$

We write a Pohozaev identity for $u_{\varepsilon}$ on $\Omega_{\varepsilon}$: let $X$ be any smooth vector field in $M$. Integrating (1.3) by parts against $\langle X, \nabla u_{\varepsilon} \rangle$ classically shows that

$$ \begin{align} & \int_{\Omega_{\varepsilon}} h_{\varepsilon} u_{\varepsilon} \langle X, \nabla u_{\varepsilon} \rangle dv_g + \int_{\Omega_{\varepsilon}} \left( \nabla X - \frac{1}{2} \mathrm{div}_g X g \right) (\nabla u_{\varepsilon}, \nabla u_{\varepsilon}) dv_g \tag{4.36} \\ & \qquad + \int_{\partial \Omega_{\varepsilon}} \left( \frac{1}{2} \langle X, \nu \rangle |\nabla u_{\varepsilon}|_g^2 - \langle X, \nabla u_{\varepsilon} \rangle \partial_{\nu} u_{\varepsilon} - \frac{1}{p_{\varepsilon}} \langle X, \nu \rangle |u_{\varepsilon}|^{p_{\varepsilon}} \right) d\sigma_g \nonumber \\ & \qquad = - \frac{1}{p_{\varepsilon}} \int_{\Omega_{\varepsilon}} \mathrm{div}_g X |u_{\varepsilon}|^{p_{\varepsilon}} dv_g \nonumber \end{align} $$

(see e.g. [16, Proposition 6.2]). Let now $X_{\varepsilon}$ be the smooth vector field whose coordinates in the exponential chart at $x_{\varepsilon}$ are $(X_{\varepsilon}(x))^i = x^i$. For $x \in \Omega_{\varepsilon}$ and any $1 \le i, j \le n$ we have $\nabla_i X_{\varepsilon}^j(x) = \delta_i^j + O(d_g(x_{\varepsilon}, x)^2)$. As a consequence,

$$ \int_{\Omega_{\varepsilon}} \left( \nabla X - \frac{1}{2} \mathrm{div}_g X g \right) (\nabla u_{\varepsilon}, \nabla u_{\varepsilon}) dv_g = \int_{\Omega_{\varepsilon}} \left( 1 - \frac{n}{2} + O(d_g(x_{\varepsilon}, \cdot)^2) \right) |\nabla u_{\varepsilon}|_g^2 dv_g $$

and

$$ -\frac{1}{p_{\varepsilon}} \int_{\Omega_{\varepsilon}} \mathrm{div}_g X |u_{\varepsilon}|^{p_{\varepsilon}} dv_g = \int_{\Omega_{\varepsilon}} \left( -\frac{n}{p_{\varepsilon}} + O(d_g(x_{\varepsilon}, \cdot)^2) \right) |u_{\varepsilon}|^{p_{\varepsilon}} dv_g. $$

Integrating (1.3) by parts shows independently that

$$ \int_{\Omega_{\varepsilon}} |\nabla u_{\varepsilon}|_g^2 dv_g = \int_{\Omega_{\varepsilon}} |u_{\varepsilon}|^{p_{\varepsilon}} dv_g - \int_{\Omega_{\varepsilon}} h_{\varepsilon} u_{\varepsilon}^2 dv_g + \int_{\partial \Omega_{\varepsilon}} u_{\varepsilon} \partial_{\nu} u_{\varepsilon} d\sigma_g, $$

so that (4.36) becomes

$$ \begin{align*} & n \left( \frac{1}{2^*} - \frac{1}{p_{\varepsilon}} \right) \int_{\Omega_{\varepsilon}} |u_{\varepsilon}|^{p_{\varepsilon}} dv_g \\ & = \int_{\Omega_{\varepsilon}} \left( \frac{n-2}{2} u_{\varepsilon}^2 + h_{\varepsilon} u_{\varepsilon} \langle X, \nabla u_{\varepsilon} \rangle \right) dv_g \\ & \quad + O \left( \int_{\Omega_{\varepsilon}} d_g(x_{\varepsilon}, \cdot)^2 |\nabla u_{\varepsilon}|_g^2 dv_g + \int_{\Omega_{\varepsilon}} d_g(x_{\varepsilon}, \cdot)^2 |u_{\varepsilon}|^{p_{\varepsilon}} dv_g \right) \\ & \quad + \int_{\partial \Omega_{\varepsilon}} \left( \frac{1}{2} \langle X, \nu \rangle |\nabla u_{\varepsilon}|_g^2 - \langle X, \nabla u_{\varepsilon} \rangle \partial_{\nu} u_{\varepsilon} - \frac{n-2}{2} u_{\varepsilon} \partial_{\nu} u_{\varepsilon} \right. \\ & \qquad \left. - \frac{1}{p_{\varepsilon}} \langle X, \nu \rangle |u_{\varepsilon}|^{p_{\varepsilon}} \right) d\sigma_g. \end{align*} $$

Straightforward computations using (4.27) and (4.34) show that

$$ \int_{\Omega_{\varepsilon}} d_g(x_{\varepsilon}, \cdot)^2 |\nabla u_{\varepsilon}|_g^2 dv_g + \int_{\Omega_{\varepsilon}} d_g(x_{\varepsilon}, \cdot)^2 |u_{\varepsilon}|^{p_{\varepsilon}} dv_g \\ + \int_{\Omega_{\varepsilon}} |h_{\varepsilon} u_{\varepsilon} \langle X, \nabla u_{\varepsilon} \rangle| dv_g + \int_{\Omega_{\varepsilon}} |h_{\varepsilon}| u_{\varepsilon}^2 dv_g = \begin{cases} O(\mu_{\varepsilon}) & \text{if } n = 3 \\ O(\mu_{\varepsilon}^2 \ln \frac{1}{\mu_{\varepsilon}}) & \text{if } n = 4. \\ O(\mu_{\varepsilon}^2) & \text{if } n \ge 5 \end{cases}. $$
as $\varepsilon \to 0$. In particular $\lim_{\varepsilon \to 0} \mu_{\varepsilon}^{2^* - p_{\varepsilon}} = 1$.

*Proof of (4.35).* Let $0 < \delta < i_g(M)$ be fixed and let

$$ \Omega_{\varepsilon} = B_g(x_{\varepsilon}, \delta). $$

We write a Pohozaev identity for $u_{\varepsilon}$ on $\Omega_{\varepsilon}$: let $X$ be any smooth vector field in $M$. Integrating (1.3) by parts against $\langle X, \nabla u_{\varepsilon} \rangle$ classically shows that

$$ \begin{align} \tag{4.36} & \int_{\Omega_{\varepsilon}} h_{\varepsilon} u_{\varepsilon} \langle X, \nabla u_{\varepsilon} \rangle dv_g + \int_{\Omega_{\varepsilon}} \left( \nabla X - \frac{1}{2} \mathrm{div}_g X g \right) (\nabla u_{\varepsilon}, \nabla u_{\varepsilon}) dv_g \nonumber \\ & \qquad + \int_{\partial \Omega_{\varepsilon}} \left( \frac{1}{2} \langle X, \nu \rangle |\nabla u_{\varepsilon}|_g^2 - \langle X, \nabla u_{\varepsilon} \rangle \partial_{\nu} u_{\varepsilon} - \frac{1}{p_{\varepsilon}} \langle X, \nu \rangle |u_{\varepsilon}|^{p_{\varepsilon}} \right) d\sigma_g \nonumber \\ & = - \frac{1}{p_{\varepsilon}} \int_{\Omega_{\varepsilon}} \mathrm{div}_g X |u_{\varepsilon}|^{p_{\varepsilon}} dv_g \end{align} $$

(see e.g. [16, Proposition 6.2]). Let now $X_{\varepsilon}$ be the smooth vector field whose coordinates in the exponential chart at $x_{\varepsilon}$ are $(X_{\varepsilon}(x))^i = x^i$. For $x \in \Omega_{\varepsilon}$ and any $1 \le i, j \le n$ we have $\nabla_i X_{\varepsilon}^j(x) = \delta_i^j + O(d_g(x_{\varepsilon}, x)^2)$. As a consequence,

$$ \int_{\Omega_{\varepsilon}} \left( \nabla X - \frac{1}{2} \mathrm{div}_g X g \right) (\nabla u_{\varepsilon}, \nabla u_{\varepsilon}) dv_g = \int_{\Omega_{\varepsilon}} \left( 1 - \frac{n}{2} + O(d_g(x_{\varepsilon}, \cdot)^2) \right) |\nabla u_{\varepsilon}|_g^2 dv_g $$

and

$$ -\frac{1}{p_{\varepsilon}} \int_{\Omega_{\varepsilon}} \mathrm{div}_g X |u_{\varepsilon}|^{p_{\varepsilon}} dv_g = \int_{\Omega_{\varepsilon}} \left( -\frac{n}{p_{\varepsilon}} + O(d_g(x_{\varepsilon}, \cdot)^2) \right) |u_{\varepsilon}|^{p_{\varepsilon}} dv_g. $$

Integrating (1.3) by parts shows independently that

$$ \int_{\Omega_{\varepsilon}} |\nabla u_{\varepsilon}|_g^2 dv_g = \int_{\Omega_{\varepsilon}} |u_{\varepsilon}|^{p_{\varepsilon}} dv_g - \int_{\Omega_{\varepsilon}} h_{\varepsilon} u_{\varepsilon}^2 dv_g + \int_{\partial \Omega_{\varepsilon}} u_{\varepsilon} \partial_{\nu} u_{\varepsilon} d\sigma_g, $$

so that (4.36) becomes

$$ \begin{align*} & n \left( \frac{1}{2^*} - \frac{1}{p_{\varepsilon}} \right) \int_{\Omega_{\varepsilon}} |u_{\varepsilon}|^{p_{\varepsilon}} dv_g \\ &= \int_{\Omega_{\varepsilon}} \left( \frac{n-2}{2} u_{\varepsilon}^2 + h_{\varepsilon} u_{\varepsilon} \langle X, \nabla u_{\varepsilon} \rangle \right) dv_g \\ &\quad + O \left( \int_{\Omega_{\varepsilon}} d_g(x_{\varepsilon}, \cdot)^2 |\nabla u_{\varepsilon}|_g^2 dv_g + \int_{\Omega_{\varepsilon}} d_g(x_{\varepsilon}, \cdot)^2 |u_{\varepsilon}|^{p_{\varepsilon}} dv_g \right) \\ &\quad + \int_{\partial \Omega_{\varepsilon}} \left( \frac{1}{2} \langle X, \nu \rangle |\nabla u_{\varepsilon}|_g^2 - \langle X, \nabla u_{\varepsilon} \rangle \partial_{\nu} u_{\varepsilon} - \frac{n-2}{2} u_{\varepsilon} \partial_{\nu} u_{\varepsilon} \right. \\ &\qquad \left. - \frac{1}{p_{\varepsilon}} \langle X, \nu \rangle |u_{\varepsilon}|^{p_{\varepsilon}} \right) d\sigma_g. \end{align*} $$

Straightforward computations using (4.27) and (4.34) show that

$$ \int_{\Omega_{\varepsilon}} d_g(x_{\varepsilon}, \cdot)^2 |\nabla u_{\varepsilon}|_g^2 dv_g + \int_{\Omega_{\varepsilon}} d_g(x_{\varepsilon}, \cdot)^2 |u_{\varepsilon}|^{p_{\varepsilon}} dv_g \\ + \int_{\Omega_{\varepsilon}} |h_{\varepsilon} u_{\varepsilon} \langle X, \nabla u_{\varepsilon} \rangle| dv_g + \int_{\Omega_{\varepsilon}} |h_{\varepsilon}| u_{\varepsilon}^2 dv_g = \begin{cases} O(\mu_{\varepsilon}) & \text{if } n = 3 \\ O(\mu_{\varepsilon}^2 \ln \frac{1}{\mu_{\varepsilon}}) & \text{if } n = 4. \\ O(\mu_{\varepsilon}^2) & \text{if } n \ge 5 \end{cases}. $$
By (4.27) and (4.34) we have

$$ \int_{\partial \Omega_{\varepsilon}} \left( \frac{1}{2} \langle X, \nu \rangle |\nabla u_{\varepsilon}|_{g}^{2} - \langle X, \nabla u_{\varepsilon} \rangle \partial_{\nu} u_{\varepsilon} - \frac{n-2}{2} u_{\varepsilon} \partial_{\nu} u_{\varepsilon} - \frac{1}{p_{\varepsilon}} \langle X, \nu \rangle |u_{\varepsilon}|^{p_{\varepsilon}} \right) d\sigma_{g} = O(\mu_{\varepsilon}^{n-2}). $$

Plugging the latter computations in (4.36) thus shows that

$$ (4.37) \qquad n \left( \frac{1}{2^*} - \frac{1}{p_\varepsilon} \right) \int_{\Omega_\varepsilon} |u_\varepsilon|^{p_\varepsilon} dv_g = \begin{cases} O(\mu_\varepsilon) & \text{if } n=3 \\ O(\mu_\varepsilon^2 \ln \frac{1}{\mu_\varepsilon}) & \text{if } n=4. \\ O(\mu_\varepsilon^2) & \text{if } n \ge 5 \end{cases} $$

Finally, using (4.7) and (4.9), Fatou's lemma shows that

$$ \int_{\Omega_{\varepsilon}} |u_{\varepsilon}|^{p_{\varepsilon}} dv_g \geq c_0 $$

for some positive $c_0 > 0$ independent of $\varepsilon$. Going back to (4.37) thus proves (4.35). $\square$

**Step 6: End of the proof of Proposition 4.1.** We are now in position to conclude the proof of Proposition 4.1. The equality $c = 1$ follows from (4.6) and (4.35) and we thus only need to prove (4.5). We proceed by contradiction and assume that (4.5) is false: up to passing to a subsequence $(\varepsilon_k)_{k \ge 0}$ with $\lim_{k \to +\infty} \varepsilon_k = 0$, there exists a sequence $(y_k)_{k \ge 0}$ of points of $M$ such that

$$ (4.38) \qquad |u_{\varepsilon_k}(y_k) - B^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}(y_k)| \ge \eta_0 B_{\varepsilon_k}^{+}(y_k) $$

for some $\eta_0 > 0$, where $B^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}$ and $B_{\varepsilon_k}^+$ are as in (4.2) and (3.5). Since by (4.35) we have $c = 1$, (4.9) shows with (4.38) that

$$ (4.39) \qquad \frac{d_g(x_{\varepsilon_k}, y_{\varepsilon_k})}{\mu_{\varepsilon_k}} \to +\infty $$

as $k \to +\infty$. We write a representation formula for $u_{\varepsilon_k}$ in $M$: if $G_{h_{\varepsilon_k}}$ again denotes the Green's function for $\Delta_g + h_{\varepsilon_k}$, we have with (1.3):

$$ (4.40) \qquad u_{\varepsilon_k}(y_k) = \int_M G_{h_{\varepsilon_k}}(y_k, y) |u_{\varepsilon_k}(y)|^{p_{\varepsilon_k}-2} u_{\varepsilon_k}(y) dv_g(y). $$

Let $R > 1$ be fixed and let $x \in B(0, R)$. By (4.1) and (4.39) we have

$$ (n-2)\omega_{n-1} d_g(x_{\varepsilon_k}, y_k)^{n-2} G_{h_{\varepsilon_k}} \left( \exp_{x_{\varepsilon_k}} (\mu_{\varepsilon_k} x), y_k \right) \to 1 \quad \text{if } d_g(x_{\varepsilon_k}, y_k) \to 0 $$

uniformly in $B(0, R)$ as $k \to +\infty$. Since $h_{\varepsilon_k} \to h$ in $C^1(M)$ we also have

$$ G_{h_{\varepsilon_k}}(y_k, \exp_{x_{\varepsilon_k}}(\mu_{\varepsilon_k} x)) = G_h(x_{\varepsilon_k}, y_k) + o(1) \quad \text{if } d_g(x_{\varepsilon_k}, y_k) \not\to 0 $$
By (4.27) and (4.34) we have

$$ \int_{\partial \Omega_{\varepsilon}} \left( \frac{1}{2} \langle X, \nu \rangle |\nabla u_{\varepsilon}|_{g}^{2} - \langle X, \nabla u_{\varepsilon} \rangle \partial_{\nu} u_{\varepsilon} - \frac{n-2}{2} u_{\varepsilon} \partial_{\nu} u_{\varepsilon} - \frac{1}{p_{\varepsilon}} \langle X, \nu \rangle |u_{\varepsilon}|^{p_{\varepsilon}} \right) d\sigma_{g} = O(\mu_{\varepsilon}^{n-2}). $$

Plugging the latter computations in (4.36) thus shows that

$$ (4.37) \qquad n \left( \frac{1}{2^*} - \frac{1}{p_{\varepsilon}} \right) \int_{\Omega_{\varepsilon}} |u_{\varepsilon}|^{p_{\varepsilon}} dv_g = \begin{cases} O(\mu_{\varepsilon}) & \text{if } n=3 \\ O(\mu_{\varepsilon}^2 \ln \frac{1}{\mu_{\varepsilon}}) & \text{if } n=4 \\ O(\mu_{\varepsilon}^2) & \text{if } n \ge 5 \end{cases}. $$

Finally, using (4.7) and (4.9), Fatou's lemma shows that

$$ \int_{\Omega_{\varepsilon}} |u_{\varepsilon}|^{p_{\varepsilon}} dv_g \ge c_0 $$

for some positive $c_0 > 0$ independent of $\varepsilon$. Going back to (4.37) thus proves (4.35). $\square$

**Step 6: End of the proof of Proposition 4.1.** We are now in position to conclude the proof of Proposition 4.1. The equality $c=1$ follows from (4.6) and (4.35) and we thus only need to prove (4.5). We proceed by contradiction and assume that (4.5) is false: up to passing to a subsequence $(\varepsilon_k)_{k \ge 0}$ with $\lim_{k \to +\infty} \varepsilon_k = 0$, there exists a sequence $(y_k)_{k \ge 0}$ of points of $M$ such that

$$ (4.38) \qquad |u_{\varepsilon_k}(y_k) - B^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}(y_k)| \ge \eta_0 B_{\varepsilon_k}^{+}(y_k) $$

for some $\eta_0 > 0$, where $B^{\mu_{\varepsilon}, x_{\varepsilon}}$ and $B_{\varepsilon}^{+}$ are as in (4.2) and (3.5). By (4.35) we have that $c=1$. It follows from (4.9) that for any sequence $(z_k)_{k \in M}$ such that $d_g(x_{\varepsilon_k}, z_k) = O(\mu_{\varepsilon_k})$ for all $k$, we have that $|u_{\varepsilon_k}(z_k) - B^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}(z_k)| = o(B_{\varepsilon_k}^{+}(z_k))$ as $k \to +\infty$. Therefore, with (4.38), we get that

$$ (4.39) \qquad \frac{d_g(x_{\varepsilon_k}, y_{\varepsilon_k})}{\mu_{\varepsilon_k}} \to +\infty $$

as $k \to +\infty$. We write a representation formula for $u_{\varepsilon_k}$ in $M$: if $G_{h_{\varepsilon_k}}$ again denotes the Green's function for $\Delta_g + h_{\varepsilon_k}$, we have with (1.3):

$$ (4.40) \qquad u_{\varepsilon_k}(y_k) = \int_M G_{h_{\varepsilon_k}}(y_k, y) |u_{\varepsilon_k}(y)|^{p_{\varepsilon_k}-2} u_{\varepsilon_k}(y) dv_g(y). $$

Let $R > 1$ be fixed and let $x \in B(0, R)$. By (4.1) and (4.39) we have

$$ (n-2)\omega_{n-1} d_g(x_{\varepsilon_k}, y_k)^{n-2} G_{h_{\varepsilon_k}} \left( \exp_{x_{\varepsilon_k}} (\mu_{\varepsilon_k} x), y_k \right) \to 1 \quad \text{if } d_g(x_{\varepsilon_k}, y_k) \to 0 $$

uniformly in $B(0, R)$ as $k \to +\infty$. Since $h_{\varepsilon_k} \to h$ in $C^1(M)$ we also have

$$ G_{h_{\varepsilon_k}}(y_k, \exp_{x_{\varepsilon_k}}(\mu_{\varepsilon_k} x)) = G_h(x_{\varepsilon_k}, y_k) + o(1) \quad \text{if } d_g(x_{\varepsilon_k}, y_k) \not\to 0 $$
as $k \to +\infty$. Since $y_k$ satisfies (4.39), and by (4.1), the dominated convergence theorem gives with the definition (4.4)

$$
\begin{align}
& \int_{B_g(x_{\varepsilon_k}, R\mu_{\varepsilon_k})} G_{h_{\varepsilon_k}}(y_k, \cdot) |u_{\varepsilon_k}|^{p_{\varepsilon_k}-2} u_{\varepsilon_k} dv_g \nonumber \\
&= (1+o(1)) G_{h_{\varepsilon_k}}(x_{\varepsilon_k}, y_k) \mu_{\varepsilon_k}^{\frac{n-2}{2}} \int_{B(0,R)} |V|^{2^*-2} V dx \nonumber \\
&= (1+o(1)) G_{h_{\varepsilon_k}}(x_{\varepsilon_k}, y_k) (n-2) \omega_{n-1} \mu_{\varepsilon_k}^{\frac{n-2}{2}} (\lambda(V) + O(R^{-2})) \nonumber \\
&= B^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}(y_k) + o(B_{\varepsilon_k}(y_k)) + O(R^{-2} B_{\varepsilon_k}(y_k)) \tag{4.41}
\end{align}
$$

where $\lambda(V)$ is as in (2.4), where the third line follows from (2.8) and (2.9) and the last one from (4.2) and (4.39). Independently, straightforward computations show with (4.27) that

$$
\left| \int_{M \setminus B_g(x_{\varepsilon_k}, R\mu_{\varepsilon_k})} G_{h_{\varepsilon_k}}(y_k, \cdot) |u_{\varepsilon_k}|^{p_{\varepsilon_k}-2} u_{\varepsilon_k} dv_g \right| \lesssim R^{-2} B_{\varepsilon_k}^{+}(y_k) + o(B_{\varepsilon_k}^{+}(y_k))
$$

holds. Plugging the latter and (4.41) in (4.40) shows that, for any fixed $R > 0$ and for any $k$ large enough we have

$$
u_{\varepsilon_k}(y_k) = B^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}(y_k) + o(B_{\varepsilon_k}^{+}(y_k)) + O(R^{-2} B_{\varepsilon_k}^{+}(y_k)).
$$

By choosing $R$ large enough and passing to a subsequence we get a contradiction with (4.38). This concludes the proof of Proposition 4.1. $\square$

We conclude this section by showing an improvement of (4.34) when $k=1$ that will be needed in the next section. Let $(h_\epsilon)_{0<\epsilon\le 1}$, $(p_\epsilon)_{0<\epsilon\le 1}$ be as in the statement of Proposition 4.1 and let $(u_\epsilon)_{0<\epsilon\le 1}$ be a family of solutions of (1.3) satisfying (1.4) for some $c>0$. In particular Proposition 4.1 applies to $(u_\epsilon)_{0<\epsilon\le 1}$ and (4.5) holds. We claim that there exists a family of positive real numbers $(\sigma_\epsilon)_{0<\epsilon\le 1}$, with $\sigma_\epsilon \to 0$ as $\epsilon \to 0$, and a positive constant $C$ such that for any $x \in M$ and any $0 < \epsilon \le 1$ we have

$$
(4.42) \quad \left| \nabla u_{\epsilon}(x) - \nabla B^{\mu_{\epsilon}, x_{\epsilon}}(x) \right| \le \sigma_{\epsilon} \frac{\mu_{\epsilon}^{\frac{n-2}{2}}}{(\mu_{\epsilon} + d_g(x_{\epsilon}, x))^{n-1}} + C \mu_{\epsilon}^{\frac{n-2}{2}}
$$

where $B^{\mu_{\epsilon}, x_{\epsilon}}$ is as in (4.2).

Proof of (4.42). Let $0 < 2\delta < i_g(M)$ be fixed. We define $w_\epsilon(x) = u_\epsilon(\exp_{x_\epsilon}(x))$ for $x \in B(0, 2\delta)$. We first prove that there exists a family $(\sigma_\epsilon)_{0<\epsilon\le 1}$, with $\sigma_\epsilon \to 0$ as $\epsilon \to 0$, and a constant $C > 0$ such that for any $x \in B(0, 2\delta)$ and any $0 < \epsilon \le 1$ we have

$$
(4.43) \quad \left| \partial_t w_\epsilon(x) - \mu_\epsilon^{-\frac{n}{2}} (\partial_t V) \left( \frac{x}{\mu_\epsilon} \right) \right| \le \sigma_\epsilon \frac{\mu_\epsilon^{\frac{n-2}{2}}}{(\mu_\epsilon + |x|)^{n-1}} + C \mu_\epsilon^{\frac{n-2}{2}}.
$$

Let $x \in B(0, 2\delta)$. The representation formula (4.40) can be differentiated with respect to $x$: this gives, for $1 \le i \le n$,

$$
\partial_i w_\epsilon(x) = \int_M \partial_i G_{h_\epsilon} (\exp_{x_\epsilon}(x), y) |u_\epsilon(y)|^{p_\epsilon-2} u_\epsilon(y) dv_g(y),
$$

where we have let $\partial_i G_{h_\epsilon}(\exp_{x_\epsilon}(x), y) = \frac{\partial}{\partial x_i} (G_{h_\epsilon}(\exp_{x_\epsilon}(\cdot), y))(x)$. Standard results on the Green's functions (see e.g. [32]) show that there exists $C > 0$ such that for
as $k \to +\infty$. Since $y_k$ satisfies (4.39), and by (4.1), the dominated convergence theorem gives with the definition (4.4)

$$
\begin{align}
& \int_{B_g(x_{\varepsilon_k}, R\mu_{\varepsilon_k})} G_{h_{\varepsilon_k}}(y_k, \cdot) |u_{\varepsilon_k}|^{p_{\varepsilon_k}-2} u_{\varepsilon_k} dv_g \nonumber \\
&= (1+o(1)) G_{h_{\varepsilon_k}}(x_{\varepsilon_k}, y_k) \mu_{\varepsilon_k}^{\frac{n-2}{2}} \int_{B(0,R)} |V|^{2^*-2} V dx \nonumber \\
&= (1+o(1)) G_{h_{\varepsilon_k}}(x_{\varepsilon_k}, y_k) (n-2) \omega_{n-1} \mu_{\varepsilon_k}^{\frac{n-2}{2}} (\lambda(V) + O(R^{-2})) \nonumber \\
&= B^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}(y_k) + o(B_{\varepsilon_k}(y_k)) + O(R^{-2} B_{\varepsilon_k}(y_k)) \tag{4.41}
\end{align}
$$

where $\lambda(V)$ is as in (2.4), where the third line follows from (2.8) and (2.9) and the last one from (4.2) and (4.39). Independently, straightforward computations show with (4.27) that

$$
\left| \int_{M \setminus B_g(x_{\varepsilon_k}, R\mu_{\varepsilon_k})} G_{h_{\varepsilon_k}}(y_k, \cdot) |u_{\varepsilon_k}|^{p_{\varepsilon_k}-2} u_{\varepsilon_k} dv_g \right| \lesssim R^{-2} B_{\varepsilon_k}^{+}(y_k) + o(B_{\varepsilon_k}^{+}(y_k))
$$

holds. Plugging the latter and (4.41) in (4.40) shows that, for any fixed $R > 0$ and for any $k$ large enough we have

$$
u_{\varepsilon_k}(y_k) = B^{\mu_{\varepsilon_k}, x_{\varepsilon_k}}(y_k) + o(B_{\varepsilon_k}^{+}(y_k)) + O(R^{-2} B_{\varepsilon_k}^{+}(y_k)).
$$

By choosing $R$ large enough and passing to a subsequence we get a contradiction with (4.38). This concludes the proof of Proposition 4.1. $\square$

We conclude this section by showing an improvement of (4.34) when $k=1$ that will be needed in the next section. Let $(h_\varepsilon)_{0<\varepsilon\le 1}$, $(p_\varepsilon)_{0<\varepsilon\le 1}$ be as in the statement of Proposition 4.1 and let $(u_\varepsilon)_{0<\varepsilon\le 1}$ be a family of solutions of (1.3) satisfying (1.4) for some $c>0$. In particular Proposition 4.1 applies to $(u_\varepsilon)_{0<\varepsilon\le 1}$ and (4.5) holds. We claim that there exists a family of positive real numbers $(\sigma_\varepsilon)_{0<\varepsilon\le 1}$, with $\sigma_\varepsilon \to 0$ as $\varepsilon \to 0$, and a positive constant $C$ such that for any $x \in M$ and any $0 < \varepsilon \le 1$ we have

$$
(4.42) \quad \left| \nabla u_{\varepsilon}(x) - \nabla B^{\mu_{\varepsilon}, x_{\varepsilon}}(x) \right| \le \sigma_{\varepsilon} \frac{\mu_{\varepsilon}^{\frac{n-2}{2}}}{(\mu_{\varepsilon} + d_g(x_{\varepsilon}, x))^{n-1}} + C \mu_{\varepsilon}^{\frac{n-2}{2}}
$$

where $B^{\mu_{\varepsilon}, x_{\varepsilon}}$ is as in (4.2).

*Proof of* (4.42). Let $0 < 2\delta < i_g(M)$ be fixed. We define $w_\varepsilon(x) = u_\varepsilon(\exp_{x_\varepsilon}(x))$ for $x \in B(0, 2\delta)$. We first prove that there exists a family $(\sigma_\varepsilon)_{0<\varepsilon\le 1}$, with $\sigma_\varepsilon \to 0$ as $\varepsilon \to 0$, and a constant $C > 0$ such that for any $x \in B(0, 2\delta)$ and any $0 < \varepsilon \le 1$ we have

$$
(4.43) \quad \left| \partial_t w_\varepsilon(x) - \mu_\varepsilon^{-\frac{n}{2}} (\partial_t V)\left(\frac{x}{\mu_\varepsilon}\right) \right| \le \sigma_\varepsilon \frac{\mu_\varepsilon^{\frac{n-2}{2}}}{(\mu_\varepsilon + |x|)^{n-1}} + C \mu_\varepsilon^{\frac{n-2}{2}}.
$$

Let $x \in B(0, 2\delta)$. The representation formula (4.40) can be differentiated with respect to $x$: this gives, for $1 \le i \le n$,

$$
\partial_i w_\varepsilon(x) = \int_M \partial_i G_{h_\varepsilon} (\exp_{x_\varepsilon}(x), y) |u_\varepsilon(y)|^{p_\varepsilon-2} u_\varepsilon(y) dv_g(y),
$$

where we have let $\partial_i G_{h_\varepsilon}(\exp_{x_\varepsilon}(x), y) = \frac{\partial}{\partial x_i} (G_{h_\varepsilon}(\exp_{x_\varepsilon}(\cdot), y))(x)$. Standard results on the Green's functions (see e.g. [32]) show that there exists $C > 0$ such that for
all $\varepsilon$ small enough and all $x, y \in B(0, 2\delta)$,

$$ (4.44) \quad \left| \partial_i G_{h_\varepsilon} (\exp_{x_\varepsilon}(x), \exp_{x_\varepsilon}(y)) + \frac{1}{\omega_{n-1}} \frac{(x-y)_i}{|x-y|^n} \right| \le C|x-y|^{2-n} + C|x| |x-y|^{1-n} $$

holds. Let now, for $x \in \mathbb{R}^n$, $\tilde{V}_\varepsilon(x) = \mu_\varepsilon^{-\frac{n-2}{2}} V(\frac{x}{\mu_\varepsilon})$. Let $(z_\varepsilon)_{0<\varepsilon\le 1}$ be any family of points in $B(0, 2\delta)$. Straightforward computations using (4.35) and (4.5) show that

$$ \begin{aligned} \partial_i w_\varepsilon(z_\varepsilon) &= - \int_{B(0,2\delta)} \frac{1}{\omega_{n-1}} \frac{(z_\varepsilon - y)_i}{|z_\varepsilon - y|^n} |\tilde{V}_\varepsilon|^{p_\varepsilon - 2} \tilde{V}_\varepsilon dy + o\left(\frac{\mu_\varepsilon^{\frac{n-2}{2}}}{(\mu_\varepsilon + |z_\varepsilon|)^{n-1}}\right) + O(\mu_\varepsilon^{\frac{n-2}{2}}) \\ &= - \int_{B(0,2\delta)} \frac{1}{\omega_{n-1}} \frac{(z_\varepsilon - y)_i}{|z_\varepsilon - y|^n} |\tilde{V}_\varepsilon|^{2^* - 2} \tilde{V}_\varepsilon dy + o\left(\frac{\mu_\varepsilon^{\frac{n-2}{2}}}{(\mu_\varepsilon + |z_\varepsilon|)^{n-1}}\right) + O(\mu_\varepsilon^{\frac{n-2}{2}}) \\ &\quad + \int_{B(0,2\delta)} \frac{1}{\omega_{n-1}} \frac{(z_\varepsilon - y)_i}{|z_\varepsilon - y|^n} \left( |\tilde{V}_\varepsilon|^{2^* - 2} - |\tilde{V}_\varepsilon|^{p_\varepsilon - 2} \right) \tilde{V}_\varepsilon dy. \end{aligned} $$

Let us estimate the last term. Using the control (4.6), we get that

$$ \begin{aligned} & \left| \int_{B(0,2\delta)} \frac{1}{\omega_{n-1}} \frac{(z_\varepsilon - y)_i}{|z_\varepsilon - y|^n} \left( |\tilde{V}_\varepsilon|^{2^*-2} - |\tilde{V}_\varepsilon|^{p_\varepsilon-2} \right) \tilde{V}_\varepsilon dy \right| \\ & \le \int_{B(0,2\delta)} |z_\varepsilon - y|^{1-n} \left( |\tilde{V}_\varepsilon|^{2^*-p_\varepsilon} - 1 \right) \cdot |\tilde{V}_\varepsilon|^{p_\varepsilon-1} dy \\ & \le \int_{B(0,2\delta) \cap \{|\tilde{V}_\varepsilon| < \mu_\varepsilon^{\frac{n-2}{2}}\}} |z_\varepsilon - y|^{1-n} \left( |\tilde{V}_\varepsilon|^{2^*-p_\varepsilon} - 1 \right) \cdot |\tilde{V}_\varepsilon|^{p_\varepsilon-1} dy \\ & \quad + \int_{B(0,2\delta) \cap \{|\tilde{V}_\varepsilon| \ge \mu_\varepsilon^{\frac{n-2}{2}}\}} |z_\varepsilon - y|^{1-n} \left( |\tilde{V}_\varepsilon|^{2^*-p_\varepsilon} - 1 \right) \cdot |\tilde{V}_\varepsilon|^{p_\varepsilon-1} dy \\ & \le C \mu_\varepsilon^{\frac{n-2}{2}} (p_\varepsilon - 1) + \int_{B(0,2\delta) \cap \{|\tilde{V}_\varepsilon| \ge \mu_\varepsilon^{\frac{n-2}{2}}\}} |z_\varepsilon - y|^{1-n} \left( |\tilde{V}_\varepsilon|^{2^*-p_\varepsilon} - 1 \right) \cdot |\tilde{V}_\varepsilon|^{p_\varepsilon-1} dy \end{aligned} $$

The control (2.8) and (4.6) yield

$$ \{|\tilde{V}_\varepsilon| \ge \mu_\varepsilon^{\frac{n-2}{2}}\} \Rightarrow |\tilde{V}_\varepsilon|^{2^*-p_\varepsilon} - 1 = o(1) $$

as $\varepsilon \to 0$, so that we get

$$ \begin{aligned} \partial_i w_\varepsilon(z_\varepsilon) &= - \int_{\mathbb{R}^n} \frac{1}{\omega_{n-1}} \frac{(z_\varepsilon - y)_i}{|z_\varepsilon - y|^n} |\tilde{V}_\varepsilon|^{2^*-2} \tilde{V}_\varepsilon dy + o\left(\frac{\mu_\varepsilon^{\frac{n-2}{2}}}{(\mu_\varepsilon + |z_\varepsilon|)^{n-1}}\right) + O(\mu_\varepsilon^{\frac{n-2}{2}}) \\ &= \partial_i \tilde{V}_\varepsilon(z_\varepsilon) + o\left(\frac{\mu_\varepsilon^{\frac{n-2}{2}}}{(\mu_\varepsilon + |z_\varepsilon|)^{n-1}}\right) + O(\mu_\varepsilon^{\frac{n-2}{2}}), \end{aligned} $$

where the last equality follows from a representation formula for (1.2) in $\mathbb{R}^n$, since $V \in \Sigma$. This proves (4.43). Estimate (4.42) now follows from (4.1), (4.43), (4.44) and the explicit expression of $B^{\mu_\varepsilon, x_\varepsilon}$ in (4.2). $\square$

**Remark 4.2.** We assumed for simplicity here that $\Delta_g + h$ is coercive. If $\Delta_g + h$ had a kernel and, more generally, if $(u_\varepsilon)_{0<\varepsilon\le 1}$ was a finite-energy family of solutions of (1.3) (without the one-bubble assumption (1.4)), an analogue of Proposition 4.1 would still hold. We refer for instance to Premoselli [28] where a generalisation of Proposition 4.1 was proven when $p_\varepsilon = 2^*$ for all $0 < \varepsilon \le 1$, in the general
all $\varepsilon$ small enough and all $x, y \in B(0, 2\delta)$,

$$ (4.44) \quad \left| \partial_i G_{h_\varepsilon} (\exp_{x_\varepsilon}(x), \exp_{x_\varepsilon}(y)) + \frac{1}{\omega_{n-1}} \frac{(x-y)_i}{|x-y|^n} \right| \le C|x-y|^{2-n} + C|x| |x-y|^{1-n} $$

holds. Let now, for $x \in \mathbb{R}^n$, $\tilde{V}_\varepsilon(x) = \mu_\varepsilon^{-\frac{n-2}{2}} V(\frac{x}{\mu_\varepsilon})$. Let $(z_\varepsilon)_{0<\varepsilon\le 1}$ be any family of points in $B(0, 2\delta)$. Straightforward computations using (4.35) and (4.5) show that

$$ \begin{aligned} \partial_i w_\varepsilon(z_\varepsilon) &= - \int_{B(0,2\delta)} \frac{1}{\omega_{n-1}} \frac{(z_\varepsilon - y)_i}{|z_\varepsilon - y|^n} |\tilde{V}_\varepsilon|^{p_\varepsilon - 2} \tilde{V}_\varepsilon dy + o\left(\frac{\mu_\varepsilon^{\frac{n-2}{2}}}{(\mu_\varepsilon + |z_\varepsilon|)^{n-1}}\right) + O(\mu_\varepsilon^{\frac{n-2}{2}}) \\ &= - \int_{B(0,2\delta)} \frac{1}{\omega_{n-1}} \frac{(z_\varepsilon - y)_i}{|z_\varepsilon - y|^n} |\tilde{V}_\varepsilon|^{2^* - 2} \tilde{V}_\varepsilon dy + o\left(\frac{\mu_\varepsilon^{\frac{n-2}{2}}}{(\mu_\varepsilon + |z_\varepsilon|)^{n-1}}\right) + O(\mu_\varepsilon^{\frac{n-2}{2}}) \\ &\quad + \int_{B(0,2\delta)} \frac{1}{\omega_{n-1}} \frac{(z_\varepsilon - y)_i}{|z_\varepsilon - y|^n} \left( |\tilde{V}_\varepsilon|^{2^* - 2} - |\tilde{V}_\varepsilon|^{p_\varepsilon - 2} \right) \tilde{V}_\varepsilon dy. \end{aligned} $$

Let us estimate the last term. Using the control (4.6), we get that

$$ \begin{aligned} & \left| \int_{B(0,2\delta)} \frac{1}{\omega_{n-1}} \frac{(z_\varepsilon - y)_i}{|z_\varepsilon - y|^n} \left( |\tilde{V}_\varepsilon|^{2^*-2} - |\tilde{V}_\varepsilon|^{p_\varepsilon-2} \right) \tilde{V}_\varepsilon dy \right| \\ & \le \int_{B(0,2\delta)} |z_\varepsilon - y|^{1-n} \left( |\tilde{V}_\varepsilon|^{2^*-p_\varepsilon} - 1 \right) \cdot |\tilde{V}_\varepsilon|^{p_\varepsilon-1} dy \\ & \le \int_{B(0,2\delta) \cap \{|\tilde{V}_\varepsilon| < \mu_\varepsilon^{\frac{n-2}{2}}\}} |z_\varepsilon - y|^{1-n} \left( |\tilde{V}_\varepsilon|^{2^*-p_\varepsilon} - 1 \right) \cdot |\tilde{V}_\varepsilon|^{p_\varepsilon-1} dy \\ & \quad + \int_{B(0,2\delta) \cap \{|\tilde{V}_\varepsilon| \ge \mu_\varepsilon^{\frac{n-2}{2}}\}} |z_\varepsilon - y|^{1-n} \left( |\tilde{V}_\varepsilon|^{2^*-p_\varepsilon} - 1 \right) \cdot |\tilde{V}_\varepsilon|^{p_\varepsilon-1} dy \\ & \le C \mu_\varepsilon^{\frac{n-2}{2}} (p_\varepsilon - 1) + \int_{B(0,2\delta) \cap \{|\tilde{V}_\varepsilon| \ge \mu_\varepsilon^{\frac{n-2}{2}}\}} |z_\varepsilon - y|^{1-n} \left( |\tilde{V}_\varepsilon|^{2^*-p_\varepsilon} - 1 \right) \cdot |\tilde{V}_\varepsilon|^{p_\varepsilon-1} dy \end{aligned} $$

The control (2.8) and (4.6) yield

$$ \{|\tilde{V}_\varepsilon| \ge \mu_\varepsilon^{\frac{n-2}{2}}\} \Rightarrow |\tilde{V}_\varepsilon|^{2^*-p_\varepsilon} - 1 = o(1) $$

as $\varepsilon \to 0$, so that we get

$$ \begin{aligned} \partial_i w_\varepsilon(z_\varepsilon) &= - \int_{\mathbb{R}^n} \frac{1}{\omega_{n-1}} \frac{(z_\varepsilon - y)_i}{|z_\varepsilon - y|^n} |\tilde{V}_\varepsilon|^{2^*-2} \tilde{V}_\varepsilon dy + o\left(\frac{\mu_\varepsilon^{\frac{n-2}{2}}}{(\mu_\varepsilon + |z_\varepsilon|)^{n-1}}\right) + O(\mu_\varepsilon^{\frac{n-2}{2}}) \\ &= \partial_i \tilde{V}_\varepsilon(z_\varepsilon) + o\left(\frac{\mu_\varepsilon^{\frac{n-2}{2}}}{(\mu_\varepsilon + |z_\varepsilon|)^{n-1}}\right) + O(\mu_\varepsilon^{\frac{n-2}{2}}), \end{aligned} $$

where the last equality follows from a representation formula for (1.2) in $\mathbb{R}^n$, since $V \in \Sigma$. This proves (4.43). Estimate (4.42) now follows from (4.1), (4.43), (4.44) and the explicit expression of $B^{\mu_\varepsilon, x_\varepsilon}$ in (4.2). $\square$

**Remark 4.2.** We assumed for simplicity here that $\Delta_g + h$ is coercive. If $\Delta_g + h$ had a kernel and, more generally, if $(u_\varepsilon)_{0<\varepsilon\le 1}$ was a finite-energy family of solutions of (1.3) (without the one-bubble assumption (1.4)), an analogue of Proposition 4.1 would still hold. We refer for instance to Premoselli [28] where a generalisation of Proposition 4.1 was proven when $p_\varepsilon = 2^*$ for all $0 < \varepsilon \le 1$, in the general
multi-bubble case and when $\Delta_g + h$ is allowed to have a kernel. The proof of [28] allows to deal with any multi-bubble configuration and would still adapt to the asymptotically critical setting $p_\varepsilon \le 2^*$, $p_\varepsilon \to 2^*$ as a careful inspection of the proof shows. In the one-bubble case (1.4) that we consider here, the proof of Proposition 4.1 that we gave in this section is shorter and more direct.

5. PROOF OF THEOREMS 1.1 AND 1.2

Throughout this section we let $(M^n, g)$, $n \ge 3$ be a smooth, connected and closed manifold, $(h_\varepsilon)_{0<\varepsilon\le 1}$ be a family of $C^1$ functions in $M$ converging in $C^1(M)$ towards $h$ and let $(p_\varepsilon)_{0<\varepsilon\le 1}$ be a family of real numbers satisfying $2 < p_\varepsilon \le 2^*$ for all $0 < \varepsilon \le 1$ and such that $\lim_{\varepsilon\to 0} p_\varepsilon = 2^*$. We assume that $\Delta_g + h$ is coercive. We assume again, as in Section 4, that $(u_\varepsilon)_{0<\varepsilon\le 1}$ is a family of possibly sign-changing solutions of (1.3) satisfying (1.4):

$$u_{\varepsilon} = cB_{\varepsilon} + o(1) \quad \text{in } H^1(M),$$

for some $c > 0$, where $B_\varepsilon$ is a family of bubbles as in Definition 3.1 centered at $x_\varepsilon$ with radius $\mu_\varepsilon$ and modeled on $V \in \Sigma$. In particular the analysis of Section 4 applies: Proposition 4.1, (4.34) and (4.42) show that there exists a family $(\sigma_\varepsilon)_{0<\varepsilon\le 1}$ with $\sigma_\varepsilon \to 0$ as $\varepsilon \to 0$ and a positive constant $C$ such that, for any $x \in M$ and any $\varepsilon > 0$,

$$|u_{\varepsilon}(x) - B^{\mu_{\varepsilon}, x_{\varepsilon}}(x)| \le \sigma_{\varepsilon} B_{\varepsilon}^{+}(x),$$

$$ (5.1) \qquad \left| \nabla u_{\varepsilon}(x) - \nabla B^{\mu_{\varepsilon}, x_{\varepsilon}}(x) \right| \le \sigma_{\varepsilon} \frac{\mu_{\varepsilon}^{\frac{n-2}{2}}}{(\mu_{\varepsilon} + d_g(x_{\varepsilon}, x))^{n-1}} + C \mu_{\varepsilon}^{\frac{n-2}{2}} \quad \text{and} \\ (\mu_{\varepsilon} + d_g(x_{\varepsilon}, x))^{k} \left| \nabla^k u_{\varepsilon}(x) \right| \le C B_{\varepsilon}^{+}(x) \quad \text{for } k=1,2 $$

hold, where $B_\varepsilon^+$ and $B^{\mu_\varepsilon, x_\varepsilon}$ are as in (3.5) and (4.2). Equation (4.35) also applies and shows that $2^* - p_\varepsilon = O(\mu_\varepsilon)$. In this section we will prove Theorem 1.1 using (5.1). We first recall the classical conformal normal coordinates result of Lee-Parker [18]. It states the existence of a positive function $\Lambda \in C^\infty(M \times M)$ such that, letting $\Lambda_y(x) = \Lambda(y, x)$ for any $y, x \in M$ and $g_y := \Lambda_y^{\frac{n-2}{2}} g$ we have

$$ \Lambda_y(x) = 1 + O(d_{g_y}(y, x)^2), \\ \text{Ric}_{g_y}(x) = O(d_{g_y}(y, x)) \\ (5.2) \qquad S_{g_y} = O(d_{g_y}(y, x)^2) \quad \text{and} \\ \sqrt{|g_y|}(x) = 1 + O(d_{g_y}(y, x))^N $$

as $d_{g_y}(y, x) \to 0$, for a fixed $N$ large enough, where Ric$_{g_y}$ denotes the Ricci tensor of $g_y$ and where it is understood that $\sqrt{|g_y|}$ is computed with respect to the exponential map of $g_y$ at $y$. For $\varepsilon > 0$ we let $g_{x_\varepsilon} = \Lambda_{x_\varepsilon}^{\frac{n-2}{2}} g$ and we let $\exp_{x_\varepsilon}^{g_{x_\varepsilon}}$ be the exponential chart of $g_{x_\varepsilon}$ at $x_\varepsilon$. Let $\delta > 0$ be such that $2\delta < \inf_{\xi \in M} i_{g_\xi}(M)$. We define $g_\varepsilon$ and $v_\varepsilon$ in $B(0, 2\delta) \subset \mathbb{R}^n$ as follows:

$$ (5.3) \qquad g_{\varepsilon} = (\exp_{x_{\varepsilon}}^{g_{x_{\varepsilon}}})^{*} g_{x_{\varepsilon}} \quad \text{and} \\ v_{\varepsilon}(x) = \frac{u_{\varepsilon}}{\Lambda_{x_{\varepsilon}}} \circ \exp_{x_{\varepsilon}}^{g_{x_{\varepsilon}}}. $$
multi-bubble case and when $\Delta_g + h$ is allowed to have a kernel. The proof of [28] allows to deal with any multi-bubble configuration and would still adapt to the asymptotically critical setting $p_\varepsilon \le 2^*$, $p_\varepsilon \to 2^*$ as a careful inspection of the proof shows. In the one-bubble case (1.4) that we consider here, the proof of Proposition 4.1 that we gave in this section is shorter and more direct.

### 5. PROOF OF THEOREMS 1.1 AND 1.2

Throughout this section we let $(M^n, g)$, $n \ge 3$ be a smooth, connected and closed Riemannian manifold, $(h_\varepsilon)_{0<\varepsilon\le 1}$ be a family of $C^1$ functions in $M$ converging in $C^1(M)$ towards $h$ and let $(p_\varepsilon)_{0<\varepsilon\le 1}$ be a family of real numbers satisfying $2 < p_\varepsilon \le 2^*$ for all $0 < \varepsilon \le 1$ and such that $\lim_{\varepsilon\to 0} p_\varepsilon = 2^*$. We assume that $\Delta_g + h$ is coercive. We assume again, as in Section 4, that $(u_\varepsilon)_{0<\varepsilon\le 1}$ is a family of possibly sign-changing solutions of (1.3) satisfying (1.4):

$$u_\varepsilon = cB_\varepsilon + o(1) \quad \text{in } H^1(M),$$

for some $c > 0$, where $B_\varepsilon$ is a family of bubbles as in Definition 3.1 centered at $x_\varepsilon$ with radius $\mu_\varepsilon$ and modeled on $V \in \Sigma$. In particular the analysis of Section 4 applies: Proposition 4.1, (4.34) and (4.42) show that there exists a family $(\sigma_\varepsilon)_{0<\varepsilon\le 1}$ with $\sigma_\varepsilon \to 0$ as $\varepsilon \to 0$ and a positive constant $C$ such that, for any $x \in M$ and any $\varepsilon > 0$,

$$|u_\varepsilon(x) - B^{\mu_\varepsilon, x_\varepsilon}(x)| \le \sigma_\varepsilon B_\varepsilon^+(x),$$

$$ (5.1) \qquad \left| \nabla u_\varepsilon(x) - \nabla B^{\mu_\varepsilon, x_\varepsilon}(x) \right| \le \sigma_\varepsilon \frac{\mu_\varepsilon^{\frac{n-2}{2}}}{(\mu_\varepsilon + d_g(x_\varepsilon, x))^{n-1}} + C \mu_\varepsilon^{\frac{n-2}{2}} \quad \text{and} \\ (\mu_\varepsilon + d_g(x_\varepsilon, x))^k |\nabla^k u_\varepsilon(x)| \le C B_\varepsilon^+(x) \quad \text{for } k=1,2 $$

hold, where $B_\varepsilon^+$ and $B^{\mu_\varepsilon, x_\varepsilon}$ are as in (3.5) and (4.2). Equation (4.35) also applies and shows that $2^* - p_\varepsilon = O(\mu_\varepsilon)$. In this section we will prove Theorem 1.1 using (5.1). We first recall the classical conformal normal coordinates result of Lee-Parker [18]. It states the existence of a positive function $\Lambda \in C^\infty(M \times M)$ such that, letting $\Lambda_y(x) = \Lambda(y, x)$ for any $y, x \in M$ and $g_y := \Lambda_y^{\frac{n-2}{2}} g$ we have

$$ (5.2) \qquad \begin{aligned} \Lambda_y(x) &= 1 + O(d_{g_y}(y, x)^2), \\ \mathrm{Ric}_{g_y}(x) &= O(d_{g_y}(y, x)) \\ S_{g_y} &= O(d_{g_y}(y, x)^2) \quad \text{and} \\ \sqrt{|g_y|}(x) &= 1 + O(d_{g_y}(y, x))^N \end{aligned} $$

as $d_{g_y}(y, x) \to 0$, for a fixed $N$ large enough, where $\mathrm{Ric}_{g_y}$ denotes the Ricci tensor of $g_y$ and where it is understood that $\sqrt{|g_y|}$ is computed with respect to the exponential map of $g_y$ at $y$. For $\varepsilon > 0$ we let $g_{x_\varepsilon} = \Lambda_{x_\varepsilon}^{\frac{n-2}{2}} g$ and we let $\exp_{x_\varepsilon}^{g_{x_\varepsilon}}$ be the exponential chart of $g_{x_\varepsilon}$ at $x_\varepsilon$. Let $\delta > 0$ be such that $2\delta < \inf_{\varepsilon \in M} i_{g_\varepsilon}(M)$. We define $g_\varepsilon$ and $v_\varepsilon$ in $B(0, 2\delta) \subset \mathbb{R}^n$ as follows:

$$ (5.3) \qquad \begin{aligned} g_\varepsilon &= (\exp_{x_\varepsilon}^{g_{x_\varepsilon}})^* g_{x_\varepsilon} \quad \text{and} \\ v_\varepsilon(x) &= \frac{u_\varepsilon}{\Lambda_{x_\varepsilon}} \circ \exp_{x_\varepsilon}^{g_{x_\varepsilon}}. \end{aligned} $$
First, (5.2) shows that, for any fixed $x \in M$, we have

$$ (5.4) \qquad d_{g_x}(x, y) = d_g(x, y) + O(d_g(x, y)^3) \quad \text{as } d_g(x, y) \to 0. $$

The latter with (2.8), (4.2) and (5.1) implies that

$$ |v_{\varepsilon}(x) - \tilde{V}_{\varepsilon}(x)| \le \sigma_{\varepsilon} \tilde{B}_{\varepsilon}^{+}(x) + C \mu_{\varepsilon}^{\frac{n-2}{2}}, $$

$$ (5.5) \qquad \left| \nabla v_{\varepsilon}(x) - \nabla \tilde{V}_{\varepsilon}(x) \right| \le \sigma_{\varepsilon} \frac{\mu_{\varepsilon}^{\frac{n-2}{2}}}{(\mu_{\varepsilon} + |x|)^{n-1}} + C \mu_{\varepsilon}^{\frac{n-2}{2}} \quad \text{and} \\ (\mu_{\varepsilon} + d_{g_{\varepsilon}}(x_{\varepsilon}, x))^{k} |\nabla^{k} v_{\varepsilon}(x)| \le C \tilde{B}_{\varepsilon}^{+}(x) \quad \text{for } k = 1, 2 $$

for any $x \in B(0, 2\delta)$, where $C > 0$ is independent of $\varepsilon$ and where we have let

$$ \tilde{B}_{\varepsilon}^{+}(x) = \frac{\mu_{\varepsilon}^{\frac{n-2}{2}}}{\left(\mu_{\varepsilon}^{2} + \frac{|x|^2}{n(n-2)}\right)^{\frac{n-2}{2}}} \quad \text{and} \quad \tilde{V}_{\varepsilon}(x) = \mu_{\varepsilon}^{-\frac{n-2}{2}} V\left(\frac{x}{\mu_{\varepsilon}}\right). $$

Moreover, it follows from (4.9) that

$$ (5.6) \qquad \lim_{\varepsilon \to 0} \hat{v}_{\varepsilon} = V \quad \text{in } C_{\text{loc}}^{2}(\mathbb{R}^{n}), \quad \text{where } \hat{v}_{\varepsilon}(x) := \mu_{\varepsilon}^{\frac{n-2}{2}} v_{\varepsilon}(\mu_{\varepsilon}x). $$

Using (1.3) and the conformal invariance property of the conformal Laplacian it is
easily seen that $v_{\varepsilon}$ satisfies

$$ (5.7) \qquad \Delta_{\xi} v_{\varepsilon} + A_{\varepsilon}(v_{\varepsilon}) + \hat{h}_{\varepsilon} v_{\varepsilon} = \Lambda_{\varepsilon}^{p_{\varepsilon}-2*} |v_{\varepsilon}|^{p_{\varepsilon}-2} v_{\varepsilon} \quad \text{in } B(0, 2\delta), $$

where $\Lambda_{\varepsilon} := \Lambda_{x_{\varepsilon}} \circ \exp_{x_{\varepsilon}}^{g_{x_{\varepsilon}}}$, $\xi$ denotes the Euclidean metric and where we have let
$\Delta_{\xi} = - \sum_{i=1}^{n} \partial_{i}^{2}$,

$$ (5.8) \qquad \hat{h}_{\varepsilon}(x) = \left[ \Lambda_{x_{\varepsilon}}^{2-2*} \left( h_{\varepsilon} - \frac{n-2}{4(n-1)} S_g \right) \right] (\exp_{x_{\varepsilon}}^{g_{x_{\varepsilon}}}(x)) $$

and

$$ (5.9) \qquad A_{\varepsilon}(v_{\varepsilon}) = (\Delta_{g_{\varepsilon}} - \Delta_{\xi}) v_{\varepsilon} + \frac{n-2}{4(n-1)} S_{g_{\varepsilon}} v_{\varepsilon}. $$

Integrating (5.7) against $\langle x, \nabla v_{\varepsilon} \rangle + \frac{n-2}{2} v_{\varepsilon}$ in $B(0, \delta)$ with respect to the Euclidean metric yields the following Pohozaev identity (see Marques [21, Lemma 2.1]):

$$ \begin{align} & \int_{\partial B(0, \delta)} \left( \frac{n-2}{2} v_{\varepsilon} \partial_{\nu} v_{\varepsilon} - \frac{\delta}{2} |\nabla v_{\varepsilon}|_{\xi}^{2} + \delta(\partial_{\nu} v_{\varepsilon})^{2} + \frac{1}{p_{\varepsilon}} \Lambda_{\varepsilon}^{p_{\varepsilon}-2*} \delta |v_{\varepsilon}|^{p_{\varepsilon}} \right) d\sigma \nonumber \\ & = \frac{1}{p_{\varepsilon}} \int_{B(0, \delta)} (p_{\varepsilon} - 2^{*}) x^{k} \partial_{k} \Lambda_{\varepsilon} \Lambda_{x_{\varepsilon}}^{p_{\varepsilon}-1-2^{*}} |v_{\varepsilon}|^{p_{\varepsilon}} dx \nonumber \\ & \quad + n \left( \frac{1}{p_{\varepsilon}} - \frac{1}{2^{*}} \right) \int_{B(0, \delta)} \Lambda_{\varepsilon}^{p_{\varepsilon}-2^{*}} |v_{\varepsilon}|^{p_{\varepsilon}} dx \nonumber \\ & \quad + \int_{B(0, \delta)} \left( x^{k} \partial_{k} v_{\varepsilon} + \frac{n-2}{2} v_{\varepsilon} \right) A_{\varepsilon}(v_{\varepsilon}) dx \nonumber \\ & \quad + \int_{B(0, \delta)} \left( x^{k} \partial_{k} v_{\varepsilon} + \frac{n-2}{2} v_{\varepsilon} \right) \hat{h}_{\varepsilon} v_{\varepsilon} dx. \tag{5.10} \end{align} $$

In all the computations of this section it is intended that all the terms involving $p_{\varepsilon} - 2^{*}$ disappear in the case where $p_{\varepsilon} = 2^{*}$ for any $\varepsilon > 0$. We now estimate all the
First, (5.2) shows that, for any fixed $x \in M$, we have

$$ (5.4) \qquad d_{g_x}(x, y) = d_g(x, y) + O(d_g(x, y)^3) \quad \text{as } d_g(x, y) \to 0. $$

The latter with (2.8), (4.2) and (5.1) implies that

$$ |v_{\varepsilon}(x) - \tilde{V}_{\varepsilon}(x)| \le \sigma_{\varepsilon} \tilde{B}_{\varepsilon}^{+}(x) + C \mu_{\varepsilon}^{\frac{n-2}{2}}, $$

$$ (5.5) \qquad \begin{aligned} \left| \nabla v_{\varepsilon}(x) - \nabla \tilde{V}_{\varepsilon}(x) \right| &\le \sigma_{\varepsilon} \frac{\mu_{\varepsilon}^{\frac{n-2}{2}}}{(\mu_{\varepsilon} + |x|)^{n-1}} + C \mu_{\varepsilon}^{\frac{n-2}{2}} \quad \text{and} \\ (\mu_{\varepsilon} + d_{g_{\varepsilon}}(x_{\varepsilon}, x))^{k} |\nabla^{k} v_{\varepsilon}(x)| &\le C \tilde{B}_{\varepsilon}^{+}(x) \quad \text{for } k = 1, 2 \end{aligned} $$

for any $x \in B(0, 2\delta)$, where $C > 0$ is independent of $\varepsilon$ and where we have let

$$ \tilde{B}_{\varepsilon}^{+}(x) = \frac{\mu_{\varepsilon}^{\frac{n-2}{2}}}{\left(\mu_{\varepsilon}^{2} + \frac{|x|^2}{n(n-2)}\right)^{\frac{n-2}{2}}} \quad \text{and} \quad \tilde{V}_{\varepsilon}(x) = \mu_{\varepsilon}^{-\frac{n-2}{2}} V\left(\frac{x}{\mu_{\varepsilon}}\right). $$

Moreover, it follows from (4.9) that

$$ (5.6) \qquad \lim_{\varepsilon \to 0} \hat{v}_{\varepsilon} = V \quad \text{in } C_{\text{loc}}^{2}(\mathbb{R}^{n}), \text{ where } \hat{v}_{\varepsilon}(x) := \mu_{\varepsilon}^{\frac{n-2}{2}} v_{\varepsilon}(\mu_{\varepsilon} x). $$

Using (1.3) and the conformal invariance property of the conformal Laplacian it is
easily seen that $v_{\varepsilon}$ satisfies

$$ (5.7) \qquad \Delta_{\xi} v_{\varepsilon} + A_{\varepsilon}(v_{\varepsilon}) + \hat{h}_{\varepsilon} v_{\varepsilon} = \Lambda_{\varepsilon}^{p_{\varepsilon}-2*} |v_{\varepsilon}|^{p_{\varepsilon}-2} v_{\varepsilon} \quad \text{in } B(0, 2\delta), $$

where $\Lambda_{\varepsilon} := \Lambda_{x_{\varepsilon}} \circ \exp_{x_{\varepsilon}}^{g_{x_{\varepsilon}}}$, $\xi$ denotes the Euclidean metric and where we have let
$\Delta_{\xi} = - \sum_{i=1}^{n} \partial_{i}^{2}$,

$$ (5.8) \qquad \hat{h}_{\varepsilon}(x) = \left[ \Lambda_{x_{\varepsilon}}^{2-2*} \left( h_{\varepsilon} - \frac{n-2}{4(n-1)} S_g \right) \right] (\exp_{x_{\varepsilon}}^{g_{x_{\varepsilon}}}(x)) $$

and

$$ (5.9) \qquad A_{\varepsilon}(v_{\varepsilon}) = (\Delta_{g_{\varepsilon}} - \Delta_{\xi}) v_{\varepsilon} + \frac{n-2}{4(n-1)} S_{g_{\varepsilon}} v_{\varepsilon}. $$

Integrating (5.7) against $\langle x, \nabla v_{\varepsilon} \rangle + \frac{n-2}{2} v_{\varepsilon}$ in $B(0, \delta)$ with respect to the Euclidean
metric yields the following Pohozaev identity (see Marques [21, Lemma 2.1]):

$$ \begin{aligned} & \int_{\partial B(0, \delta)} \left( \frac{n-2}{2} v_{\varepsilon} \partial_{\nu} v_{\varepsilon} - \frac{\delta}{2} |\nabla v_{\varepsilon}|_{\xi}^{2} + \delta(\partial_{\nu} v_{\varepsilon})^{2} + \frac{1}{p_{\varepsilon}} \Lambda_{\varepsilon}^{p_{\varepsilon}-2*} \delta |v_{\varepsilon}|^{p_{\varepsilon}} \right) d\sigma \\ &= \frac{1}{p_{\varepsilon}} \int_{B(0, \delta)} (p_{\varepsilon} - 2^{*}) x^{k} \partial_{k} \Lambda_{\varepsilon} \Lambda_{x_{\varepsilon}}^{p_{\varepsilon}-1-2^{*}} |v_{\varepsilon}|^{p_{\varepsilon}} dx \\ & \quad + n \left( \frac{1}{p_{\varepsilon}} - \frac{1}{2^{*}} \right) \int_{B(0, \delta)} \Lambda_{\varepsilon}^{p_{\varepsilon}-2^{*}} |v_{\varepsilon}|^{p_{\varepsilon}} dx \\ & \quad + \int_{B(0, \delta)} \left( x^{k} \partial_{k} v_{\varepsilon} + \frac{n-2}{2} v_{\varepsilon} \right) A_{\varepsilon}(v_{\varepsilon}) dx \\ & \quad + \int_{B(0, \delta)} \left( x^{k} \partial_{k} v_{\varepsilon} + \frac{n-2}{2} v_{\varepsilon} \right) \hat{h}_{\varepsilon} v_{\varepsilon} dx. \end{aligned} \quad (5.10) $$

In all the computations of this section it is intended that all the terms involving $p_{\varepsilon} - 2^{*}$ disappear in the case where $p_{\varepsilon} = 2^{*}$ for any $\varepsilon > 0$. We now estimate all the
integrals appearing in (5.10). First, straightforward computations using (5.2) and (5.5) show that

$$
\begin{align}
& \frac{1}{p_{\varepsilon}} \int_{B(0,\delta)} (p_{\varepsilon} - 2^{*}) x^{k} \partial_{k} \Lambda_{\varepsilon} \Lambda_{\varepsilon}^{p_{\varepsilon}-1-2^{*}} |v_{\varepsilon}|^{p_{\varepsilon}} dx \nonumber \\
& \qquad + n \left( \frac{1}{p_{\varepsilon}} - \frac{1}{2^{*}} \right) \int_{B(0,\delta)} \Lambda_{\varepsilon}^{p_{\varepsilon}-2^{*}} |v_{\varepsilon}|^{p_{\varepsilon}} dx \nonumber \\
&= \frac{(n-2)^{2}}{4n} (2^{*} - p_{\varepsilon}) \int_{\mathbb{R}^{n}} |V|^{2^{*}} dx + o(|p_{\varepsilon} - 2^{*}|) \tag{5.11}
\end{align}
$$

and that

$$
\begin{equation}
\begin{split}
& \int_{\partial B(0, \delta)} \left( \frac{n-2}{2} v_{\varepsilon} \partial_{\nu} v_{\varepsilon} - \frac{\delta}{2} |\nabla v_{\varepsilon}|^{2}_{\xi} + \delta(\partial_{\nu} v_{\varepsilon})^{2} + \frac{1}{p_{\varepsilon}} \Lambda_{\varepsilon}^{p_{\varepsilon}-2^{*}} \delta |v_{\varepsilon}|^{p_{\varepsilon}} \right) d\sigma \\
& \quad = O(\mu_{\varepsilon}^{n-2}) = \begin{cases} o(\mu_{\varepsilon}^{2}) & \text{if } n \ge 5 \\ O(\mu_{\varepsilon}^{2}) & \text{if } n = 4 \end{cases}
\end{split}
\tag{5.12}
\end{equation}
$$

as $\epsilon \to 0$. A simple integration by parts shows that

$$
\begin{align}
& \int_{B(0,\delta)} \left( x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon \right) \hat{h}_\epsilon v_\epsilon dx \tag{5.13} \\
&= \frac{1}{2} \int_{\partial B(0,\delta)} \delta \hat{h}_\epsilon v_\epsilon^2 d\sigma - \int_{B(0,\delta)} \left( \hat{h}_\epsilon + \frac{1}{2} x^k \partial_k \hat{h}_\epsilon \right) v_\epsilon^2 dx, \nonumber
\end{align}
$$

so that with (5.5), we get that

$$
\begin{equation}
\begin{split}
& \int_{B(0,\delta)} \left( x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon \right) \hat{h}_\epsilon v_\epsilon dx \\
& \quad = - \int_{B(0,\delta)} \left( \hat{h}_\epsilon + \frac{1}{2} x^k \partial_k \hat{h}_\epsilon \right) v_\epsilon^2 dx + \begin{cases} o(\mu_\epsilon^2) & \text{if } n \ge 5 \\ O(\mu_\epsilon^2) & \text{if } n = 4 \end{cases} \\
& \quad = - \int_{B(0,\delta)} \hat{h}_\epsilon v_\epsilon^2 dx + \begin{cases} o(\mu_\epsilon^2) & \text{if } n \ge 5 \\ O(\mu_\epsilon^2) & \text{if } n = 4 \end{cases}
\end{split}
\tag{5.14}
\end{equation}
$$

Let now $w \in C^\infty(B(0, 2\delta))$ be a smooth function. We have, for any $\epsilon > 0$ and by
(5.2),

$$
\Delta_{g_{\epsilon}} w - \Delta_{\xi} w = -(g_{\epsilon}^{ij} - \delta^{ij}) \partial_{ij}^{2} w - \partial_{i}(g_{\epsilon})^{ij} \partial_{j} w + O(|x|^{N-1} |\nabla w|_{\xi})
$$

where repeated indices are summed over. By (5.2) we have Ric$_{g_\epsilon} (x_\epsilon) = 0$. Cartan's
expansion of $g_\epsilon$ and the symmetries of the Riemann tensor show that for any $1 \le$
$j \le n$

$$
\partial_i (g_\epsilon)^{ij}(x) = -\frac{1}{3} \mathrm{Ric}_{g_\epsilon} (x_\epsilon)_{pj} x^p + O(|x|^2) = O(|x|^2).
$$

By (5.2) again we also have $S_g(x_\epsilon) = 0$, so that $\mathrm{Rm}_{g_\epsilon}(x_\epsilon) = \mathrm{Weyl}_{g_\epsilon}(x_\epsilon)$, where $\mathrm{Weyl}_{g_\epsilon}$ is the (4,0) Weyl tensor of $g_\epsilon$. Using again Cartan's expansion of $g_\epsilon$ we find that for any smooth function $w$ in $B(0, 2\delta)$, we have

$$
(5.15) \quad \Delta_{g_{\epsilon}} w - \Delta_{\xi} w = \frac{1}{3} \mathrm{Weyl}_{g_{\epsilon}}(x_{\epsilon})_{ipqj} x^p x^q \partial_{ij}^2 w + O(|x|^3 |\nabla^2 w|) + O(|x|^2 |\nabla w|).
$$
integrals appearing in (5.10). First, straightforward computations using (5.2) and (5.5) show that

$$
\begin{align}
& \frac{1}{p_{\varepsilon}} \int_{B(0,\delta)} (p_{\varepsilon} - 2^{*}) x^{k} \partial_{k} \Lambda_{\varepsilon} \Lambda_{\varepsilon}^{p_{\varepsilon}-1-2^{*}} |v_{\varepsilon}|^{p_{\varepsilon}} dx \nonumber \\
& \qquad + n \left( \frac{1}{p_{\varepsilon}} - \frac{1}{2^{*}} \right) \int_{B(0,\delta)} \Lambda_{\varepsilon}^{p_{\varepsilon}-2^{*}} |v_{\varepsilon}|^{p_{\varepsilon}} dx \nonumber \\
&= \frac{(n-2)^{2}}{4n} (2^{*} - p_{\varepsilon}) \int_{\mathbb{R}^{n}} |V|^{2^{*}} dx + o(|p_{\varepsilon} - 2^{*}|) \tag{5.11}
\end{align}
$$

and that

$$
\begin{equation}
\begin{split}
& \int_{\partial B(0, \delta)} \left( \frac{n-2}{2} v_{\varepsilon} \partial_{\nu} v_{\varepsilon} - \frac{\delta}{2} |\nabla v_{\varepsilon}|^{2}_{\xi} + \delta(\partial_{\nu} v_{\varepsilon})^{2} + \frac{1}{p_{\varepsilon}} \Lambda_{\varepsilon}^{p_{\varepsilon}-2^{*}} \delta |v_{\varepsilon}|^{p_{\varepsilon}} \right) d\sigma \\
& \quad = O(\mu_{\varepsilon}^{n-2}) = \begin{cases} o(\mu_{\varepsilon}^{2}) & \text{if } n \ge 5 \\ O(\mu_{\varepsilon}^{2}) & \text{if } n = 4 \end{cases}
\end{split}
\tag{5.12}
\end{equation}
$$

as $\epsilon \to 0$. A simple integration by parts shows that

$$
\begin{align}
& \int_{B(0,\delta)} \left( x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon \right) \hat{h}_\epsilon v_\epsilon dx \tag{5.13} \\
&= \frac{1}{2} \int_{\partial B(0,\delta)} \delta \hat{h}_\epsilon v_\epsilon^2 d\sigma - \int_{B(0,\delta)} \left( \hat{h}_\epsilon + \frac{1}{2} x^k \partial_k \hat{h}_\epsilon \right) v_\epsilon^2 dx, \nonumber
\end{align}
$$

so that with (5.5), we get that

$$
\begin{equation}
\begin{split}
& \int_{B(0,\delta)} \left( x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon \right) \hat{h}_\epsilon v_\epsilon dx \\
& \quad = - \int_{B(0,\delta)} \left( \hat{h}_\epsilon + \frac{1}{2} x^k \partial_k \hat{h}_\epsilon \right) v_\epsilon^2 dx + \begin{cases} o(\mu_\epsilon^2) & \text{if } n \ge 5 \\ O(\mu_\epsilon^2) & \text{if } n = 4 \end{cases} \\
& \quad = - \int_{B(0,\delta)} \hat{h}_\epsilon v_\epsilon^2 dx + \begin{cases} o(\mu_\epsilon^2) & \text{if } n \ge 5 \\ O(\mu_\epsilon^2) & \text{if } n = 4 \end{cases}
\end{split}
\tag{5.14}
\end{equation}
$$

Let now $w \in C^\infty(B(0, 2\delta))$ be a smooth function. We have, for any $\epsilon > 0$ and by (5.2),

$$
\Delta_{g_{\epsilon}} w - \Delta_{\xi} w = -(g_{\epsilon}^{ij} - \delta^{ij}) \partial_{ij}^2 w - \partial_i (g_{\epsilon})^{ij} \partial_j w + O(|x|^{N-1} |\nabla w|_{\xi})
$$

where repeated indices are summed over. By (5.2) we have Ric$_{g_\epsilon} (x_\epsilon) = 0$. Cartan's expansion of $g_\epsilon$ and the symmetries of the Riemann tensor show that for any $1 \le j \le n$

$$
\partial_i (g_\epsilon)^{ij} (x) = -\frac{1}{3} \mathrm{Ric}_{g_\epsilon} (x_\epsilon)_{pj} x^p + O(|x|^2) = O(|x|^2).
$$

By (5.2) again we also have $S_g(x_\epsilon) = 0$, so that $\mathrm{Rm}_{g_\epsilon}(x_\epsilon) = \mathrm{Weyl}_{g_\epsilon}(x_\epsilon)$, where $\mathrm{Weyl}_{g_\epsilon}$ is the (4,0) Weyl tensor of $g_\epsilon$. Using again Cartan's expansion of $g_\epsilon$ we find that for any smooth function $w$ in $B(0, 2\delta)$, we have

$$
(5.15) \quad \Delta_{g_{\epsilon}} w - \Delta_{\xi} w = \frac{1}{3} \mathrm{Weyl}_{g_{\epsilon}}(x_{\epsilon})_{ipqj} x^p x^q \partial_{ij}^2 w + O(|x|^3 |\nabla^2 w|) + O(|x|^2 |\nabla w|).
$$
Applying the latter to $v_\epsilon$ and using (5.5) shows that

$$A_{\epsilon}(v_{\epsilon}) = \frac{1}{3} \mathrm{Weyl}_{g_{\epsilon}}(x_{\epsilon})_{ipqj} x^{p} x^{q} \partial_{ij}^{2} v_{\epsilon} + O\left( \frac{\mu_{\epsilon}^{\frac{n-2}{2}}}{(\mu_{\epsilon} + |x|)^{n-3}} \right).$$

Integrating the latter against $x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon$ thus gives, using again (5.5), that

$$\begin{align} \int_{B(0,\delta)} \left( x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon \right) A_\epsilon(v_\epsilon) dx \tag{5.16} \\ = \frac{1}{3} \mathrm{Weyl}_{g_\epsilon}(x_\epsilon)_{ipqj} \int_{B(0,\delta)} x^p x^q \partial_{ij}^2 v_\epsilon \left( x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon \right) dx \nonumber \\ + \begin{cases} o(\mu_\epsilon^2) & \text{if } n \ge 5 \\ O(\mu_\epsilon^2) & \text{if } n = 4 \end{cases} \nonumber \end{align}$$

where again repeated indices are summed over. We now separate between the cases $n=3$, $n=4$ and $n \ge 5$ and prove Theorems 1.1 and 1.2.

5.1. The case $n \ge 5$. With $\hat{v}_\epsilon$ and $\hat{h}_\epsilon$ as in (5.6) and (5.8) a change of variables gives

$$\int_{B(0,\delta)} \hat{h}_{\epsilon} dx = \mu_{\epsilon}^{2} \int_{B(0, \mu_{\epsilon}^{-1} \delta)} \hat{h}_{\epsilon}(\mu_{\epsilon} x) \hat{v}_{\epsilon}^{2} dx.$$

It follows from (5.5) that $|\hat{v}_\epsilon(x)| \le C(1+|x|)^{2-n} \in L^2(\mathbb{R}^n)$ for $n \ge 5$. Therefore, using (5.6) and Lebesgue's convergence theorem, we get that (5.14) yields

$$\begin{align} \int_{B(0,\delta)} \left( x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon \right) \hat{h}_\epsilon v_\epsilon dx \tag{5.17} \\ = - \left( h(x_0) - \frac{n-2}{4(n-1)} S_g(x_0) \right) \int_{\mathbb{R}^n} V^2 dx \cdot \mu_\epsilon^2 + o(\mu_\epsilon^2) \nonumber \end{align}$$

as $\epsilon \to 0$, where $x_0 = \lim_{\epsilon \to 0} x_\epsilon$. This limit exists up to passing to a subsequence as $\epsilon \to 0$, which will be implicit from now on. Note that the integral $\int_{\mathbb{R}^n} V^2 dx$ is finite by (2.8) since $n \ge 5$.

We now deal with (5.16), still for $n \ge 5$. By construction $g_\epsilon$ converges towards $(\exp^{g_{x_0}}) g_{x_0}$ in $C^1(B(0, \delta))$ as $\epsilon \to 0$, so that $\mathrm{Weyl}_{g_\epsilon}(x_\epsilon)_{ipqj} \to \mathrm{Weyl}_{g_{x_0}}(x_0)_{ipqj}$ as $\epsilon \to 0$, where the components of $\mathrm{Weyl}_{g_{x_0}}$ are measured in the exponential chart of $g_{x_0}$ at $x_0$. Since $g_{x_0}$ is conformal to $g$ and since $\Lambda_{x_0}(x_0) = 1$ we have $g_{x_0}(x_0) = g(x_0)$, and by the conformal invariance of the Weyl tensor we also have

$$\mathrm{Weyl}_{g_{x_0}}(x_0)_{ipqj} = \mathrm{Weyl}_g(x_0)_{ipqj},$$

where the components of $\mathrm{Weyl}_g(x_0)$ in the right-hand side are measured in the exponential chart of $g$ at $x_0$. By (5.5) and the dominated convergence theorem we therefore obtain, with (5.16), that

$$\begin{align} \int_{B(0,\delta)} \left( x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon \right) A_\epsilon(v_\epsilon) dx \tag{5.18} \\ = \frac{1}{3} \mathrm{Weyl}_g(x_0)_{ipqj} \int_{\mathbb{R}^n} x^p x^q \partial_{ij}^2 V \left( x^k \partial_k V + \frac{n-2}{2} V \right) dx \cdot \mu_\epsilon^2 + o(\mu_\epsilon^2) \nonumber \\ = - \frac{(n-2)^2}{4n} \left( \int_{\mathbb{R}^n} |V|^{2*} dx \right) \mathrm{Weyl}_g \otimes B \cdot \mu_\epsilon^2 + o(\mu_\epsilon^2) \nonumber \end{align}$$
Applying the latter to $v_\epsilon$ and using (5.5) shows that

$$A_{\epsilon}(v_{\epsilon}) = \frac{1}{3} \mathrm{Weyl}_{g_{\epsilon}}(x_{\epsilon})_{ipqj} x^{p} x^{q} \partial_{ij}^{2} v_{\epsilon} + O\left( \frac{\mu_{\epsilon}^{\frac{n-2}{2}}}{(\mu_{\epsilon} + |x|)^{n-3}} \right).$$

Integrating the latter against $x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon$ thus gives, using again (5.5), that

$$\begin{align} \int_{B(0,\delta)} \left( x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon \right) A_\epsilon(v_\epsilon) dx \tag{5.16} \\ = \frac{1}{3} \mathrm{Weyl}_{g_\epsilon}(x_\epsilon)_{ipqj} \int_{B(0,\delta)} x^p x^q \partial_{ij}^2 v_\epsilon \left( x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon \right) dx \nonumber \\ + \begin{cases} o(\mu_\epsilon^2) & \text{if } n \ge 5 \\ O(\mu_\epsilon^2) & \text{if } n = 4 \end{cases} \nonumber \end{align}$$

where again repeated indices are summed over. We now separate between the cases $n=3$, $n=4$ and $n \ge 5$ and prove Theorems 1.1 and 1.2.

5.1. The case $n \ge 5$. With $\hat{v}_\epsilon$ and $\hat{h}_\epsilon$ as in (5.6) and (5.8) a change of variables gives

$$\int_{B(0,\delta)} \hat{h}_{\epsilon} dx = \mu_{\epsilon}^{2} \int_{B(0, \mu_{\epsilon}^{-1} \delta)} \hat{h}_{\epsilon}(\mu_{\epsilon} x) \hat{v}_{\epsilon}^{2} dx.$$

It follows from (5.5) that $|\hat{v}_\epsilon(x)| \le C(1+|x|)^{2-n} \in L^2(\mathbb{R}^n)$ for $n \ge 5$. Therefore, using (5.6) and Lebesgue's convergence theorem, we get that (5.14) yields

$$\begin{align} \int_{B(0,\delta)} \left( x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon \right) \hat{h}_\epsilon v_\epsilon dx \tag{5.17} \\ = - \left( h(x_0) - \frac{n-2}{4(n-1)} S_g(x_0) \right) \int_{\mathbb{R}^n} V^2 dx \cdot \mu_\epsilon^2 + o(\mu_\epsilon^2) \nonumber \end{align}$$

as $\epsilon \to 0$, where $x_0 = \lim_{\epsilon \to 0} x_\epsilon$. This limit exists up to passing to a subsequence as $\epsilon \to 0$, which will be implicit from now on. Note that the integral $\int_{\mathbb{R}^n} V^2 dx$ is finite by (2.8) since $n \ge 5$.

We now deal with (5.16), still for $n \ge 5$. By construction $g_\epsilon$ converges towards $(\exp^{g_{x_0}})^* g_{x_0}$ in $C^1(B(0, \delta))$ as $\epsilon \to 0$, so that $\mathrm{Weyl}_{g_\epsilon}(x_\epsilon)_{ipqj} \to \mathrm{Weyl}_{g_{x_0}}(x_0)_{ipqj}$ as $\epsilon \to 0$, where the components of $\mathrm{Weyl}_{g_{x_0}}$ are measured in the exponential chart of $g_{x_0}$ at $x_0$. Since $g_{x_0}$ is conformal to $g$ and since $\Lambda_{x_0}(x_0) = 1$ we have $g_{x_0}(x_0) = g(x_0)$, and by the conformal invariance of the Weyl tensor we also have

$$\mathrm{Weyl}_{g_{x_0}}(x_0)_{ipqj} = \mathrm{Weyl}_g(x_0)_{ipqj},$$

where the components of $\mathrm{Weyl}_g(x_0)$ in the right-hand side are measured in the exponential chart of $g$ at $x_0$. By (5.5) and the dominated convergence theorem we therefore obtain, with (5.16), that

$$\begin{align} \int_{B(0,\delta)} \left( x^k \partial_k v_\epsilon + \frac{n-2}{2} v_\epsilon \right) A_\epsilon(v_\epsilon) dx \tag{5.18} \\ = \frac{1}{3} \mathrm{Weyl}_g(x_0)_{ipqj} \int_{\mathbb{R}^n} x^p x^q \partial_{ij}^2 V \left( x^k \partial_k V + \frac{n-2}{2} V \right) dx \cdot \mu_\epsilon^2 + o(\mu_\epsilon^2) \nonumber \\ = - \frac{(n-2)^2}{4n} \left( \int_{\mathbb{R}^n} |V|^{2*} dx \right) \mathrm{Weyl}_g \otimes B \cdot \mu_\epsilon^2 + o(\mu_\epsilon^2) \nonumber \end{align}$$