# Complete metrics with constant fractional higher order $Q$-curvature on the
punctured sphere

João Henrique Andrade  ,  Juncheng Wei  and  Zikai Ye  Department of
Mathematics, University of British Columbia
V6T 1Z2, Vancouver-BC, Canada
and
Institute of Mathematics and Statistics, University of São Paulo
05508-090, São Paulo-SP, Brazil<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Mathematics, University of British Columbia
V6T 1Z2, Vancouver-BC, Canada<EMAIL_ADDRESS>Department of Mathematics,
University of British Columbia
V6T 1Z2, Vancouver-BC, Canada<EMAIL_ADDRESS>
###### Abstract.

This manuscript is devoted to constructing complete metrics with constant
higher fractional curvature on punctured spheres with finitely many isolated
singularities. Analytically, this problem is reduced to constructing singular
solutions for a conformally invariant integro-differential equation that
generalizes the critical GJMS problem. Our proof follows the earlier
construction in Ao et al. [36], based on a gluing method, which we briefly
describe. Our main contribution is to provide a unified approach for
fractional and higher order cases. This method relies on proving Fredholm
properties for the linearized operator around a suitably chosen approximate
solution. The main challenge in our approach is that the solutions to the
related blow-up limit problem near isolated singularities need to be fully
classified; hence we are not allowed to use a simplified ODE method. To
overcome this issue, we approximate solutions near each isolated singularity
by a family of half-bubble tower solutions. Then, we reduce our problem to
solving an (infinite-dimensional) Toda-type system arising from the
interaction between the bubble towers at each isolated singularity. Finally,
we prove that this system’s solvability is equivalent to the existence of a
balanced configuration.

###### Key words and phrases:

Fractional poly-Laplacian, Higher order PDEs, GJMS operators, Critical
exponent, Gluing technique, Toda systems

###### 2020 Mathematics Subject Classification:

35J60, 35B09, 35J30, 35B40, 35R11

This work was partially supported by São Paulo Research Foundation (FAPESP)
#2020/07566-3 and #2021/15139-0 and Natural Sciences and Engineering Research
Council of Canada (NSERC)

###### Contents

  1. 1 Introduction
  2. 2 Notations
  3. 3 Dual representation formula
  4. 4 Delaunay-type solutions
    1. 4.1 Integral Emden–Fowler coordinates
    2. 4.2 Asymptotic classification of Delaunay-type solutions
  5. 5 Approximate solution
    1. 5.1 Local asymptotic behavior
    2. 5.2 Balanced configurations
    3. 5.3 Admissible perturbation parameters
    4. 5.4 Generalized Delaunay solutions
    5. 5.5 Normalized approximate kernels
    6. 5.6 Weighted functional spaces
    7. 5.7 Perturbation of the approximate solution
      1. 5.7.1 Quantitative estimates
      2. 5.7.2 Finite-dimensional reduction
      3. 5.7.3 Fixed-point argument
  6. 6 Estimates for the projections on the approximate null space
    1. 6.1 Projection on the normalized approximate kernels
    2. 6.2 Derivative of the projection on the normalized approximate kernels
  7. 7 Gluing technique
    1. 7.1 Infinite-dimensional Toda-system
    2. 7.2 Proof of the main result
  8. A Estimates on the bubble-towers interactions
  9. B Nondegeneracy of the bubble solution

##  1\. Introduction

The problem of constructing complete metrics on punctured spheres with
prescribed fractional higher order curvature is longstanding in differential
geometry. In [32], Graham, Jenne, Mason, and Sparling constructed conformally
covariant differential operators $P_{2m}(g)$ on a given compact
$n$-dimensional Riemannian manifold $(M^{n},g)$ for any $m\in\mathbb{N}$ such
that the leading order term of $P_{2m}(g)$ is $(-\Delta_{g})^{m}$ with $n>2m$.
One can then construct the associated $Q$-curvature of order $2m$ by
$Q_{2m}(g)=P_{2m}(g)(1)$. When $m=1$, one recovers the conformal Laplacian

$P_{2}(g)=-\Delta_{g}+\frac{n-2}{4(n-1)}R_{g}\quad{\rm with}\quad
Q_{2}(g)=\frac{n-2}{4(n-1)}R_{g},$

where $\Delta_{g}$ is the Laplace-Beltrami operator of $g$ and $R_{g}$ is its
scalar curvature. We also refer to [3, Appendix A] for the explicit formulae
for $P_{2}(g)$, $P_{4}(g)$ and $P_{6}(g)$. Subsequently, Grahan and Zworski
[33] and Chang and González [22] extended these definitions in the case the
background metric is the round metric on the sphere to obtain (nonlocal)
operators $P_{2\sigma}(g)$ of any order $\sigma\in(0,\frac{n}{2})$ as well as
its corresponding $Q$-curvature. Once again, the leading order part of
$P_{2\sigma}(g)$ is $(-\Delta_{g})^{\sigma}$, understood as the principal
value of a singular integral operator.

Nevertheless, the expressions for $P_{2\sigma}(g)$ and $Q_{2\sigma}(g)$ for a
general $\sigma\in\mathbb{R}_{+}$ are far more complicated. Namely, the
fractional curvature $Q_{2\sigma}(g)$ is defined from the conformal fractional
Laplacian $P_{2\sigma}(g)$ as $Q_{2\sigma}(g)=P_{2\sigma}(g)(1)$. It is a
nonlocal version of the scalar curvature (corresponding to the local case
$\sigma=1$). The conformal higher order fractional Laplacian $P_{2\sigma}(g)$
is a (nonlocal) pseudo-differential operator of order $2\sigma$, which can be
constructed from scattering theory on the conformal infinity $M^{n}$ of a
conformally compact Einstein manifold $(X^{n+1},g^{+})$ as a generalized
Dirichlet-to-Neumann operator for the eigenvalue problem

$-\Delta_{g^{+}}U-\frac{(n+2\sigma)^{2}}{4}U=0\quad\text{ in }\quad X,$

where $U\in\mathcal{C}^{\infty}(X)$ is the respective extension of
$u\in\mathcal{C}^{\infty}(M)$ This construction is a natural one from the
point of view of the AdS/CFT correspondence in Physics, also known as
Maldacena’s duality [45]. We refer the reader to [1, 55] for more details.

In this manuscript, we are restricted to the $n$-dimensional sphere
$\mathbb{S}^{n}\subset\mathbb{R}^{n+1}$, where $n>2\sigma$ and
$\sigma\in(1,+\infty)$ equipped with the standard round metric $g_{0}$, which
is given by the pullback of the usual Euclidean metric $\delta$ under the
stereographic projection
$\Pi:\mathbb{S}^{n}\setminus\\{\mathrm{e}_{1}\\}\rightarrow\mathbb{R}^{n}\setminus\\{0\\}$
with $\mathrm{e}_{1}=(1,0,\dots,0)\in\mathbb{S}^{n}$ denoting its north pole.
For any $k\in\mathbb{R}$ with $0\leqslant k\leqslant n$, we seek complete
metrics on $\mathbb{S}^{n}\setminus\Lambda^{k}$ of the form
$g=\mathrm{u}^{{4}/{(n-2\sigma)}}g_{0}$, where $\Lambda\subset\mathbb{S}^{n}$
is such that $\\#\Sigma=N$. In order to $g$ to be complete on
$\mathbb{S}^{n}\setminus\Lambda$, one has to impose
$\liminf_{\mathrm{d}(p,\Lambda)}\mathrm{u}(p)=+\infty$. Also, we prescribe the
resulting metric to have constant $Q_{2\sigma}$-curvature, which we normalize
to be

$Q_{n,\sigma}=Q_{2\sigma}({g_{0}})=\Gamma\left(\frac{n+2s}{2}\right)\Gamma\left(\frac{n-2s}{2}\right)^{-1},$

where $\Gamma(z)=\int_{0}^{\infty}\tau^{z-1}e^{-\tau}\mathrm{d}\tau$ is the
standard Gamma function.

Let us now introduce some standard terminology. For any $\sigma\in(1,+\infty]$
with $n>2\sigma$ and $N\geqslant 2$, we respectively denote by

$\displaystyle\mathcal{M}_{2\sigma,\Lambda}(g_{0})=\left\\{g\in[g_{0}]:\mbox{$g$
is complete on
$\mathbb{S}^{n}\setminus\Lambda$}\;\mbox{and}\;Q_{2\sigma}(g)\equiv
Q_{n,\sigma}\right\\}$

and

$\displaystyle\mathcal{M}_{2\sigma,N}(g_{0})=\left\\{g\in[g_{0}]:\mbox{$g$ is
complete on $\mathbb{S}^{n}\setminus\Lambda$ with
$\\#\Lambda=N$}\;\mbox{and}\;Q_{2\sigma}(g)\equiv Q_{n,\sigma}\right\\}$ (1)

the marked and unmarked moduli spaces of complete constant higher order
fractional $Q$-curvature metrics with isolated singularities. We also denote
by ${\rm sing}({g})=\Lambda$ its respective singular set.

In this fashion, our main theorem in this paper is the following

######  Theorem 1.

Let $\sigma\in(1,+\infty)$ with $n>2\sigma$. For any configuration
$\Lambda\subset\mathbb{S}^{n}$ such that $\\#\Lambda=N$ with $N\geqslant 2$,
there exists a metric $g\in\mathcal{M}_{2\sigma,N}(g_{0})$ satisfying ${\rm
sing}(g)=\Lambda$ and is unmarked nondegenerate. For a generic set of
$\Lambda=\\{p_{1},\ldots,p_{N}\\}$, this solution is marked nondegenerate, and
for such a metric
$(p_{1},\ldots,p_{N},\varepsilon_{1},\ldots,\varepsilon_{N})\in\mathbb{R}^{N(n+1)}$
constitute a full set of coordinates in $\mathcal{M}_{2\sigma,N}(g_{0})$ near
$g_{0}$. In particular, one has
$\mathcal{M}_{2\sigma,N}(g_{0})\neq\varnothing$.

Let us derive an analytical formulation for our main result. The family of
higher order fractional curvatures transform nicely under a conformal change.
Indeed, for any $\bar{g}\in[g]$, one has

$Q_{2\sigma}({\bar{g}})=\frac{2}{n-2\sigma}u^{-\frac{n+2\sigma}{n-2\sigma}}P_{2\sigma}({g_{0}})u,$

where
$P_{2\sigma}({g_{0}}):\mathcal{C}^{\infty}(M)\rightarrow\mathcal{C}^{\infty}(M)$
is the fractional higher order GJMS operator on the sphere

$P_{2\sigma}({g_{0}}):=\Gamma\left(\sqrt{-\Delta_{g_{0}}+\frac{(n-1)^{2}}{4}}+2\sigma+\frac{1}{2}\right)\Gamma\left(\sqrt{-\Delta_{g_{0}}+\frac{(n-1)^{2}}{4}}-2\sigma+\frac{1}{2}\right)^{-1},$

where $\Delta_{g_{0}}$ is the Laplace–Beltrami operator and
$[g]=\\{\bar{g}=u^{{4}/{(n-2\sigma)}}g:u\in\mathcal{C}_{+}^{\infty}(M)\\}$ is
the conformal class of $g$, where $u\in\mathcal{C}_{+}^{\infty}(M)$ if and
only if $u\in\mathcal{C}^{\infty}(M)$ and $u>0$. Furthermore, one has the
transformation law

$P_{2\sigma}({g})\phi=\mathrm{u}^{-\frac{n+2\sigma}{n-2\sigma}}P_{2\sigma}({g_{0}})(\mathrm{u}\phi)\quad\mbox{for
all}\quad\phi\in\mathcal{C}^{\infty}(\mathbb{S}^{n}\setminus\Lambda),$

which means that GJMS operators are conformally covariant. Hence, finding
conformal complete metrics $g=\mathrm{u}^{{4}/{(n-2\sigma)}}g_{0}$ with
prescribed curvature $Q_{2\sigma}(g)=Q_{n,\sigma}$ on
$\mathbb{S}^{n}\setminus\Lambda$ is equivalent to finding smooth positive
solutions $\mathrm{u}\in\mathcal{C}^{\infty}(\mathbb{S}^{n}\setminus\Lambda)$
to the nonlocal higher order geometric PDE

$\begin{cases}P_{2\sigma}({g_{0}})\mathrm{u}=c_{n,\sigma}\mathrm{u}^{\frac{n+2\sigma}{n-2\sigma}}\quad\mbox{on}\quad\mathbb{S}^{n}\setminus\Lambda,\\\
\liminf_{\mathrm{d}(p,\Lambda)\rightarrow 0}\mathrm{u}(p)=+\infty,\end{cases}$
($\mathcal{Q}_{2\sigma,\Lambda,g_{0}}$)

where $c_{n,\sigma}>0$ is a normalizing constant and ${\rm
sing}(\mathrm{u}):=\Lambda$ denotes the singular set.

Next, it will be convenient to transfer the PDE
($\mathcal{Q}_{2\sigma,\Lambda,g_{0}}$) to Euclidean space, which we can do
using the standard stereographic projection. In these coordinates, our
conformal metric takes the form
$g=\mathrm{u}^{{4}/{(n-2\sigma)}}g_{0}=(\mathrm{u}\cdot u_{\rm
sph})^{{4}/{(n-2\sigma)}}\delta$. Thus,
$u\in\mathcal{C}^{\infty}(\mathbb{R}^{n}\setminus\Sigma)$ given by
$u=\mathrm{u}\cdot u_{\rm sph}$ is a positive singular solution to
($\mathcal{Q}_{2\sigma,\Sigma}$). As a notational shorthand, we adopt the
convention that $\mathrm{u}$ refers to a conformal factor relating the metric
$g$ to the round metric, i.e. $g=\mathrm{u}^{{4}/{(n-2\sigma)}}g_{0}$, while
$u$ refers to a conformal factor relating the metric $g$ to the Euclidean
metric, i.e. $g=u^{{4}/{(n-2\sigma)}}\delta$, with the two related by
$u=\mathrm{u}\cdot u_{\rm sph}$. Hence, we aim to construct positive singular
solutions $u\in\mathcal{C}^{\infty}(\mathbb{R}^{n}\setminus\Sigma)$ to the
following higher order fractional Yamabe equation with prescribed isolated
singularities

$\begin{cases}(-\Delta)^{\sigma}u=f_{\sigma}(u)\quad{\rm
in}\quad\mathbb{R}^{n}\setminus\Sigma,\\\
u(x)=\mathcal{O}(|x|^{2\sigma-n})\quad{\rm
as}\quad|x|\rightarrow+\infty,\end{cases}$ ($\mathcal{Q}_{2\sigma,\Sigma}$)

where $\sigma\in(1,+\infty)$ with $n>2\sigma$. The subset ${\rm
sing}(u):=\Sigma\subset\mathbb{R}^{n}$ is called the singular set, which is
assumed to be $\Sigma=\\{x_{1},\cdots,x_{N}\\}$ for some $N\in\mathbb{N}$ and
such that

$\liminf_{\mathrm{d}(x,\Sigma)\rightarrow 0}u(x)=+\infty.$

We are interested in fast-decaying solutions; we assume the following
condition $\lim_{|x|\rightarrow+\infty}u(x)=0$. The integral operator on the
right-hand side of ($\mathcal{Q}_{2\sigma,\Sigma}$) is the so-called higher
order fractional Laplacian which is defined as

$(-\Delta)^{\sigma}:=(-\Delta)^{s}\circ(-\Delta)^{m},$

where $m:=[\sigma]$ and $s:=\sigma-[\sigma]$.

Here $(-\Delta)^{m}=(-\Delta)\circ\cdots\circ(-\Delta)$ denotes the poly-
Laplacian and $(-\Delta)^{s}$ denotes the fractional Laplacian defined as

$(-\Delta)^{s}u(x):={\rm
p.v.}\int_{\mathbb{R}^{n}}\mathcal{K}_{s}(x-y)[u(x)-u(y)]\mathrm{d}y,$

where
$\mathcal{K}_{s}:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a
singular potential given by

$\mathcal{K}_{s}(x-y):=\kappa_{n,s}|x-y|^{-(n+2s)}$ (2)

with

$\kappa_{n,s}=\pi^{-\frac{n}{2}}2^{2s}s{\Gamma\left(\frac{n}{2}+s\right)}{\Gamma(1-s)^{-1}}.$

The nonlinearity $f_{\sigma}:\mathbb{R}\rightarrow\mathbb{R}$ in the left-hand
side of ($\mathcal{Q}_{2\sigma,\Sigma}$) is given by

$f_{\sigma}(\xi)=c_{n,\sigma}|\xi|^{\frac{n+2\sigma}{n-2\sigma}},$

where

$c_{n,\sigma}:=2^{2\sigma}{\Gamma\left(\frac{n+2\sigma}{4}\right)^{2}}{\Gamma\left(\frac{n-2\sigma}{4}\right)^{-2}}.$

is a normalizing constant. We remark that this nonlinearity has critical
growth in the sense of the Sobolev embedding
$H^{\sigma}(\mathbb{R}^{n})\hookrightarrow
L^{2^{*}_{\sigma}}(\mathbb{R}^{n})$, where
$2^{*}_{\sigma}:=\frac{2n}{n-2\sigma}$.

Our main result in this manuscript extends this result for the remaining
cases. We are based on the unified approach given by Ao et al.[36] and Jin and
Xiong [37] to prove the existence of solutions to our integral equation, which
can be stated as follows

######  Theorem 2.

Let $\sigma\in(1,+\infty]$ with $n>2\sigma$. For any configuration
$\Sigma=\\{x_{1},\dots,x_{N}\\}$ with $N\geqslant 2$, one can find a smooth
positive singular solution to ($\mathcal{Q}_{2\sigma,\Sigma}$) such that ${\rm
sing}(u)=\Sigma$.

In [37], the authors use a dual representation and maximization methods to
study the existence of Emden–Fowler solution on the range
$\sigma\in(0,\frac{n}{2})$. Although this representation is enough to prove
the existence of blow-up limit solutions by direct maximization methods, it is
unsuitable for our gluing technique. This paper follows the approach in [36]
with the dual equation (${\mathcal{Q}}^{\prime}_{2\sigma,\Sigma}$).
Nevertheless, we need to give an alternative proof to describe the local
behavior near each isolated singularity in terms of the bubble tower solution
(see Lemma 4.10). This alternative proof is the main feature of this paper
since it enables us to extend the techniques in [36] for integral equations
that cannot be realized as the dual formulation of a differential equation,
which is undoubted of independent interest.

Instead, we notice that ($\mathcal{Q}_{2\sigma,\Sigma}$) has a dual
counterpart, which is given by

$\begin{cases}u=(-\Delta)^{-\sigma}(f_{\sigma}\circ u)\quad{\rm
in}\quad\mathbb{R}^{n}\setminus\Sigma,\\\
u(x)=\mathcal{O}(|x|^{2\sigma-n})\quad{\rm
as}\quad|x|\rightarrow+\infty,\end{cases}$
(${\mathcal{Q}}^{\prime}_{2\sigma,\Sigma}$)

where $(-\Delta)^{-\sigma}$ denotes the inverse operator of the standard
higher order fractional Laplacian, namely

$(-\Delta)^{-\sigma}f_{\sigma}(u(x)):=(\mathcal{R}_{\sigma}\ast
f_{\sigma}(u))(x)={\rm
p.v.}\int_{\mathbb{R}^{n}}\mathcal{R}_{\sigma}(x-y)f_{\sigma}(u(y))\mathrm{d}y,$

where
$\mathcal{R}_{\sigma}:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}$
is the Riesz potential given by

$\mathcal{R}_{\sigma}(x-y):=C_{n,\sigma}|x-y|^{-(n-2\sigma)}$ (3)

with $C_{n,\sigma}>0$ a normalizing constant. Our starting point in this paper
will be to prove that ($\mathcal{Q}_{2\sigma,\Sigma}$) and
(${\mathcal{Q}}^{\prime}_{2\sigma,\Sigma}$) are equivalent (see Lemma 3.3).

When $\sigma\in\mathbb{N}$ is integer, that is $\sigma=m$, Eq.
($\mathcal{Q}_{2\sigma,\Sigma}$) becomes the poly-harmonic equation

$\begin{cases}(-\Delta)^{m}u=f_{m}(u)\quad{\rm
in}\quad\mathbb{R}^{n}\setminus\Sigma,\\\
u(x)=\mathcal{O}(|x|^{2\sigma-n})\quad{\rm
as}\quad|x|\rightarrow+\infty.\end{cases}$ ($\mathcal{P}_{2m,\Sigma}$)

The most natural case of ($\mathcal{P}_{2m,\Sigma}$) is when $m=1$. In this
situation, this equation becomes the classical Lane–Emden equation. On this
subject, Mazzeo and Pacard [46, Theorem 2] based on a gluing technique via ODE
theory to prove an existence theorem. Furthermore, when $\sigma\in(0,1)$, we
arrive at

$\begin{cases}(-\Delta)^{s}u=f_{s}(u)\quad{\rm
in}\quad\mathbb{R}^{n}\setminus\Sigma,\\\
u(x)=\mathcal{O}(|x|^{2\sigma-n})\quad{\rm
as}\quad|x|\rightarrow+\infty.\end{cases}$ ($\mathcal{F}_{2s,\Sigma}$)

Recently, Ao et al. [8, Theorem 1.1] extended the earlier existence results
for this case. Their construction is substantially different from the previous
one and relies on the concept of a bubble tower (or Half-Dancer solutions). We
can summarize these results in the following statement

######  Theorem A.

Let $\sigma\in(0,1]$ with $n>2\sigma$. For any configuration
$\Sigma=\\{x_{1},\dots,x_{N}\\}$ with $N\geqslant 2$, one can find a smooth
positive singular solution to ($\mathcal{Q}_{2\sigma,\Sigma}$) such that ${\rm
sing}(u)=\Sigma$.

Let us briefly explain our strategy for the proof. We are based Schoen’s [53]
tactic, which consists in finding an explicit infinite set of functions that
span an approximate nullspace, such that the linearized nonlinear nonlocal
operator around this infinite-dimensional family of solutions is invertible on
its orthogonal complement. He first solves the equation on the complement.
Then he provides a set of balancing conditions to ensure that the solution to
this restricted problem is a solution to the original problem. This method was
recently extended for fractional operators [6].

This technique set differs substantially from the one in [46]. In their
construction, the authors obtain an one-parameter family of solutions which
blows up quickly enough near the singular set. These solutions are different
in spirit from the ones in the nonlocal case. Since blow-limit Delaunay
solutions for the scalar curvature problem are classified to depend only on
two parameters in [18]. Then, by linearizing the problem around these
solutions families, the resulting linear operator is proved to be surjective
on some reasonable space of functions, at least when the neck size parameter
is sufficiently small. A standard iteration argument may be used to obtain an
exact solution to ($\mathcal{Q}_{2\sigma,\Sigma}$) with a suitable blow-up
rate. This strategy derives from its connections with the earlier
constructions of the CMC with Delaunay-type ends [48]. Compared with the
fractional case $\sigma\in(0,1)$, the main difference in our strategy is the
proof of the refined asymptotics near half-bubble towers solutions, which
holds in a more general setting.

Using this approach, we can perturb each bubble within the tower separately
and construct a bubble tower at each singularity, and as an appropriate
approximate solution to ($\mathcal{Q}_{2\sigma,\Sigma}$). However, it is
essential to note that the linearization of this approximate solution is not
injective, as there is an infinite-dimensional kernel. As a result, an
infinite-dimensional Lyapunov-Schmidt reduction procedure is employed. This
approach is similar to Kapouleas’ CMC construction [38], which Malchiodi
adapted in [44] to produce new entire solutions for a semilinear equation with
a subcritical exponent that differ from the well-known spike solutions since
they decay to zero when moving away from three half lines and do not tend to
zero at infinity. For this, he constructed a half-Dancer solution along each
half-line. Whence, to solve the original problem from the perturbed one, an
infinite-dimensional system of Toda-type needs to be solved, which arises from
studying the interactions between the different bubbles in the tower. The most
robust interactions occur in the zero-mode level and turn into some
compatibility conditions (see Definition 5.5). In this fashion, a
configuration satisfying such conditions is called balanced, related to well-
known balancing properties enjoyed by the sum of the Pohozaev invariants.
Nonetheless, the remaining interactions can be made small and are dealt with
through a fixed-point argument.

These compatibility conditions do not restrict the location of the singularity
points but only affect the Delaunay parameter (neck size) at each end. We also
note that due to the heavy influence of the underlying geometry, the first
compatibility condition is similar to the ones found in [47] for the local
case $\sigma=1$. However, the rest of the configuration depends on the Toda-
type system. In the local setting, a similar procedure to remove the
resonances of the linearized problem was considered in [10], but the Toda-type
system was finite-dimensional in their case.

On the technical level, our strategy is to employ the gluing method and
Lyapunov–Schmidt reduction method. First, we find a suitable approximate
solution: a perturbation of the summation of half-Delaunay solutions with a
singularity at each puncture. Then, use the reduction method to find a
perturbed solution that satisfies the associated linearized problem with the
right-hand side given by some Lagrangian multiplier containing the approximate
kernels of the linearized operator. This family of kernels spans an infinite-
dimensional set called the approximate null space. The last step is to
determine the infinitely-dimensional free parameter set such that all the
coefficients of the projection on approximate null space vanish. This problem
is reduced to the solvability of some infinite-dimensional Toda system around
each singular point. A fundamental property in the proof is to have a
sufficiently good approximate solution (a half-Dancer) so that all the
estimates are exponentially decreasing in terms of the bubble tower parameter.
A fixed point argument in suitable weighted sequence spaces then solves the
problem of adjusting the parameters to have all equal to zero.

We remark that instead of relying on the well-known extension problem for the
fractional Laplacian [17], we are inspired by the approach in the approach
given by Delatorre et al. [26] to rewrite the fractional Laplacian in radial
coordinates in terms of a new integro-differential operator in logarithmic
cylindrical coordinates. In our case, such an extension does not exist in
general. We emphasize that our proof is written solely in the dual
formulation, and it can be extended to general integral equations not arising
as the dual representation of a differential equation.

In light of the seminal result of Mazzeo, Pollack, and Uhlenbeck [49] (see
also [41]), it is natural to wonder if the marked moduli space in (1) can be
furnished with more structure. It is believed that the result below should
hold

######  Conjecture 1.

Let $\sigma\in(0,+\infty]$ with $n>2\sigma$. For any singular set
$\Lambda\subset\mathbb{S}^{n}$ such that $\\#\Lambda=N$ with $N\geqslant 2$,
the marked moduli space of complete constant higher order fractional singular
$Q$-curvature metrics on the punctured round sphere
$\mathcal{M}_{2\sigma,N}(g_{0})$ is an analytic manifold with formal dimension
equal the number of isolated singularities, that is, ${\rm
dim}(\mathcal{M}_{2\sigma,N}(g_{0}))=N$.

Another possible development is to study the case in which a singular set is a
disjoint union of smooth submanifolds with possible distinct positive
Hausdorff dimensions. In this situation, it would be interesting to prove that
the moduli space defined in (1) is still non-empty and, in strong contrast
with the case of isolated singularities, is infinite-dimensional; this will be
the topic of a forthcoming paper.

Let us explain this case in more detail. It is well-known that the character
of the analysis required to prove the existence of solutions when $R(g)<0$,
which dates back to the work of Loewner and Nirenberg [43] (see also [11,
29]), is fundamentally different than in the positive scalar curvature case.
Therefore, most of the literature is concentrated on the positive scalar
curvature case $R(g)>0$. In this setting, it is natural to have a solution one
needs to impose some necessary conditions on the dimension. More challenging
it would be to construct solutions to ($\mathcal{Q}_{2\sigma,\Sigma}$) with
uncountable isolated singularities, for instance, in the lattice
$\Sigma=\mathbb{Z}^{n}$. The existence of weak solutions with larger dimension
singular set for the singular Yamabe equation has been studied by Mazzeo and
Smale [50] and by Mazzeo and Pacard [46] for the scalar curvature case. As
well as by Hyder and Sire [35] for the (fourth order) $Q$-curvature metrics,
and by Ao et al. [6] for the (fractional order) $Q$-curvature metrics, based
in the construction of entire solutions from [7].

More generally, such solutions may be constructed on an arbitrary compact
manifold $(M^{n},g)$ of nonnegative scalar curvature $R(g)\geqslant 0$
whenever the singular set is a finite disjoint union of submanifolds with
positive bounded Hausdorff dimension, which we describe as follows. Given
$\sigma\in(0,+\infty]$ with $n>2\sigma$ and $N\geqslant 2$, we let
$\Lambda\subset\mathbb{S}^{n}$ be a finite disjoint union of submanifolds
$\Lambda=\Lambda^{0}\cup\Lambda^{+}$, where
$\Lambda_{+}=\cup_{\ell=1}^{d}\Lambda_{\ell}^{k_{\ell}}$ with $k_{\ell}:={\rm
dim}_{\mathcal{H}}(\Lambda_{\ell})$ denoting its Hausdorff dimension.
Furthermore, we denote by

$\displaystyle\mathcal{M}^{k}_{2\sigma,\Lambda}(g_{0})=\left\\{g\in[g_{0}]:\mbox{$g$
is complete on
$\mathbb{S}^{n}\setminus\Lambda^{k}$}\;\mbox{and}\;Q_{2\sigma}(g)\equiv
Q_{n,\sigma}\right\\}$

the moduli space of complete constant higher order fractional $Q$-curvature
metrics with higher dimensional singularities. Notice that we simply denote
$\mathcal{M}^{0}_{2\sigma,\Lambda}(g_{0})=\mathcal{M}_{2\sigma,\Lambda}(g_{0})$.

To summarize this discussion, we have the following statement

######  Theorem B.

Let $\sigma\in(0,+\infty]$ with $n>2\sigma$. Assume that
$\Lambda=\Lambda^{0}\cup\Lambda^{+}$ is a finite disjoint union of
submanifolds satisfying $\Lambda^{0}=\varnothing$ and
$\Lambda_{+}=\cup_{\ell=1}^{d}\Lambda_{\ell}^{k_{\ell}}$ with
$0<k_{\ell}<\frac{n-2\sigma}{2}$. Then, there exists a metric
$g\in\mathcal{M}^{k}_{2\sigma,\Lambda}(g_{0})$ that ${\rm sing}(g)=\Lambda$.
In particular, one has
$\mathcal{M}^{k}_{2\sigma,\Lambda}(g_{0})\neq\varnothing$ and it is an
infinite-dimensional analytic manifold.

With a mind on this statement, it would be natural to prove a similar result
as below

######  Conjecture 2.

Let $\sigma\in(0,+\infty]$ with $n>2\sigma$. Assume that
$\Lambda=\Lambda^{0}\cup\Lambda^{+}$ such that $\\#\Lambda^{0}=N$ with
$N\geqslant 2$ and $\Lambda_{+}=\cup_{\ell=1}^{d}\Lambda_{\ell}^{k_{\ell}}$
with $0<k_{\ell}<\frac{n-2\sigma}{2}$. Then, there exists a metric
$g\in\mathcal{M}^{k}_{2\sigma,\Lambda}(g_{0})$ satisfying that ${\rm
sing}(g)=\Lambda$. In particular, one has
$\mathcal{M}^{k}_{2\sigma,\Lambda}(g_{0})\neq\varnothing$ and it is an
infinite-dimensional analytic manifold.

As usual, we need to prove the existence of positive solutions to the GJMS
equation on the conformally flat case
$\mathbb{S}^{n}\setminus\mathbb{S}^{k}\simeq\mathbb{R}^{n}\setminus\mathbb{R}^{k}$
with fast-decay away from the singular set. Moreover, the dimension estimate
above is sharp in the same sense of González, Mazzeo, and Sire [31]. Namely,
if a complete metric blows up at a smooth $k$-dimensional submanifold and is
polyhomogeneous, then $k\in\mathbb{R}_{+}$ must satisfy the restriction above.
All the analysis for this type of equation comes from its conformal
properties, which produce a geometric interpretation of scattering theory and
conformally covariant operators. Exploiting the conformal equivalence
$\mathbb{R}^{n}\setminus\mathbb{R}^{k}\simeq\mathbb{S}^{n-k-1}\times\mathbb{H}^{k+1}$,
where $\mathbb{R}_{+}^{n+1}$ is replaced by anti-de Sitter (AdS) space, but
the arguments run in parallel. In the same direction but with another flavor,
we quote the multiplicity results in [13, 14, 12, 4, 21], which also exploit
this conformal invariance and are based on a topological bifurcation
technique; this is believed to be true in the much broader case of conformally
variational invariants (c.f. [20, 19]).

One could extend this construction in a more geometric direction for not-round
metrics. On this subject, we cite [52, 15] for non-flat gluing constructions
for the constant curvature equation. Recently, in [2], a similar gluing
construction is used to prove existence results for fourth order constant
$Q$-curvature nondegenerate metrics with suitable growth conditions on the
Weyl tensor.

We now describe the plan for the rest of the paper. In Section 2, we establish
some terminology that will be used throughout the paper. In Section 3, we
prove the dual representation formula relating
($\mathcal{Q}_{2\sigma,\Sigma}$) with
(${\mathcal{Q}}^{\prime}_{2\sigma,\Sigma}$). In Section 4, we classify
Delaunay-type solutions as bubble towers. In Section 5, we provide balancing
equations. Next, we define balanced configuration parameters and admissible
perturbation sequences. We use this to define approximate solutions and prove
some estimates for the linearized operator around this approximating family.
In Section 6, we summarize some estimates involving the coefficients of the
projection on the approximate null space. In Section 7, we reduce the proof of
Theorem 2 to solving an infinite dimensional Toda system. We prove that under
admissibility conditions, this system can be solved. In Appendix A, we recall
some estimates concerning the interaction between two spherical solutions with
different centers and radii.

##  2\. Notations

We establish some notations and definitions that we will use frequently
throughout the text for easy reference.

  * •

$m:=\lfloor\sigma\rfloor$ is the integer part of $\sigma$, that is, be the
greatest integer that does not exceed $\sigma$;

  * •

$s:=\\{\sigma\\}$ is the fractional part of $\sigma$, that is,
$s:=\sigma-\lfloor\sigma\rfloor$;

  * •

$0<\xi,\nu,\zeta_{1}\ll 1$ are small constants;

  * •

$C>0$ is a universal constant that may vary from line to line and even in the
same line.

  * •

$a_{1}\lesssim a_{2}$ if $a_{1}\leqslant Ca_{2}$, $a_{1}\gtrsim a_{2}$ if
$a_{1}\geqslant Ca_{2}$, and $a_{1}\simeq a_{2}$ if $a_{1}\lesssim a_{2}$ and
$a_{1}\gtrsim a_{2}$.

  * •

$u=\mathcal{O}(f)$ as $x\rightarrow x_{0}$ for
$x_{0}\in\mathbb{R}\cup\\{\pm\infty\\}$, if $\limsup_{x\rightarrow
x_{0}}(u/f)(x)<\infty$ is the Big-O notation;

  * •

$u=\mathrm{o}(f)$ as $x\rightarrow x_{0}$ for
$x_{0}\in\mathbb{R}\cup\\{\pm\infty\\}$, if $\lim_{x\rightarrow
x_{0}}(u/f)(x)=0$ is the little-o notation;

  * •

$u\simeq\widetilde{u}$, if $u=\mathcal{O}(\widetilde{u})$ and
$\widetilde{u}=\mathcal{O}(u)$ as $x\rightarrow x_{0}$ for
$x_{0}\in\mathbb{R}\cup\\{\pm\infty\\}$;

  * •

$\mathcal{C}^{j,\alpha}(\mathbb{R}^{n})$, where $j\in\mathbb{N}$ and
$\alpha\in(0,1)$, is the Hölder space; we simply denote
$\mathcal{C}^{j}(\mathbb{R}^{n})$ if $\alpha=0$.

  * •

$W^{m,q}(\mathbb{R}^{n})$ is the classical Sobolev space, where
$m\in\mathbb{N}$ and $q\in[1,+\infty]$; when $m=0$ we simply denote
$L^{q}(\mathbb{R}^{n})$ when $q=2$, we simply denote
$H^{\ell}(\mathbb{R}^{n})$;

  * •

$\mathcal{C}^{2\sigma}(\mathbb{R}^{n})={\mathcal{C}^{2m,2s}}(\mathbb{R}^{n})$
is the classical Hölder space;

  * •

$\gamma_{\sigma}=\frac{n-2\sigma}{2}$ is the Fowler rescaling exponent with
$\gamma_{\sigma}^{\prime}=\frac{n+2\sigma}{2}$ its dual;

  * •

$2^{*}_{\sigma}=\frac{2n}{n-2\sigma}$ is the critical Sobolev exponent with
$2^{*\prime}_{\sigma}=\frac{2n}{n+2\sigma}$ its dual;

  * •

$A_{1},A_{2}>0,A_{3}<0$ are constant defined by (A.1), (A.2), and (A.3),
respectively;

  * •

$\mathcal{I}_{\infty}:=\\{1,\dots,N\\}\times\mathbb{N}\times\\{0,\dots,n\\}\simeq\ell^{\infty}(\mathbb{R}^{(n+1)N})$
is total index set;

  * •

$\boldsymbol{p}:=(p_{1},\dots,p_{N})\in\mathbb{S}^{nN}$ with
$\Lambda:=\\{p_{1},\dots,p_{N}\\}\subset\mathbb{S}^{n}$;

  * •

$\boldsymbol{x}:=(x_{1},\dots,x_{N})\in\mathbb{R}^{nN}$ with
$\Sigma:=\\{x_{1},\dots,x_{N}\\}\subset\mathbb{R}^{n}$;

  * •

$\boldsymbol{L}=(L_{1},\dots,L_{N})\in\mathbb{R}^{N}_{+}$ is a vector of
periods such that $|\boldsymbol{L}|\gg 1$ is large enough arising from
Proposition 5.3. Equivalently,
$\boldsymbol{\varepsilon}=(\varepsilon_{1},\dots,\varepsilon_{N})\in\mathbb{R}^{N}_{+}$
is a vector of necksizes such that $0<|\boldsymbol{\varepsilon}|\ll 1$ is
small enough;

  * •

$(\boldsymbol{x},\boldsymbol{L})\in\mathbb{R}^{(n+1)N}$ are the moduli space
parameters and $\Upsilon_{\rm
conf}:\mathbb{R}^{(n+1)N}\rightarrow\mathbb{R}^{(n+2)N}$ is the configuration
map;

  * •

$\boldsymbol{q}=(q_{1},\dots,q_{N})\in\mathbb{R}^{N}_{+}$ is a vector of
comparable periods such that $|\boldsymbol{q}|\gg 1$ is also large enough and
satisfy (5.10), $\boldsymbol{R}=(R^{1},\dots,R^{N})\in\mathbb{R}^{N}$ and
$\boldsymbol{a}_{0}=(a_{0}^{1},\dots,a_{0}^{N})\in\mathbb{R}^{nN}$ are the
deformation parameters;

  * •

$(\boldsymbol{q},\boldsymbol{a}_{0},\boldsymbol{R})\in\mathbb{R}^{(n+2)N}$ are
the configurations parameters and $\Upsilon_{\rm
per}:\mathbb{R}^{(n+2)N}\rightarrow\ell^{\infty}_{\tau}(\mathbb{R}^{(n+1)N})$
is the configuration map;

  * •

$(\boldsymbol{q}^{b},\boldsymbol{a}_{0}^{b},\boldsymbol{R}^{b})\in{\rm
Bal}_{\sigma}(\Sigma)$ denotes a balanced configuration, that is, it satisfies
($\mathscr{B}_{1}$) and ($\mathscr{B}_{2}$).

  * •

$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\in\ell^{\infty}_{\tau}(\mathbb{R}^{(n+1)N})$
(or
$(\boldsymbol{a}_{j},\boldsymbol{r}_{j})\in\ell^{\infty}_{\tau}(\mathbb{R}^{(n+1)N})$)
are the perturbation sequences and $\Upsilon_{\rm
per}:\mathbb{R}^{(n+2)N}\rightarrow\ell^{\infty}_{\tau}(\mathbb{R}^{(n+1)N})$
is the perturbation map;

  * •

$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\in{\rm Adm}_{\sigma}(\Sigma)$
denotes the admissible perturbation sequences, that is, it satisfies
($\mathscr{A}_{0}$) and ($\mathscr{A}_{1}$); equivalently $\Upsilon_{\rm
per}^{-1}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\in{\rm
Bal}_{\sigma}(\Sigma)$;

  * •

$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}\in\mathcal{C}_{*,\tau}(\mathbb{R}^{n}\setminus\Sigma)$
denotes a Delaunay solution with associated error denoted by
$\phi_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}\in\mathcal{C}_{*,\tau}(\mathbb{R}^{n}\setminus\Sigma)$
and $\Upsilon_{\rm
sol}:\ell^{\infty}_{\tau}(\mathbb{R}^{(n+1)N})\rightarrow\mathcal{C}_{*,\tau}(\mathbb{R}^{n}\setminus\Sigma)$
is the solution map;

  * •

$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}\in{\rm
Apx}_{\sigma}(\Sigma)$ is an approximate solution, that is, $\Upsilon_{\rm
sol}^{-1}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})})\in{\rm
Adm}_{\sigma}(\Sigma)$;

  * •

$\\{Z_{j\ell}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\\}_{(i,j,\ell)\in\mathcal{I}_{\infty}}\subset\mathcal{C}^{0}(\mathbb{R}^{n}\setminus\Sigma)$
is the associated family of approximating normalized cokernels;

  * •

$\\{c_{j\ell}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\\}_{(i,j,\ell)\in\mathcal{I}_{\infty}}\subset\mathcal{C}^{0}(\mathbb{R}^{n}\setminus\Sigma)$
is the associated family of coefficients of the projections on approximating
normalized cokernels;

  * •

$\\{\beta_{j\ell}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\\}_{(i,j,\ell)\in\mathcal{I}_{\infty}}\subset\mathcal{C}^{0}(\mathbb{R}^{n}\setminus\Sigma)$
is the associated family of projections on approximating normalized cokernels.

##  3\. Dual representation formula

This section shows that our equation and its dual are correspondents. We are
based on the removable singularity result from [9, Theorem 1.1]. We also refer
to [51, Proposition 4.1] for the integer cases $\sigma\in\mathbb{N}$. In what
follows, we consider the space

$L_{s}(\mathbb{R}^{n}):=\left\\{u\in L^{1}_{\rm
loc}(\mathbb{R}^{n}):\int_{\mathbb{R}^{n}}\frac{|u(x)|}{1+|x|^{n+2s}}\mathrm{d}x<+\infty\right\\}$

with $s\in(0,1)$.

We first introduce the notation of distributional solutions to
($\mathcal{Q}_{2\sigma,\Sigma}$).

######  Definition 3.1.

Let $\sigma\in\mathbb{R}_{+}$ and $n>2\sigma$. We say that a smooth solution
$u\in\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\Sigma)\cap L^{1}_{\rm
loc}(\mathbb{R}^{n})$ to ($\mathcal{Q}_{2\sigma,\Sigma}$) is a solution in the
distributional sense to ($\mathcal{Q}_{2\sigma,\Sigma}$) if the equality below
holds

$\int_{\mathbb{R}^{n}}u(-\Delta)^{\sigma}\varphi\mathrm{d}x=\int_{\mathbb{R}^{n}}f_{\sigma}(u)\varphi\mathrm{d}x\quad{\rm
in}\quad\mathbb{R}^{n}\setminus\Sigma$ (3.1)

for all $\varphi\in\mathcal{C}_{c}^{\infty}(\mathbb{R}^{n}\setminus\Sigma)$.

######  Remark 3.2.

One can check that smooth solution to ($\mathcal{Q}_{2\sigma,\Sigma}$) are
indeed distributional solutions.

We need the following auxiliary result to prove the equivalence: a combination
of [9, Theorem 1.1 and Lemma 5.4].

######  Lemma A.

Let $\sigma\in\mathbb{R}_{+}$ and $n>2\sigma$. If
$u\in\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\Sigma)\cap L^{1}_{\rm
loc}(\mathbb{R}^{n})$ is a distributional solution to
($\mathcal{Q}_{2\sigma,\Sigma}$), then $f_{\sigma}\circ u\in L^{1}_{\rm
loc}(\mathbb{R}^{n})$ and $u\in L^{1}_{\rm loc}(\mathbb{R}^{n})$ is a
distributional solution in $\mathbb{R}^{n}$, that is, the distributional
equation (3.1) holds. Moreover, one has

$\int_{\mathbb{R}^{n}}\frac{f_{\sigma}(u(x))}{1+|x|^{n-2\sigma}}\mathrm{d}x<+\infty.$
(3.2)

Consequently, we obtain that
$w\in\mathcal{C}^{\infty}(\mathbb{R}^{n}\setminus\Sigma)$ is defined as

$w(x):=\int_{\mathbb{R}^{n}}\mathcal{R}_{\sigma}(x-y)f_{\sigma}(u(y))\mathrm{d}y$
(3.3)

is well-defined and belongs to $L_{s}(\mathbb{R}^{n})$ for every $s>0$.

Finally, we also recall the Liouville theorem from [34, Lemma 2.4].

######  Lemma B.

Let $\sigma\in\mathbb{R}_{+}$ and $n>2\sigma$. Assume that $w\in
L_{s}(\mathbb{R}^{n})$ for some $s\geqslant 0$ and

$(-\Delta)^{\sigma}w=0\quad{\rm in}\quad\mathbb{R}^{n},$

for some $\sigma\geqslant s$. Then, one has that $w$ is a polynomial of degree
at most $\lfloor 2s\rfloor$.

With the lemmas above, we have our main result in this section.

######  Proposition 3.3.

Let $\sigma\in\mathbb{R}_{+}$ and $n>2\sigma$. It holds that
($\mathcal{Q}_{2\sigma,\Sigma}$) and
(${\mathcal{Q}}^{\prime}_{2\sigma,\Sigma}$) are equivalents.

###### Proof of Proposition 3.3.

Let $u\in\mathcal{C}^{\infty}(\mathbb{R}^{n}\setminus\Sigma)$ be a positive
singular fast-decaying solution to ($\mathcal{Q}_{2\sigma,\Sigma}$). From
(3.2), we have that $w\in L_{s}(\mathbb{R}^{n})$ for every $s>0$, $s\neq
2\sigma$. Hence, if we define $\widehat{w}=u-w$, then $\widehat{w}\in
L_{s}(\mathbb{R}^{n})$ for all $s>0$ with $s\neq 2\sigma$. In addition, since
$(-\Delta)^{\sigma}\widehat{w}=0$ in $\mathbb{R}^{n}$, we conclude that
$\widehat{w}$ is a polynomial of degree at most $2m$, thanks to the Liouville
theorem in Lemma B. Recall that we are considering solutions satisfying
$\lim_{|x|\rightarrow+\infty}u(x)=0$. Consequently, $\widehat{w}\equiv 0$, and
the dual representation holds. ∎

##  4\. Delaunay-type solutions

This section is devoted to the construction of solutions for the case of a
single isolated singularity, that is, $\Sigma=\\{0\\}$. We are inspired in
[37], which is an adaption of the earlier constructions in [8, 26, 27] for the
cases $\sigma\in(0,1)$ and $\sigma=1$.

###  4.1. Integral Emden–Fowler coordinates

As a matter of fact, when $\Sigma=\\{0\\}$, Eq.
($\mathcal{Q}_{2\sigma,\Sigma}$) can be rewritten as

$\begin{cases}(-\Delta)^{\sigma}u=f_{\sigma}(u)\quad{\rm
in}\quad\mathbb{R}^{n}\setminus\\{0\\},\\\
\displaystyle\lim_{|x|\rightarrow+\infty}u(x)=0,\end{cases}$
($\mathcal{Q}_{2\sigma,\infty}$)

or into its dual form

$\begin{cases}u=(-\Delta)^{-\sigma}(f_{\sigma}\circ u)\quad{\rm
in}\quad\mathbb{R}^{n}\setminus\\{0\\}\\\
\displaystyle\lim_{|x|\rightarrow+\infty}u(x)=0.\end{cases}$
($\mathcal{Q}^{\prime}_{2\sigma,\infty}$)

It is straightforward to see from Proposition 3.3 that
($\mathcal{Q}_{2\sigma,\infty}$) are ($\mathcal{Q}^{\prime}_{2\sigma,\infty}$)
equivalents.

######  Remark 4.1.

For any $\sigma\in(1,+\infty]$ and $n>2\sigma$, there are two distinguished
solutions to ($\mathcal{Q}^{\prime}_{2\sigma,\infty}$), which we describe as
follows:

  * (a)

The cylindrical solution

$u_{\rm cyl}(|x|)=a_{n,\sigma}|x|^{-\gamma_{\sigma}},$ (4.1)

which is singular at the origin.

  * (b)

The standard spherical solution $($also known as “bubble” solution$)$

$u_{\rm
sph}(|x|)=\left(\frac{2}{1+\left|x\right|^{2}}\right)^{\gamma_{\sigma}},$
(4.2)

which is non-singular at the origin.

We remark that all non-singular solutions to the blow-up limit problem were
classified in [24], which are given by deformations of the standard bubble
solution. This reflects the invariance of equation
($\mathcal{Q}^{\prime}_{2\sigma,\infty}$) with respect to the entire Euclidean
group with translations and dilations.

######  Proposition A.

Let $\sigma\in(1,+\infty]$ and $n>2\sigma$. If
$u\in\mathcal{C}^{2\sigma}(\mathbb{R}^{n})$ is a positive smooth non-singular
solution to ($\mathcal{Q}^{\prime}_{2\sigma,\infty}$), then there exists
$\lambda\in\mathbb{R}$ and $x_{0}\in\mathbb{R}^{n}$ such that

$u\equiv U_{\lambda,x_{0}},$ (4.3)

where

$U_{\lambda,x_{0}}(x)=\left(\frac{2\lambda}{\lambda^{2}+\left|x-x_{0}\right|^{2}}\right)^{\gamma_{\sigma}}$
(4.4)

for some $\lambda>0$ and $x_{0}\in\mathbb{R}^{n}$. This family of solutions
will be called spherical or bubble solutions.

The problem of classifying the complete set of positive smooth singular
solutions to ($\mathcal{Q}^{\prime}_{2\sigma,\infty}$) is much more
challenging and only accomplished for a few cases. On this subject, Chen, Li,
and Ou proved that all solutions are radially symmetric with respect to the
origin. In addition, Jin and Xiong [37] only proved the existence of such a
solution by a direct maximization method. Furthermore, they also study the
local asymptotic behavior of positive singular solutions to

$(-\Delta)^{\sigma}u=f_{\sigma}(u)\quad{\rm in}\quad B_{R}^{*},$
($\mathcal{Q}_{2\sigma,R}$)

or into its dual form

$u=(-\Delta)^{-\sigma}(f_{\sigma}\circ u)\quad{\rm in}\quad B_{R}^{*},$
($\mathcal{Q}^{\prime}_{2\sigma,R}$)

where $B_{R}^{*}\subset\mathbb{R}^{n}\setminus\\{0\\}$ given by
$B_{R}^{*}=B_{R}(0)\setminus\\{0\\}$ is the punctured ball of radius $R>0$.

To study this class of equations, we define an important change of variables
that turns ($\mathcal{Q}_{2\sigma,\infty}$) into an integral one-dimensional
problem.

######  Definition 4.2.

Let $\sigma\in(1,+\infty]$ and $n>2\sigma$ We define the integral Emden–Fowler
change of variables $($or cylindrical logarithm coordinates$)$ given by

$\mathfrak{F}_{\sigma}:\mathcal{C}_{c}^{\infty}(B_{R}^{*})\rightarrow\mathcal{C}_{c}^{\infty}(\mathcal{C}_{L})\quad\mbox{given
by}\quad\mathfrak{F}_{\sigma}(u)=e^{-\gamma_{\sigma}t}u(e^{-t},\theta),$ (4.5)

where $t=-\ln R$, $\theta=x/|x|$, $\mathcal{C}_{L}:=(L,+\infty)$ with
$L=-\ln|x|$ and $\gamma_{\sigma}:=\frac{n-2\sigma}{2}$. The inverse of this
isomorphism is

$\left(\mathfrak{F}_{\sigma}\right)^{-1}:\mathcal{C}_{c}^{\infty}(\mathcal{C}_{L})\rightarrow\mathcal{C}_{c}^{\infty}(B_{R}^{*})\quad\mbox{given
by}\quad\left(\mathfrak{F}_{\sigma}\right)^{-1}(v)=|x|^{\gamma_{\sigma}}v(-\ln|x|,\theta).$
(4.6)

The quantity $\gamma_{\sigma}>0$ will be referred to as the Fowler rescaling
exponent. From now on, we denote by $v(t,\theta):=\mathfrak{F}_{\sigma}(u(x))$
and $u(x):=(\mathfrak{F}_{\sigma})^{-1}(v(t,\theta))$, conversely.

Using this change of variable, Eq. ($\mathcal{Q}_{2\sigma,R}$) can be
reformulated as the following one-dimensional problem

$\displaystyle\begin{cases}(-\Delta)_{\rm
cyl}^{\sigma}v=f_{\sigma}(v)\quad{\rm in}\quad\mathcal{C}_{L},\\\
\displaystyle\lim_{t\rightarrow+\infty}v(t)=0.\end{cases}$
($\mathcal{O}_{2\sigma,L}$)

Here $(-\Delta)_{\rm
cyl}^{\sigma}:\mathcal{C}^{2\sigma}(\mathcal{C}_{L})\rightarrow\mathcal{C}^{0}(\mathcal{C}_{L})$
is the operator higher-order operator given by

$(-\Delta)^{\sigma}_{\rm cyl}:=(-\Delta)^{s}_{\rm cyl}\circ(-\Delta)^{m}_{\rm
cyl},$ (4.7)

where $(-\Delta)^{m}_{\rm cyl}$ and $(-\Delta)^{s}_{\rm cyl}$ denote the
cylindrical poly-Laplacian and the fractional Laplacian, respectively, defined
as

$(-\Delta)^{m}_{\rm
cyl}:=\sum_{\ell=0}^{2m}\sum_{j=0}^{2m}K^{(\ell)}_{2m,j}\partial_{t}^{(j)}(-\Delta_{\theta})^{\ell},$

where $K^{(\ell)}_{2m,j}=K^{(\ell)}_{2m,j}(n)>0$ for
$j,\ell\in\\{0,\dots,2m\\}$ are dimensional constants, and

$(-\Delta)_{\rm
cyl}^{s}v(t,\theta):=\int_{-L}^{+L}\widehat{\mathcal{K}}_{\sigma}(t-\tau,\theta-\varsigma)[v(t,\theta)-v(\tau,\varsigma)]\mathrm{d}\tau\mathrm{d}\varsigma,$

where ${\mathcal{K}}_{\sigma,{\rm
cyl}}:\mathcal{C}_{L}\times\mathcal{C}_{L}\rightarrow\mathbb{R}$ is the kernel
(2) written in Emden–Fowler coordinates. As usual, the dual form of this
equation is given by

$\displaystyle\begin{cases}v=(-\Delta)_{\rm cyl}^{-\sigma}(f_{\sigma}\circ
v)\quad{\rm in}\quad\mathcal{C}_{L},\\\
\displaystyle\lim_{t\rightarrow+\infty}v(t)=0.\end{cases}$
($\mathcal{O}_{2\sigma,L}^{\prime}$)

Here $(-\Delta)_{\rm cyl}^{-\sigma}$ is the integral linear operator defined
by

$\displaystyle(-\Delta)_{\rm cyl}^{-\sigma}(f_{\sigma}\circ
v)(t,\theta):=(\widehat{\mathcal{R}}_{\sigma}\ast(f_{\sigma}\circ
v))(t,\theta)=\int_{-\infty}^{+\infty}\widehat{\mathcal{R}}_{\sigma}(t-\tau,\theta-\varsigma)f_{\sigma}(v(\tau,\varsigma))\mathrm{d}\tau,$

where
$\widehat{\mathcal{R}}_{\sigma}:\mathcal{C}_{L}\times\mathcal{C}_{L}\rightarrow\mathbb{R}$
is the Riesz kernel (3) written in Emden–Fowler coordinates. Henceforth, we
keep the notation ${\mathcal{K}}_{\sigma,{\rm
cyl}}=\widehat{\mathcal{K}}_{\sigma}$ and ${\mathcal{R}}_{\sigma,{\rm
cyl}}=\widehat{\mathcal{R}}_{\sigma}$ for the sake of simplicity.

######  Remark 4.3.

Notice that $(-\Delta)_{\rm cyl}^{-\sigma}$ is an abuse of notation, which we
keep for simplicity. In the geometric language, this change of variables
corresponds to a restriction of the conformal diffeomorphism between the
entire cylinder and the punctured space. In other words, one has

$(-\Delta)_{\rm cyl}^{-\sigma}=P_{2\sigma}(g_{\rm cyl}),$

where $g_{{\rm cyl}}=\mathrm{d}t^{2}+\mathrm{d}\theta^{2}$ stands for the
cylindrical metric and ${\mathrm{d}\theta}=e^{-2t}\delta$, where $\delta$ is
the standard flat metric.

Notice that since in the blow-up limit situation ($R=+\infty$), solutions to
($\mathcal{Q}_{2\sigma,\infty}$) are rotationally invariant, that is,
$u(x)=u(r)$ with $r=|x|$. Using this change of variable, Eq.
($\mathcal{Q}_{2\sigma,\infty}$) can be reformulated as the following one-
dimensional problem

$\displaystyle\begin{cases}(-\Delta)_{\rm
cyl}^{\sigma}v=f_{\sigma}(v)\quad{\rm in}\quad\mathbb{R},\\\
\displaystyle\lim_{t\rightarrow+\infty}v(t)=0.\end{cases}$
($\mathcal{O}_{2\sigma,\infty}$)

Here $(-\Delta)_{\rm cyl}^{\sigma}$ represents the operator higher-order
operator (written in Emden–Fowler coordinates (4.5)), namely

$(-\Delta)_{\rm
cyl}^{\sigma}v(t):=\int_{-\infty}^{+\infty}\widehat{\mathcal{K}}_{\sigma}(t-\tau)[v(t)-v(\tau)]\mathrm{d}\tau,$
(4.8)

where
$\widehat{\mathcal{K}}_{\sigma}:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$
is a kernel given by

$\displaystyle\widehat{\mathcal{K}}_{\sigma}(t)=2^{-\gamma^{\prime}_{\sigma}}\int_{\mathbb{S}^{n-1}}|\cosh(t)-\langle\theta,\tau\rangle|^{-\gamma_{\sigma}^{\prime}}\mathrm{d}\tau=\int_{\mathbb{S}^{n-1}}e^{-\gamma_{\sigma}^{\prime}t}\left(1+e^{-2t}-2e^{-t}\langle\theta,\tau\rangle\right)^{-\gamma_{\sigma}^{\prime}}\mathrm{d}\tau.$
(4.9)

As before, the dual form of this equation is given by

$\displaystyle\begin{cases}v=(-\Delta)_{\rm cyl}^{-\sigma}(f_{\sigma}\circ
v)\quad{\rm in}\quad\mathbb{R},\\\
\displaystyle\lim_{t\rightarrow+\infty}v(t)=0.\end{cases}$
($\mathcal{O}_{2\sigma,\infty}^{\prime}$)

Here $(-\Delta)_{\rm cyl}^{-\sigma}$ is the integral linear operator defined
by

$\displaystyle(-\Delta)_{\rm cyl}^{-\sigma}(f_{\sigma}\circ
v)(t):=(\widehat{\mathcal{R}}_{\sigma}\ast(f_{\sigma}\circ
v))(t)=\int_{-\infty}^{+\infty}\widehat{\mathcal{R}}_{\sigma}(t-\tau)f_{\sigma}(v(\tau))\mathrm{d}\tau,$

where
$\widehat{\mathcal{R}}_{\sigma}:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$
is a kernel given by

$\displaystyle\widehat{\mathcal{R}}_{\sigma}(t)={2^{-\gamma_{\sigma}}}\omega_{n-2}\int_{-1}^{1}{\left(1-\zeta_{1}^{2}\right)^{\frac{n-3}{2}}}{\left|\cosh(t)-\zeta_{1}\right|^{-\gamma_{\sigma}}}\mathrm{d}\zeta_{1}.$
(4.10)

######  Remark 4.4.

It is possible to express this kernel in terms of hypergeometric functions. We
also observe $\widehat{\mathcal{R}}_{\sigma}(\xi)\sim 1$ is bounded and Hölder
continuous, whereas $\widehat{\mathcal{K}}_{\sigma}(\xi)\sim|\xi|^{1-2s}$ when
$\sigma\in(1,+\infty)$. Furthermore, they behave qualitatively as

$\widehat{\mathcal{K}}_{\sigma}(\xi)\sim
e^{-\gamma_{\sigma}^{\prime}|\xi|}\quad{\rm as}\quad|\xi|\rightarrow+\infty$
(4.11)

and

$\widehat{\mathcal{R}}_{\sigma}(\xi)\sim e^{-\gamma_{\sigma}|\xi|}\quad{\rm
as}\quad|\xi|\rightarrow+\infty,$ (4.12)

where $\xi:=|t-\tau|$. We refer to [8, 37] for proof of these facts.

Using this new formulation, one has the following

######  Remark 4.5.

As before, there are two distinguished solutions to
($\mathcal{O}_{2\sigma,L}^{\prime}$), which we describe as follows:

  * (a)

The cylindrical solution, which is

$v_{\rm cyl}(t)\equiv a_{n,\sigma},$

where $v_{\rm cyl}=\mathfrak{F}_{\sigma}(u_{\rm
cyl})\in\mathcal{C}^{2\sigma}(\mathbb{R})$ with $u_{\rm
cyl}\in\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\\{0\\})$ given by (4.1).

  * (b)

The standard spherical solution $($also known as “bubble”$)$ which is

$v_{\rm sph}(t)=\cosh(t)^{\gamma_{\sigma}},$ (4.13)

where $v_{\rm sph}=\mathfrak{F}_{\sigma}(u_{\rm
sph})\in\mathcal{C}^{2\sigma}(\mathbb{R})$ with $u_{\rm
sph}\in\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\\{0\\})$ given by (4.2).

###  4.2. Asymptotic classification of Delaunay-type solutions

Now we prove the existence of even solutions to
($\mathcal{O}_{2\sigma,L}^{\prime}$) with large periods which are close to the
standard bubble tower solution given by (5.5) in a suitable weighted Hölder
norm.

First, for the standard bubble solution, we have the following nondegeneracy
result, which is based on [42, Lemma 5.1] and [25, Lemma 5.1]. In our
situation, this is proved in [39, Lemma A.1]. Nevertheless, we include a
sketch of the proof in Appendix B for completeness.

######  Lemma 4.6.

Let $\sigma\in(1,\infty)$ and $n>2\sigma$. The standard bubble solution
$u_{\rm sph}\in\mathcal{C}^{2\sigma}({\mathbb{R}^{n}})$ given by (4.2)
satisfying ($\mathcal{Q}_{2\sigma,\infty}$) is nondegenerate in a sense, the
set of bounded solutions to the linearized equation

$\phi-(-\Delta)^{-\sigma}(f^{\prime}_{\sigma}(u_{\rm sph})\phi)=0\quad{\rm
in}\quad\mathbb{R}^{n}$ (4.14)

are spanned by the functions

$\gamma_{\sigma}u_{\rm sph}+x\cdot\nabla u_{\rm sph}\quad\text{ and
}\quad\partial_{x_{i}}u_{\rm sph}\quad{\rm for}\quad i\in\\{1,\dots,n\\}.$

###### Proof.

See Appendix B. ∎

One can also reformulate the last result as follows

######  Lemma 4.7.

Let $\sigma\in(1,\infty)$ and $n>2\sigma$. The standard bubble solution
$v_{\rm sph}\in\mathcal{C}^{2\sigma}({\mathbb{R}})$ given by (4.13) satisfying
($\mathcal{O}_{2\sigma,L}^{\prime}$) is nondegenerate in the sense that all
bounded solutions of the linearized equation

$\psi-(-\Delta)_{\rm cyl}^{-\sigma}(f^{\prime}_{\sigma}(v_{\rm
sph})\psi)=0\quad{\rm in}\quad\mathbb{R}$

are spanned by the translations $v_{\rm sph}(\cdot-T)$ with $T>0$.

###### Proof.

It follows by undoing the Emden–Fowler change of variables in (4.5). ∎

Second, we restrict ourselves to the open interval $(-L,L)$ equipped with
Dirichlet boundary conditions. In what follows, we fix $L\in\mathbb{N}$. Let
$j\in\mathbb{N}$ and $\alpha\in(0,1)$, we denote by
$\mathcal{C}^{j,\alpha}_{L}(\mathbb{R})$ the classical Hölder space
$\mathcal{C}^{j,\alpha}(\mathbb{R})$ restricted to $2L$-periodic functions on
the open interval $(-L,L)$. For $\alpha=0$, we simply denote
$\mathcal{C}^{j}_{L}(\mathbb{R})$. Let $\ell\in\mathbb{N}$ and
$q\in[1,+\infty]$, we will keep the notation $W^{\ell,q}_{L}(\mathbb{R})$ for
the classical Sobolev space $W^{\ell,q}(\mathbb{R})$ restricted to
$2L$-periodic functions on the open interval $(-L,L)$. For $q=2$, we simply
denote $H^{\ell}_{L}(\mathbb{R})$.

To seek $2L$-periodic solutions, we consider the following periodic problem

$\begin{cases}v=(-\Delta)_{\rm cyl}^{-\sigma,L}({f}_{\sigma}\circ v)\quad{\rm
in}\quad\mathbb{R},\\\ \lim_{t\rightarrow+\infty}v(t)=0,\end{cases}$
($\mathcal{O}^{\prime}_{2\sigma,L}$)

where $(-\Delta)_{\rm
cyl}^{-\sigma,L}:\mathcal{C}_{L}^{0}(\mathbb{R})\rightarrow\mathcal{C}_{L}^{2\sigma}(\mathbb{R})$
is the integral periodic linear operator defined by

$(-\Delta)_{\rm cyl}^{-\sigma,L}(f_{\sigma}\circ v)(t):={\rm
p.v.}\int_{-L}^{L}f_{\sigma}(v(\tau))\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)\mathrm{d}\tau.$

For this, we shall work with the norm given by

$\|v\|_{H_{L,0}^{\sigma}(\mathbb{R})}:=\left([v^{(m)}]_{L_{L}^{s}(\mathbb{R})}+\sum_{\ell=0}^{m}\|v^{(\ell)}\|^{2}_{L_{L}^{2}(\mathbb{R})}\right)^{1/2},$

where

$[v^{(m)}]_{L_{L}^{s}(\mathbb{R})}:=\int_{-L}^{L}\int_{-L}^{L}[v^{(m)}(t)-v^{(m)}(\tau)]^{2}\widehat{\mathcal{K}}_{s,L}(t-\tau)\mathrm{d}\tau\mathrm{d}t.$

We also define the following higher-order functional space

$\displaystyle
H_{L}^{\sigma}(\mathbb{R})=\\{v\in\mathcal{C}^{2\sigma}_{L}(\mathbb{R}):\
\|v\|_{H_{L}^{\sigma}(\mathbb{R})}<\infty\\}.$

Furthermore, the evenness and periodicity

$\displaystyle H_{L,*}^{\sigma}(\mathbb{R})=\\{v\in
H_{L}^{\sigma}(\mathbb{R}):v(t)=v(-t)\ {\rm and}\ v(t+2{L})=v(t)\ {\rm for\
all}\ t\in\mathbb{R}\\}.$

As well as, taking into consideration the boundary condition

$\displaystyle H_{L,0}^{\sigma}(\mathbb{R})=\\{v\in
H_{L}^{\sigma}(\mathbb{R}):v^{(\ell)}(-L)=v^{(\ell)}(L)=0\ {\rm for}\
\ell\in\\{1,\dots,m\\}\\}.$

Finally, a suitable space to work is

$\displaystyle H_{L,0,*}^{\sigma}(\mathbb{R})=H_{L,0}^{\sigma}(\mathbb{R})\cap
H_{L,*}^{\sigma}(\mathbb{R}).$

Here $\widehat{\mathcal{K}}_{s,L}:(-L,L)\times(-L,L)\rightarrow\mathbb{R}$
given by

$\widehat{\mathcal{K}}_{s,L}(t-\tau)=\sum_{j\in\mathbb{Z}}\widehat{\mathcal{K}}_{s}(t-\tau-
jL)$ (4.15)

is a periodic Kernel, where
$\widehat{\mathcal{K}}_{s}:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ is
defined as (4.9).

Now we will introduce some standard Hölder fractional from [28, Theorem 8.2].

######  Lemma C.

Let $s\in(0,1)$ and $n>2s$. Assume that $p\in[1,+\infty)$. Then, there exists
a constant $C>0$, depending on $\sigma$ and $p$, such that

$\|v\|_{\mathcal{C}^{0,\alpha}_{L}(\mathbb{R})}\leqslant
C\left(\|v\|^{p}_{L^{p}_{L}(\mathbb{R})}+\int_{-L}^{L}\int_{-L}^{L}\frac{|v(t)-v(\tau)|^{p}}{|t-\tau|^{1+sp}}\mathrm{d}t\mathrm{d}\tau\right)^{\frac{1}{p}}$
(4.16)

for any $v\in L^{p}_{L}(\mathbb{R})$, where $\alpha=s-\frac{1}{p}$.

######  Lemma 4.8.

Let $\sigma\in(1,+\infty]$ and $n>2\sigma$. Assume that $p\in[1,+\infty)$ and
$\sigma\in(1,\frac{n}{2})$ is such that $\sigma-\frac{1}{p}\notin\mathbb{Z}$.
Then, there exists a constant $C>0$, depending on $\sigma$ and $p$, such that

$\|v\|_{\mathcal{C}^{\ell,\alpha}_{L}(\mathbb{R})}\leqslant
C\left(\|v\|^{p}_{W^{m,p}_{L}(\mathbb{R})}+\int_{-L}^{L}\int_{-L}^{L}\frac{|v^{(m)}(t)-v^{(m)}(\tau)|^{p}}{|t-\tau|^{1+sp}}\mathrm{d}t\mathrm{d}\tau\right)^{\frac{1}{p}}$
(4.17)

for any $v\in W^{\sigma,p}_{L}(\mathbb{R})$, where
$\ell=\lfloor\sigma-\frac{1}{p}\rfloor$ and
$\alpha=\sigma-\frac{1}{p}-\lfloor\sigma-\frac{1}{p}\rfloor$.

###### Proof.

It is a direct consequence of Lemma C by using a standard induction argument.
∎

We also need the following strong maximum principle.

######  Lemma 4.9.

Let $\sigma\in(1,+\infty]$ and $n>2\sigma$. If $v\in
H_{L}^{\sigma}(\mathbb{R})\cap\mathcal{C}^{0}(\mathbb{R})$ is a nonnegative
solution to ($\mathcal{O}_{2\sigma,L}^{\prime}$). Then, either $v>0$ or
$v\equiv 0$.

###### Proof.

Indeed, since $v\geqslant 0$, it follows that

$v=(-\Delta)_{\rm cyl}^{-\sigma}({f}_{\sigma}\circ v)\geqslant 0.$ (4.18)

Assume that there exists a point $t_{0}\in\mathbb{R}$ with $v(t_{0})=0$, then

$\displaystyle v(t_{0})-v(t)=(-\Delta)_{\rm
cyl}^{-\sigma}[{f}_{\sigma}(v(t_{0}))-{f}_{\sigma}(v(t))]$ $\displaystyle=$
$\displaystyle\mathrm{p.v.}\int_{-\infty}^{+\infty}{f}_{\sigma}(v(t_{0}))\widehat{\mathcal{R}}_{\sigma}(t_{0}-\tau)\mathrm{d}\tau-\mathrm{p.v.}\int_{-\infty}^{+\infty}{f}_{\sigma}(v(t))\widehat{\mathcal{R}}_{\sigma}(t-\tau)\mathrm{d}\tau$
$\displaystyle=$
$\displaystyle-\mathrm{p.v.}\int_{-\infty}^{+\infty}{f}_{\sigma}(v(\tau))\widehat{\mathcal{R}}_{\sigma}(t-\tau)\mathrm{d}\tau\leqslant
0$

satisfies (4.18) only in the case $v\equiv 0$. ∎

Now we have the most important lemma in this section

######  Lemma 4.10.

Let $\sigma\in(1,+\infty)$ and $n>2\sigma$. For any $L\gg 1$ sufficiently
large, there exist there exist a sequence of periods
$(L_{j})\in\ell^{\infty}(\mathbb{R}_{+})$, an error function
$\psi_{(0,L_{j})}\in H^{\sigma}_{L}(\mathbb{R})$ and a unique positive even
solution $\bar{v}_{(0,L_{j})}\in H_{L}^{\sigma}(\mathbb{R})$ to the following
periodic boundary value problem

$\left\\{\begin{aligned} v=(-&\Delta)_{\rm cyl}^{-\sigma,L}({f}_{\sigma}\circ
v)\quad\text{ in }\quad(-L,L),\\\ v^{(\ell)}(-L)&=v^{(\ell)}(L)=0\quad\text{
for }\quad\ell=1,3,\dots,2m-1,\end{aligned}\right.$
($\overline{{\mathcal{O}}}^{\prime}_{2\sigma,L}$)

which satisfy

$\bar{v}_{(0,L_{j})}(t)=\widehat{V}^{+}_{(0,L_{j})}(t)+\psi_{(0,L_{j})}(t)$

and

$\|\psi_{(0,L_{j})}\|_{H^{\sigma}(\mathbb{R})}\rightarrow 0\quad{\rm as}\quad
L\rightarrow+\infty,$

where $\widehat{V}^{+}_{(0,L_{j})}=\sum_{j\in\mathbb{Z}}V_{(0,L_{j})}(t)$ with
$V_{(0,L_{j})}(t)=\cosh(t-L_{j})$ and $L_{j}=2jL$ for $j\in\mathbb{Z}$ is the
standard bubble tower solution $($see Definition 5.2$)$ . Moreover, we have
the following Hölder estimate

$\|\psi_{(0,L_{j})}\|_{\mathcal{C}^{2\sigma}_{L}(\mathbb{R})}\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}$ (4.19)

for some $\alpha\in(0,1)$ and $\xi>0$ independent of the period $L\gg 1$
large.

###### Proof.

First, by symmetry $\bar{v}_{L}\in\mathcal{C}^{2\sigma}_{L}(\mathbb{R})$ given
by (5.2) satisfies the boundary condition at $t=\pm L$, that is,
$\widehat{V}^{+}_{(0,L_{j})}\in H_{(0,L_{j})}^{\sigma}(\mathbb{R})$. Now
writing $v=\widehat{V}^{+}_{(0,L_{j})}+\psi$, we can reformulate
($\mathcal{O}^{\prime}_{2\sigma,L}$) as

$\mathscr{N}_{\sigma,L}(\widehat{V}^{+}_{(0,L_{j})}+\psi)=0\quad{\rm
in}\quad\mathbb{R},$

where

$\displaystyle\mathscr{N}_{\sigma,L}(v):=v-(-\Delta)_{\rm
cyl}^{-\sigma,L}(f_{\sigma}\circ v).$ (4.20)

From now on, let us fix the notation

$\displaystyle\mathscr{N}_{\sigma}(0,L_{j})(\psi):=\mathscr{N}_{\sigma,L}(\widehat{V}^{+}_{(0,L_{j})}+\psi)$

Next, by linearizing this functional around the standard bubble tower
solution, we find

$\mathscr{L}_{\sigma}(0,L_{j})(\psi)=\mathscr{E}_{\sigma}(0,L_{j})(\widehat{V}^{+}_{(0,L_{j})})+\mathscr{S}_{\sigma}(0,L_{j})(\psi),$
(4.21)

where $\mathscr{L}_{\sigma}(0,L_{j}):H_{L}^{\sigma}(\mathbb{R})\to
H_{L}^{\sigma}(\mathbb{R})$ defined as
$\mathscr{L}_{\sigma}(0,L_{j}):=\mathrm{d}\mathscr{N}_{\sigma}[\widehat{V}^{+}_{(0,L_{j})}]$
satisfies

$\mathscr{L}_{\sigma}(0,L_{j})(\psi):=\psi-\mathscr{K}_{\sigma}(0,L_{j})(\psi),$
(4.22)

where

$\displaystyle\mathscr{K}_{\sigma}(0,L_{j})(\psi):=(-\Delta)^{-\sigma}(f^{\prime}_{\sigma}\circ\widehat{V}^{+}_{(0,L_{j})})\psi=\int_{-L}^{L}f^{\prime}_{\sigma}(\widehat{V}^{+}_{(0,L_{j})})\psi\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)\mathrm{d}\tau$
(4.23)

represents the derivative of the nonlinear functional (4.20) at the standard
bubble tower solution (5.2). Also, the superlinear term
$\mathscr{L}_{\sigma}(0,L_{j}):H_{L}^{\sigma}(\mathbb{R})\to
H_{L}^{\sigma}(\mathbb{R})$ is given by

$\mathscr{S}_{\sigma}(0,L_{j})(\psi)=\int_{-L}^{L}\left[f_{\sigma}(\widehat{V}^{+}_{(0,L_{j})}+\psi)-f_{\sigma}(\widehat{V}^{+}_{(0,L_{j})})-f^{\prime}_{\sigma}(\widehat{V}^{+}_{(0,L_{j})})\psi\right]\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)\mathrm{d}\tau.$

is a superlinear term, and the remainder error term is given by

$\mathscr{E}_{\sigma}(0,L_{j})(\widehat{V}^{+}_{(0,L_{j})})=\int_{-L}^{L}\left[f_{\sigma}\left(\sum_{j\in\mathbb{Z}}V_{(0,L_{j})}(\tau)\right)-\sum_{j\in\mathbb{Z}}f_{\sigma}\left(V_{(0,L_{j})}(\tau)\right)\right]\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)\mathrm{d}\tau,$
(4.24)

which represents the error in approximating a solution to (5.2) by a standard
bubble-tower solution.

To apply classical Fredholm theory we need to prove the following claim:

Claim 1: The operator
$\mathscr{K}_{\sigma}(0,L_{j}):H_{L}^{\sigma}(\mathbb{R})\to
H_{L}^{\sigma}(\mathbb{R})$ defined in (4.23) is bounded and compact and
satisfies

$\left\|\mathscr{K}_{\sigma}(0,L_{j})(\psi)\right\|_{H^{\sigma}_{L}(\mathbb{R})}\lesssim\left\|\psi\right\|_{L^{2}_{L}(\mathbb{R})}$
(4.25)

uniformly on $L\gg 1$ large.

Initially, we prove that the operator is bounded and compact. Indeed, by its
definition,
$\widehat{\mathcal{K}}_{\sigma,L}\in\mathcal{C}^{j,\alpha}(\mathbb{R})$ for
any $j\in\\{0,\dots,m\\}$ and some $\alpha\in(0,1)$. Moreover, for each $j$,
we have

$\left|\frac{\mathrm{d}^{j}}{\mathrm{d}t^{j}}\widehat{\mathcal{K}}_{\sigma,L}(t)\right|\lesssim
e^{-c_{j}|t|}$ (4.26)

uniformly on $L\gg 1$ large. Similarly, we have

$|\widehat{V}^{+}_{(0,L_{j})}|\lesssim e^{-c|t|}$

uniformly on $L\gg 1$ large. By Hölder’s inequality, it follows directly

$\left\|\mathscr{K}_{\sigma}(0,L_{j})(\psi)^{(j)}\right\|_{L_{L}^{2}(\mathbb{R})}\lesssim\left\|\psi\right\|_{L_{L}^{2}(\mathbb{R})}\quad{\rm
for\ all}\quad j\in\\{0,\dots,m\\}$

uniformly on $L\gg 1$ large.

Also, using the Hölder continuity of
$(\widehat{\mathcal{K}}_{\sigma,L})^{(m)}\in\mathcal{C}^{0,\alpha+s}(\mathbb{R})$,
we have

$\displaystyle\left|\mathscr{K}_{\sigma}(0,L_{j})(\psi)^{(m)}(\tau)-\mathscr{K}_{\sigma}(0,L_{j})(\psi)^{(m)}(t)\right|$
$\displaystyle\leqslant$
$\displaystyle\int_{-L}^{L}f^{\prime}_{\sigma}(\widehat{V}^{+}_{(0,L_{j})}(\xi))\left|\frac{\mathrm{d}^{m}}{\mathrm{d}\xi^{m}}\left(\widehat{\mathcal{K}}_{\sigma,L}(\tau-\xi)-\widehat{\mathcal{K}}_{\sigma,L}(t-\xi)\right)\right||\psi(\xi)|\mathrm{d}\xi$
$\displaystyle\lesssim$
$\displaystyle\int_{-L}^{L}|\psi(\xi)||t-\tau|^{\alpha}\mathrm{d}\xi.$

Thus, by using the asymptotic behavior of the kernel near the origin given by
(4.11) and (4.26) combined with the last inequality, we obtain

$\displaystyle\left[\mathscr{K}_{\sigma}(0,L_{j})(\psi)\right]_{L^{s}_{L}(\mathbb{R})}$
$\displaystyle=\int_{-L}^{L}\left|\mathscr{K}_{\sigma}(0,L_{j})(\psi)^{(m)}(\tau)-\mathscr{K}_{\sigma}(0,L_{j})(\psi)^{(m)}(t)\right|^{2}\widehat{\mathcal{K}}_{\sigma,L}(t-\tau)\mathrm{d}\tau\mathrm{d}t$
$\displaystyle\lesssim\left\|\psi\right\|_{L^{2}_{L}(\mathbb{R})},$

uniformly on $L\gg 1$ large, which proves (4.25). In conclusion, by compact
embedding, the desired conclusion holds for the map
$\mathscr{K}_{\sigma}(0,L_{j}):H_{L}^{\sigma}(\mathbb{R})\to
H_{L}^{\sigma}(\mathbb{R})$.

The proof of the first claim is now finished.

Second, now in order to apply Fredholm alternative to conclude that for any
$h\in L^{2}_{L}(\mathbb{R})$, there exists a unique solution
$\psi_{(0,L_{j})}\in H_{L}^{\sigma}(\mathbb{R})$ to the linear inhomogeneous
problem

$\psi_{(0,L_{j})}-\mathscr{K}_{\sigma}(0,L_{j})(\psi_{(0,L_{j})})=h\quad{\rm
in}\quad(-L,L).$

One needs to prove the uniqueness result below:

Claim 2: The linear homogeneous equation

$\psi-\mathscr{K}_{\sigma}(0,L_{j})(\psi)=0\quad{\rm in}\quad(-L,L)$

admits only zero solutions in $L^{2}_{L}(\mathbb{R})$.

As a matter of fact, note that the equation above with Hölder’s inequality
yields directly that

$\left\|\psi\right\|_{L^{\infty}_{L}(\mathbb{R})}\lesssim\left\|\psi\right\|_{L^{2}_{L}(\mathbb{R})}.$
(4.27)

Next, we use the nondegeneracy of the standard bubble solution in Lemma 4.6 to
conclude that $\psi\equiv 0$, and thus we prove Claim 2.

Lastly, by the standard fixed-point argument, a unique solution
$\psi_{(0,L_{j})}\in H_{L}^{\sigma}(\mathbb{R})$ to (4.21) satisfying the
estimate

$\|\psi_{(0,L_{j})}\|_{H_{L,0}^{\sigma}(\mathbb{R})}\lesssim\|\mathscr{E}_{\sigma}(0,L_{j})(\widehat{V}^{+}_{(0,L_{j})})\|_{L_{L}^{2}(\mathbb{R})},$
(4.28)

to conclude the proof of the proposition, we are left to obtain estimates for
the right-hand side of the last inequality.

This is the content of our third claim.

Claim 3: It holds that $\|\psi\|_{H_{L,0}^{\sigma}(\mathbb{R})}\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}$ for some $\xi>0$ uniformly on $L\gg 1$ large.

In fact, using (5.2) it follows

$\mathscr{E}_{\sigma}(0,L_{j})(\widehat{V}^{+}_{(0,L_{j})})=c_{n,\sigma}\int_{-L}^{L}\left[\left(\sum_{j\in\mathbb{Z}}V_{(0,L_{j})}(\tau)\right)^{\frac{n+2\sigma}{n-2\sigma}}-\left(\sum_{j\in\mathbb{Z}}V_{(0,L_{j})}(\tau)^{\frac{n+2\sigma}{n-2\sigma}}\right)\right]\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)\mathrm{d}\tau.$

Since by symmetry, we have $V_{(0,-L_{j})}(t)\leqslant V_{(0,L_{j})}(t)$ for
$t\geqslant 0$, it holds

$\displaystyle|\mathscr{E}_{\sigma}(0,L_{j})(\widehat{V}^{+}_{(0,L_{j})})|$
$\displaystyle\lesssim\int_{-L}^{L}\left(V_{(0,\infty)}^{\frac{4\sigma}{n-2\sigma}}\sum_{j\in\mathbb{Z}^{*}}V_{(0,L_{j})}(\tau)+\sum_{j\in\mathbb{Z}^{*}}V_{(0,L_{j})}(\tau)^{\frac{n+2\sigma}{n-2\sigma}}\right)\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)\mathrm{d}\tau$
(4.29)
$\displaystyle\lesssim\sum_{j\in\mathbb{Z}^{*}}\int_{-L}^{L}V_{(0,\infty)}(\tau)^{\frac{4\sigma}{n-2\sigma}}V_{(0,L_{j})}(\tau)\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)\mathrm{d}\tau+\sum_{j\in\mathbb{Z}^{*}}\int_{-L}^{L}V_{(0,L_{j})}(\tau)^{\frac{n+2\sigma}{n-2\sigma}}\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)\mathrm{d}\tau.$

From (4.29), we find

$\displaystyle\int_{-L}^{L}|\mathscr{E}_{\sigma}(0,L_{j})(\widehat{V}^{+}_{(0,L_{j})})|^{2}\mathrm{d}t$
$\displaystyle\lesssim\int_{-L}^{L}\left(\sum_{j\in\mathbb{Z}^{*}}\int_{-L}^{L}V_{(0,\infty)}(\tau)^{\frac{8\sigma}{n-2\sigma}}V_{(0,L_{j})}(\tau)^{2}\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)^{2}\mathrm{d}\tau+\sum_{j\in\mathbb{Z}^{*}}\int_{-L}^{L}V_{(0,L_{j})}(\tau)^{\frac{2n+4\sigma}{n-2\sigma}}\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)^{2}\mathrm{d}\tau\right)\mathrm{d}t$
$\displaystyle\lesssim\sum_{j\in\mathbb{Z}^{*}}\int_{-L}^{L}\int_{-L}^{L}V_{(0,\infty)}(\tau)^{\frac{8\sigma}{n-2\sigma}}V_{(}0,L_{j})(\tau)^{2}\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)^{2}\mathrm{d}\tau\mathrm{d}t+\sum_{j\in\mathbb{Z}^{*}}\int_{-L}^{L}\int_{-L}^{L}V_{(0,L_{j})}(\tau)^{\frac{2n+4\sigma}{n-2\sigma}}\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)^{2}\mathrm{d}\tau$
$\displaystyle=:I_{1}+I_{2}.$ (4.30)

To estimate these two terms, we fix $\alpha\in(0,1)$ and subdivide
$\mathbb{R}=\\{|t|\leqslant\alpha
L\\}\mathbin{\dot{\cup}}\\{|t|\geqslant\alpha L\\}$. Then, we use the
exponential decay of the standard bubble solution from Proposition B to obtain

$\sum_{j\in\mathbb{Z}^{*}}V_{(0,L_{j})}\lesssim
e^{-\gamma_{\sigma}L(2-\alpha)}\quad{\rm
and}\quad\sum_{j\in\mathbb{Z}^{*}}V_{(0,L_{j})}\lesssim
e^{-\gamma_{\sigma}L}.$ (4.31)

Hence, by substituting in (4.31) into the first term in (4.2), we obtain

$\displaystyle I_{1}$
$\displaystyle=\sum_{j\in\mathbb{Z}^{*}}\int_{-L}^{L}\int_{-L}^{L}V_{(0,\infty)}(\tau)^{\frac{8\sigma}{n-2\sigma}}V_{(0,L_{j})}(\tau)^{2}\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)^{2}\mathrm{d}\tau\mathrm{d}t$
$\displaystyle\lesssim
e^{-2\gamma_{\sigma}L(2-\alpha)}+e^{-2\gamma_{\sigma}L\left(\frac{4\sigma\alpha}{n-2\sigma}\right)}e^{-2\gamma_{\sigma}L}e^{-2\gamma_{\sigma}L\left(\frac{n+2\sigma}{n-2\sigma}\right)}$
(4.32) $\displaystyle\lesssim
e^{-2\gamma_{\sigma}L(2-\alpha)}+e^{-2\gamma_{\sigma}L(1+\xi)}$

for some $\xi>0$ (depending only on $n$, $\sigma$, and $\alpha$), where we
used the asymptotic behavior of the Kernel (4.12) for $\tau\rightarrow+\infty$
large and the fact that it is bounded for $\tau\rightarrow 0$ small.

Furthermore, by substituting (4.31) into the second term in (4.2), we have

$\displaystyle
I_{2}=\sum_{j\in\mathbb{Z}^{*}}\int_{-L}^{L}\int_{-L}^{L}V_{(0,L_{j})}(\tau)^{\frac{2n+4\sigma}{n-2\sigma}}\widehat{\mathcal{R}}_{\sigma,L}(t-\tau)^{2}\mathrm{d}\tau\mathrm{d}t\lesssim
e^{-2\gamma_{\sigma}L(2-\alpha)}+e^{-2\gamma_{\sigma}L(\frac{n+2\sigma}{n-2\sigma})}.$
(4.33)

In conclusion, by substituting (4.2) and (4.33) into (4.29), we have

$\|\mathscr{E}_{\sigma}(0,L_{j})(\widehat{V}^{+}_{(0,L_{j})})\|_{L^{2}(\mathbb{R})}\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}$

for some $\xi>0$ uniformly on $L\gg 1$ large, which combined with (4.28)
proves the third claim.

Finally, by standard estimates in Lemma 4.8 combined with the regularity
lifting theorem from [23, Theorem 3.3.1] applied to (4.21), it follows that
$\psi_{(0,L_{j})}\in H_{L,0,*}^{\sigma}(\mathbb{R})$ is smooth and satisfies

$\|\psi_{(0,L_{j})}\|_{\mathcal{C}^{2\sigma+\alpha}(\mathbb{R})}\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}$

for some $\xi>0$ independent of $L\gg 1$ large.

Therefore, the maximum principle in Lemma 4.9 concludes the proof of the
proposition. ∎

######  Remark 4.11.

It is worth noticing that the last proof differs in spirit from the local
inversion technique in [36, Proposition 2.3]. Instead of using this method, we
give an alternative proof based on the dual formulation from Lemma 3.3. This
technique is of independent interest to a larger class of integral equations
not necessarily arising as the dual of a differential equation.

######  Remark 4.12.

We notice that it is straightforward to extend the local inversion method in
[36, Proposition 2.3] at least for the higher order local cases
$\sigma=m\in\mathbb{N}$. For see this fact, we write the poly-harmonic
operator in Emden–Fowler coordinates, which gives us

$(-\Delta)_{\rm cyl}^{m}:=(-\Delta)_{\rm rad}^{m}+(-\Delta)_{\rm ang}^{m}$

with

$(-\Delta)^{m}_{\rm
rad}:=\partial_{t}^{(m)}-K^{(0)}_{2m-2}\partial_{t}^{(2m-2)}+\cdots+(-1)^{m}K_{1}^{(0)}\partial_{t}^{(1)}+(-1)^{m+1}K_{0}^{(0)},$

and

$(-\Delta)^{m}_{\rm
ang}:=\sum_{\ell=1}^{2m}\sum_{j=0}^{2m}(-1)^{\frac{j+2}{2}}K^{(\ell)}_{2m,j}\partial_{t}^{(j)}(-\Delta_{\theta})^{\ell},$

where $K^{(\ell)}_{2m,j}=K^{(\ell)}_{2m,j}(n)>0$ for $j\in\\{0,\dots,2m\\}$
and $\ell\in\\{1,\dots,2m\\}$ are dimensional constants. For this computation,
we refer the interested reader to [5]. After that, we need to build on the
classification result from standard bubbles for the critical Sobolev embedding
$H^{m}(\mathbb{R}^{n})\rightarrow L^{2_{m}^{*}}(\mathbb{R}^{n})$, where
$2_{m}^{*}:=\frac{2n}{n-2m}$ from [54] and a standard nondegeneracy technique
as in Lemma 4.6.

As an immediate consequence of the last proposition, one has

######  Corollary 4.13.

Let $\sigma\in(1,+\infty)$ and $n>2\sigma$. For any $L\gg 1$ sufficiently
large, there exist there exist a sequence of periods
$(L_{j})\in\ell^{\infty}(\mathbb{R}_{+})$, an error function
$\psi_{(0,L_{j})}\in H^{\sigma}_{L}(\mathbb{R})$ and a unique positive even
periodic solution $\bar{v}_{(0,L_{j})}\in H_{L}^{\sigma}(\mathbb{R})$ to
($\mathcal{O}_{2\sigma,L}^{\prime}$) satisfying

$\bar{v}_{(0,L_{j})}(t)=\widehat{V}^{+}_{(0,L_{j})}(t)+\psi_{(0,L_{j})}(t)$

and

$\|\psi_{(0,L_{j})}\|_{H_{L}^{\sigma}(\mathbb{R})}\rightarrow 0\quad{\rm
as}\quad L\rightarrow+\infty,$

where $\widehat{V}^{+}_{(0,L_{j})}\in\mathcal{C}^{2\sigma}(\mathbb{R})$ is the
standard bubble tower solution given by (5.5). Moreover, we have the following
Hölder estimate

$\|\psi_{(0,L_{j})}\|_{\mathcal{C}^{2\sigma}_{L}(\mathbb{R})}\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}$ (4.34)

for some $\alpha\in(0,1)$ and $\xi>0$ independent of the period $L\gg 1$
large.

Since ($\mathcal{O}_{2\sigma,\infty}$) is translational invariant, we now will
use the periodic solution $\bar{v}_{(0,L_{j})}$ which attains its minimum at
the points $t=2jL$ with $j\in\mathbb{Z}$. Indeed, using Lemma 4.10, this
periodic solution can be expressed as a perturbation of a bubble tower with
estimated error.

######  Definition 4.14.

Let $\sigma\in(1,+\infty)$ and $n>2\sigma$. For any $L\gg 1$ sufficiently
large, let us define the generalized bubble tower solution

  * (i)

$($Emden–Fowler coordinates$)$

$\bar{v}_{(0,L_{j})}(t):=\widehat{V}^{+}_{(0,L_{j})}(t)+\psi_{(0,L_{j})}(t),$
(4.35)

where $\widehat{V}^{+}_{(0,L_{j})}\in\mathcal{C}^{2\sigma}(\mathbb{R})$ is the
standard half bubble tower solution given by (5.6) and
$\psi_{(0,L_{j})}\in\mathcal{C}^{2\sigma}(\mathbb{R})$ the perturbation
function constructed in Corollary 4.13. More precisely, one has

$\widehat{V}^{+}_{(0,L_{j})}(t)=\sum_{j\in\mathbb{N}}\cosh(t-L_{j}-L)^{\gamma_{\sigma}},\quad{\rm
where}\quad L_{j}=(1+2j)L\quad{\rm for}\quad j\in\mathbb{N}.$

  * (ii)

$($Spherical coordinates$)$

$\bar{u}_{(0,L_{j})}(x):=\widehat{U}^{+}_{(0,L_{j})}(x)+\phi_{(0,L_{j})}(x),$
(4.36)

where
$\widehat{U}^{+}_{(0,L_{j})}\in\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\Sigma)$
is the standard half bubble tower solution given by (5.3) and
$\phi_{(0,L_{j})}\in\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\Sigma)$ is
the perturbation function constructed in Corollary 5.1. More precisely, we
have

$\widehat{U}^{+}_{(0,L_{j})}(x)=\sum_{j\in\mathbb{N}}\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+|x|^{2}}\right)^{\gamma_{\sigma}},\quad{\rm
where}\quad\lambda_{j}=e^{-(1+2j)L}\quad{\rm for}\quad j\in\mathbb{N}.$

We find better asymptotics near the isolated singularities for the deformed
solution obtained in Lemma 4.10. These refined estimates in terms of the
bubble tower solution will be a crucial part of estimating the errors in our
approximate solution in the gluing procedure in Section 7.

######  Lemma 4.15.

The asymptotics holds

$\bar{v}_{(0,L_{j})}(t)=v_{\rm sph}(t)(1+\mathrm{o}(1))\quad{\rm as}\quad
L\rightarrow+\infty,$ (4.37)

or undoing the Emden–Fowler change of variables, it holds

$\bar{u}_{(0,L_{j})}(t)=u_{\rm sph}(|x|)(1+\mathrm{o}(1))\quad{\rm as}\quad
L\rightarrow+\infty.$

Moreover, one has

$\bar{u}_{(0,L_{j})}(x)=|x|^{n-2\sigma}e^{-\gamma_{\sigma}L}(1+\mathrm{o}(1))\quad{\rm
as}\quad L\rightarrow+\infty.$ (4.38)

and

$\varepsilon_{L}:=\bar{v}_{(0,L_{j})}(0)=e^{-\gamma_{\sigma}L}(1+\mathrm{o}(1))\quad{\rm
as}\quad L\rightarrow+\infty.$

This parameter is called the neck size or Delaunay parameter.

###### Proof.

Notice that for $t\leqslant 0$ (or $|x|\geqslant 1$), the proof of the lemma
follows as a combination of Corollary 4.13 together with exponential decay of
the standard spherical solution in Emden–Fowler coordinates to prove (4.37). ∎

##  5\. Approximate solution

This section will construct a suitable approximate solution to
(${\mathcal{Q}}^{\prime}_{2\sigma,\Sigma}$). We also prove some estimates
concerning the behavior of such a solution near the singular set. As we have
mentioned, one of the main ideas is that, although we would like the
approximate solution to have Delaunay-type singularities around each point
isolated singularity, it should have a fast decay once we are away from the
singular set to glue to the flat background manifold. To this end, we will
only take half a Delaunay solution (this is, only values $j\in\mathbb{N}$).

###  5.1. Local asymptotic behavior

In this subsection, we study the local behavior of solutions to
($\mathcal{Q}^{\prime}_{2\sigma,\infty}$) near the isolated singularity at the
origin. Namely, we show that near the origin, it can be approximated by a
bubble tower solution. In contrast with the cases $\sigma\in\\{1,2,3\\}$ on
which a complete classification of this local behavior is given in terms of
the two-parameter family of Delaunay solutions studied in [18, 30, 5], which
are inspired by the classical result of Korevaar et al. [40] for $\sigma=1$
and Caffarelli et al. $\sigma\in(0,1)$ [16]. These will be the building blocks
in constructing suitable approximate solutions to
($\mathcal{Q}^{\prime}_{2\sigma,\infty}$).

First, recall the local asymptotic classification result from [37].

######  Proposition B.

Let $\sigma\in(1,+\infty)$ and $n>2\sigma$.

  * (i)

Assume that $R=+\infty$. For any $L\gg 1$ sufficiently large, there exists a
blow-up limit solution to ($\mathcal{Q}^{\prime}_{2\sigma,\infty}$) denoted by
$u_{L}\in\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\\{0\\})$ and given by

$u_{(0,L)}(x)=\left(\mathfrak{F}_{\sigma}\right)^{-1}(v_{(0,L)})=|x|^{-\gamma_{\sigma}}v_{(0,L)}(-\ln|x|),$

where $v_{L}\in\mathcal{C}^{2\sigma}(\mathbb{R})$ is a bounded periodic even
solution to ($\mathcal{O}_{2\sigma,L}^{\prime}$). In addition, one has

$v_{(0,L)}(x)=\mathcal{O}(e^{-\gamma_{\sigma}L})\quad{\rm as}\quad
t\rightarrow+\infty.$

These will be called Delaunay solutions.

  * (ii)

Assume that $0<R<+\infty$. If $u\in\mathcal{C}^{2\sigma}(B_{R}^{*})$ is a
positive singular solution to ($\mathcal{Q}^{\prime}_{2\sigma,R}$), then there
exists a Delaunay solution with a large period, denoted by $u_{L}$, such that

$u(x)=u_{(0,L)}(x)(1+\mathrm{o}(1))\quad{\rm as}\quad|x|\rightarrow 0,$

or

$v(t)=v_{(0,L)}(t)(1+\mathrm{o}(1))\quad{\rm as}\quad t\rightarrow+\infty,$

where $L\gg 1$ is sufficiently large.

In addition, writing Lemma 4.15 into using Emden–Fowler change of variables,
we can reformulate it as an improvement for the result above.

######  Corollary 5.1.

Let $\sigma\in(1,+\infty)$ and $n>2\sigma$. For any $L\gg 1$ sufficiently
large, there exist a sequence of periods
$(L_{j})\in\ell^{\infty}(\mathbb{R}_{+})$ such that $\phi_{(0,L_{j})}\in
H^{\sigma}_{L}(\mathbb{R}^{n}\setminus\\{0\\})$ and a unique positive even
solution $\bar{u}_{(0,L_{j})}\in
H_{L}^{\sigma}(\mathbb{R}^{n}\setminus\\{0\\})$ to
($\mathcal{Q}^{\prime}_{2\sigma,\infty}$) satisfying

$\bar{u}_{(0,L_{j})}(x)=\widehat{U}^{+}_{(0,L_{j})}(x)+\phi_{(0,L_{j})}(x),$

and

$\|\phi_{(0,L_{j})}\|_{H_{L}^{\sigma}(\mathbb{R}^{n}\setminus\\{0\\})}\rightarrow
0\quad{\rm as}\quad L\rightarrow+\infty,$

where
$\widehat{U}^{+}_{(0,L_{j})}\in\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\\{0\\})$
is the standard bubble tower solution given by (5.2). Moreover, we have the
following Hölder estimate

$\|\phi_{(0,L_{j})}\|_{\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\\{0\\})}\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}$ (5.1)

for some $\alpha\in(0,1)$ and $\xi>0$ independent of $L\gg 1$ large.

Based on the definition of a spherical solution in (4.2) and (4.13), we
introduce the concept of a standard bubble tower solution. In addition, in
order to have fast decay far from the singularity $(t\rightarrow-\infty)$, we
will need only half a bubble tower. This fact motivates the following
definition

######  Definition 5.2.

Let $\sigma\in(1,+\infty)$ and $n>2\sigma$. For any $L\gg 1$ sufficiently
large, let us define the following standard bubble tower solution

  * (i)

$($Spherical coordinates$)$

$\widehat{U}_{(0,L_{j})}(x):=\sum_{j\in\mathbb{Z}}U_{(0,L_{j})}(x),$ (5.2)

and

$\widehat{U}^{+}_{(0,L_{j})}(x):=\sum_{j\in\mathbb{N}}U_{(0,L_{j})}(x),$ (5.3)

where

$U_{(0,L_{j})}(x)=\left(\frac{\lambda_{j}}{\lambda_{j}^{2}+\left|x\right|^{2}}\right)^{\gamma_{\sigma}}\quad{\rm
with}\quad\lambda_{j}=e^{-2jL}\quad{\rm for}\quad j\in\mathbb{Z}.$ (5.4)

  * (ii)

$($Emden–Fowler coordinates$)$

$\widehat{V}_{(0,L_{j})}(t):=\sum_{j\in\mathbb{Z}}V_{(0,L_{j})}(t),$ (5.5)

and

$\widehat{V}^{+}_{(0,L_{j})}(t):=\sum_{j\in\mathbb{N}}V_{(0,L_{j})}(t),$ (5.6)

where

$V_{(0,L_{j})}(t)=\cosh(t-L_{j})^{\gamma_{\sigma}}\quad{\rm with}\quad
L_{j}=2jL\quad{\rm for}\quad j\in\mathbb{Z}.$ (5.7)

These will be called the standard bubble tower solution.

As a consequence of Corollary 4.13, we will improve the last two results.
Indeed, we show that near an isolated singularity, solutions are close to some
bubble tower solution up to some controlled error.

######  Proposition 5.3.

Let $\sigma\in(1,+\infty]$ with $n>2\sigma$. If
$u\in\mathcal{C}^{2\sigma}(B_{R}^{*})$ is a positive smooth singular solution
to ($\mathcal{Q}^{\prime}_{2\sigma,R}$) with $R>0$, then there exists a
sequence of periods $(L_{j})\in\ell^{\infty}(\mathbb{R}_{+})$ and a blow-up
limit solution
$\bar{u}_{(0,L_{j})}\in\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\\{0\\})$
to ($\mathcal{Q}^{\prime}_{2\sigma,\infty}$) such that

$u(x)=\bar{u}_{(0,L_{j})}(x)(1+\mathrm{o}(1))\quad{\rm as}\quad|x|\rightarrow
0.$

More precisely, one has

$\bar{u}_{(0,L_{j})}(x)=\widehat{U}_{(0,L_{j})}(x)+\phi_{(0,L_{j})}(x),$ (5.8)

where
$\widehat{U}_{(0,L_{j})}\in\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\\{0\\})$
is the standard bubble tower solution in (5.3) and $\phi_{(0,L_{j})}\in
H^{\sigma}(\mathbb{R}^{n})$ satisfies

$\|\phi_{(0,L_{j})}\|_{\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\\{0\\})}\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}$ (5.9)

for some $\alpha\in(0,1)$ and $\xi>0$ independent of $L\gg 1$ large.

###  5.2. Balanced configurations

Here we introduce a necessary set of compatibility conditions for the
configuration parameters.

######  Definition 5.4.

Let $\sigma\in(1,+\infty)$, $n>2\sigma$, and $N\geqslant 2$. Given $L\gg 1$
large enough, we will fix the vector
$\boldsymbol{L}=(L_{1},\ldots,L_{N})\in\mathbb{R}^{N}_{+}$ to be the Delaunay
parameters, which also are related to the neck sizes of each Delaunay
solution. They will be chosen (large enough) in the proof. They will satisfy
the following conditions $|L_{i}-L|\lesssim 1$ for all $i\in\\{1,\dots,N\\}$.
More precisely, they will be related by the vector
$\boldsymbol{q}=(q_{1},\dots,q_{N})\in\mathbb{R}^{N}_{+}$, which satisfy the
following relations

$q_{i}e^{-\gamma_{\sigma}L}=e^{-\gamma_{\sigma}L_{i}}\quad{\rm for}\quad
i\in\\{1,\dots,N\\}.$ (5.10)

Next, we will give some explanation about the choice of parameters. Given the
$N(n+2)$ balancing parameters
$(\boldsymbol{q}^{b},\boldsymbol{R}^{b},\boldsymbol{\hat{a}}_{0}^{b})$
satisfying the balancing conditions ($\mathscr{B}_{1}$) and
($\mathscr{B}_{2}$), we first choose $N(n+2)$ initial perturbation parameters
$(\boldsymbol{q},\boldsymbol{R},\boldsymbol{\hat{a}}_{0})$ which are close to
the balancing parameters, i.e., (5.11) and (5.12).

######  Definition 5.5.

Let $\sigma\in(1,+\infty)$, $n>2\sigma$, and $N\geqslant 2$. For any fixed
nonnegative vector
$\boldsymbol{q}^{b}=(q_{1}^{b},\ldots,q_{N}^{b})\in\mathbb{R}^{N}_{+}$, let us
define the vector
$(\mathbf{a}_{0}^{b},\boldsymbol{R}^{b})=(a_{0}^{i,b},\dots,a_{0}^{N,b},R^{1,b},\dots
R^{N,b})\in\mathbb{R}^{(n+1)N}$ to be determined by the following balancing
conditions:

$q_{i}^{b}=A_{2}\sum_{i^{\prime}\neq
i}q_{i^{\prime}}^{b}(R^{i,b}R^{i^{\prime},b})^{\gamma_{\sigma}}|x_{i}-x_{i^{\prime}}|^{-2\gamma_{\sigma}}\quad{\rm
for}\quad i\in\\{1,\dots,N\\}.$ ($\mathscr{B}_{1}$)

and

$\frac{a_{0}^{i,b}}{(\lambda_{0}^{i,b})^{2}}=-\frac{A_{3}}{A_{1}}\sum_{i^{\prime}\neq
i}\frac{x_{i^{\prime}}-x_{i}}{|x_{i^{\prime}}-x_{i}|^{2\gamma_{\sigma}+2}}\frac{q_{i^{\prime}}^{b}}{q_{i}^{b}}(R^{i,b}R^{i^{\prime},b})^{\gamma_{\sigma}}\quad{\rm
for}\quad i\in\\{1,\dots,N\\}$ ($\mathscr{B}_{2}$)

where $\lambda_{0}^{i,b}=R^{i,b}e^{-{L_{i}^{b}}}$, and the
$L_{i}^{b}\in\mathbb{R}_{+}$ are defined from the $q_{i}^{b}\in\mathbb{R}_{+}$
by (5.10) for each $i\in\\{1,\dots,N\\}$ and the constants
$A_{1},A_{2}>0,A_{3}<0$ are defined in (A.1), (A.2), and (A.3), respectively.
We denote by
$(\boldsymbol{q}^{b},\mathbf{a}_{0}^{b},\boldsymbol{R}^{b})\in{\rm
Bal}_{\sigma}(\Sigma)$ the set of balanced configurations.

######  Remark 5.6.

We remark that it has been shown in [47, Remark 3] that for
$\boldsymbol{q}:=(q_{1}^{b},\ldots,q_{N}^{b})\in\mathbb{R}^{N}_{+}$ in the
positive octant, there exists a solution
$\boldsymbol{R}^{b}=(R^{1,b},\dots,R^{N,b})$ to equation ($\mathscr{B}_{1}$).
Once this is chosen, then we can use equation ($\mathscr{B}_{2}$) to determine
$\mathbf{a}_{0}^{b}=(a_{0}^{1,b},\dots,a_{0}^{N,b})\in(\mathbb{R}^{n}_{+})^{N}$.
In other words, the set of balanced configurations is non-empty ${\rm
Bal}_{\sigma}(\Sigma)\neq\varnothing$ for all $\sigma\in\mathbb{R}_{+}$.

Although the meaning of these compatibility conditions will become apparent in
the following sections, we have just seen that they are analogous to those of
[46] for the local case. The idea is that perturbations at the base level
should be very close to those for a single bubble. This fact also shows, in
particular, that even though our problem is nonlocal, very near the
singularity, it presents a local behavior due to the strong influence of the
underlying geometry. However, for the rest of the perturbation parameters, we
must solve an infinite-dimensional system of equations.

The last discussion motivates the definition below

######  Definition 5.7.

Let $\sigma\in(1,+\infty)$, $n>2\sigma$ and $N\geqslant 2$. We define the so-
called configuration map $\Upsilon_{\rm
conf}:\mathbb{R}^{(n+1)N}\rightarrow\mathbb{R}^{(n+2)N}$ which associates
compatible moduli space parameters $(\boldsymbol{x},\boldsymbol{L})$ with
configuration parameters $(\boldsymbol{q},\boldsymbol{a_{0}},\boldsymbol{R})$.
We say that a set moduli space parameters
$(\boldsymbol{x},\boldsymbol{L})\in\mathbb{R}^{(n+1)N}$ is compatible if its
associated set of configuration parameters
$(\boldsymbol{q},\boldsymbol{a_{0}},\boldsymbol{R})\in{\rm
Bal}_{\sigma}(\Sigma)$ is balanced.

###  5.3. Admissible perturbation parameters

We also would like to introduce some perturbation parameters
$R\in\mathbb{R},a\in\mathbb{R}^{n}$, since each standard bubble has $n+1$ free
parameters corresponding to scaling and translations, which is done for each
bubble in the bubble tower independently. Thus, we will have an infinite-
dimensional set of perturbations.

######  Definition 5.8.

Let $\sigma\in(1,+\infty)$, $n>2\sigma$ and $N\geqslant 2$. For any $L\gg 1$
sufficiently large and
$\boldsymbol{L}=(L^{1},\ldots,L^{N})\in\ell^{\infty}(\mathbb{R}^{N}_{+})$ and
$\boldsymbol{R}=(R^{1},\ldots,R^{N})\in\ell^{\infty}(\mathbb{R}^{N}_{+})$ such
that (5.10) holds, let us define the full set of perturbation parameters
$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})=(a_{j}^{1},\dots,a_{j}^{N},\lambda_{j}^{1},\dots,\lambda_{j}^{N})\in\ell^{\infty}(\mathbb{R}^{(n+1)N})$,
where

$\lambda_{j}^{i}=R_{j}^{i}e^{-{(1+2j)L_{i}}}\quad{\rm for}\quad
i\in\\{1,\dots,N\\}\quad{\rm and}\quad j\in\mathbb{Z}.$

We introduce the perturbation parameters we will use in the gluing technique

######  Definition 5.9.

Let $\sigma\in(1,+\infty)$, $n>2\sigma$ and $N\geqslant 2$. Let
$\boldsymbol{R}=(R^{1},\dots,R^{N})\in\mathbb{R}^{N}_{+}$ and
$\boldsymbol{q}\in\mathbb{R}^{N}_{+}=(q_{1},\dots,q_{N})$ be such that

$|R^{i}-R^{i,b}|\lesssim 1\quad{\rm and}\quad|q_{i}-q_{i}^{b}|\lesssim
1\quad{\rm for}\quad i\in\\{1,\dots,N\\}.$ (5.11)

Also, we let
$\boldsymbol{\lambda}^{0}_{0}=(\lambda_{0}^{1,0},\dots,\lambda_{0}^{N,0})\in\mathbb{R}^{N}_{+}$
and
$\boldsymbol{\hat{a}}_{0}^{i}=({\hat{a}}_{0}^{1},\dots,{\hat{a}}_{0}^{N})\in(\mathbb{R}^{n})^{N}$
be respectively given by

$\lambda_{0}^{i,0}=R^{i}e^{-\frac{(1+2j)}{2}L_{i}}\quad{\rm for}\quad
i\in\\{1,\dots,N\\}$

and

$\frac{a_{0}^{i,0}}{(\lambda_{0}^{i,0})^{2}}=\hat{a}_{0}^{i}\quad{\rm
for}\quad i\in\\{1,\dots,N\\}$

such that

$|\hat{a}_{0}^{i}-\hat{a}_{0}^{i,b}|\lesssim 1,$ (5.12)

where

$\hat{a}_{0}^{i,b}=\frac{a_{0}^{i,b}}{(\lambda_{0}^{i,b})^{2}}.$

For all $i\in\\{1,\dots,N\\}$ and $j\in\mathbb{N}$, let us define the sequence
of perturbation parameters
$(\mathbf{a}_{j},\boldsymbol{\lambda}_{j})=(a^{1},\dots,a^{N},R^{1},\dots,R^{N})\in\ell^{\infty}(\mathbb{R}^{(n+1)N})$
by

$R_{j}^{i}=R^{i}(1+r_{j}^{i})\quad{\rm
and}\quad\frac{a_{j}^{i}}{(\lambda_{j}^{i})^{2}}=\bar{a}_{j}^{i}=\hat{a}_{0}^{i}+\tilde{a}_{j}^{i},$
(5.13)

where
$(\tilde{\boldsymbol{a}}_{j},\boldsymbol{r}_{j})=(\tilde{a}^{i},\dots,\tilde{a}^{N},r^{1},\dots
r^{N})\in\ell^{\infty}(\mathbb{R}^{(n+1)N})$ satisfy

$|r_{j}^{i}|\lesssim e^{-\tau t_{j}^{i}}\quad{\rm
and}\quad|\tilde{a}_{j}^{i}|\lesssim e^{-\tau t_{j}^{i}}\quad{\rm for}\quad
i\in\\{1,\dots,N\\}$ (5.14)

for some $\tau>0$, where $t_{j}^{i}=(1+2j)L_{i}$.

######  Definition 5.10.

Let $\sigma\in(1,+\infty)$, $n>2\sigma$ and $N\geqslant 2$. We define the so-
called perturbation map $\Upsilon_{\rm
per}:\mathbb{R}^{(n+2)N}\rightarrow\ell^{\infty}_{\tau}(\mathbb{R}^{(n+1)N})$
such that it associates balanced configurations with a sequence of admissible
perturbations. A sequence of perturbation parameters
$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\in\ell^{\infty}(\mathbb{R}^{(n+1)N})$
or
$(\boldsymbol{a}_{j},\boldsymbol{r}_{j})\in\ell^{\infty}(\mathbb{R}^{(n+1)N})$
is said to be admissible if the parameters satisfy

  1. ($\mathscr{A}_{0}$)

For $j=0$, the configuration parameters $\Upsilon_{\rm
per}^{-1}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})=(\mathbf{q},\boldsymbol{a}_{0},\boldsymbol{R})\in\mathbb{R}^{(n+2)N}$
is a balanced, that is, $(\mathbf{q},\boldsymbol{a}_{0},\boldsymbol{R})\in{\rm
Bal}_{\sigma}(\Sigma)$;

  2. ($\mathscr{A}_{1}$)

For $j\geqslant 1$, the parameters
$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\in\ell^{\infty}(\mathbb{R}^{(N+1)n})$
satisfy the set of relations (5.10),(5.11), (5.12), (5.13), and (5.14).

We denote by $(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\in{\rm
Adm}_{\sigma}(\Sigma)$ the set of admissible configurations. Notice under
(5.13), one can work indiscriminately with either parameter. In this fashion,
we call $(\boldsymbol{0},\boldsymbol{1})\in\ell^{\infty}(\mathbb{R}^{(n+1)N})$
or $(\boldsymbol{0},\boldsymbol{0})\in\ell^{\infty}(\mathbb{R}^{(n+1)N})$ the
trivial configurations.

###  5.4. Generalized Delaunay solutions

We now define a family of approximate solutions to the problem using the
Delaunay solutions from the previous section. From now on, we denote by
$\chi:\mathbb{R}\rightarrow\mathbb{R}$ the cut-off function such that

$\displaystyle\chi(x)=\begin{cases}1,&{\rm if}\ 0<|x|\leqslant\frac{1}{2}\\\
0,&{\rm if}\ \frac{1}{2}\leqslant|x|\leqslant 1\\\ \chi(x),&{\rm if}\
|x|\geqslant 1.\end{cases}$

First, one can always assume that all the balls $B_{2}(x_{i})$ are disjoint
since we may dilate the problem by some factor $\kappa>0$ that will change the
set $\Sigma$ into $\kappa\Sigma$ and a function $u$ defined in
$\mathbb{R}^{n}\setminus\Sigma$ into
$\kappa^{-\gamma_{\sigma}}u({x}{\kappa}^{-1})$ defined in
$\mathbb{R}^{n}\setminus\kappa\Sigma$.

######  Definition 5.11.

Let $\sigma\in(1,+\infty)$ and $n>2\sigma$. For any $L\gg 1$ sufficiently
large and
$\boldsymbol{L}=(L^{1},\ldots,L^{N})\in\ell^{\infty}(\mathbb{R}^{N}_{+})$ and
$\boldsymbol{R}=(R^{1},\ldots,R^{N})\in\ell^{\infty}(\mathbb{R}^{N}_{+})$ such
that (5.10) holds let
$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\in\ell^{\infty}(\mathbb{R}^{(n+1)N})$
be its associated perturbation parameters. Fix $x_{i}\in\Sigma$ for
$i\in\\{1,\dots,N\\}$, let us define the following generalized bubble tower
solution

  * (i)

$($Spherical coordinates$)$

$\widehat{U}_{(x_{i},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}(x):=\sum_{j\in\mathbb{Z}}U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})}(x),$
(5.15)

and

$\widehat{U}^{+}_{(x_{i},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}(x):=\sum_{j\in\mathbb{N}}U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})}(x),$
(5.16)

where

$U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})}(x)=\left(\frac{\lambda^{i}_{j}}{\lambda_{j}^{2}+|x-a^{i}_{j}-x_{i}|^{2}}\right)^{\gamma_{\sigma}}$
(5.17)

with

$\lambda^{i}_{j}=R_{j}^{i}e^{-2jL^{i}_{j}}\quad{\rm for}\quad j\in\mathbb{Z}.$

and

$L^{i}_{j}=L_{i}-jL_{i}+\ln R_{j}^{i}\quad{\rm for}\quad j\in\mathbb{Z}.$

  * (ii)

$($Emden–Fowler coordinates$)$

$\widehat{V}_{(x_{i},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}(t):=\sum_{j\in\mathbb{Z}}V_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})}(t),$
(5.18)

and

$\widehat{V}_{(x_{i},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}(t):=\sum_{j\in\mathbb{N}}V_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})}(t),$
(5.19)

where

$V_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})}(t)=\cosh(-\ln|x-x_{i}-a_{j}^{i}|-L^{i}_{j})^{\gamma_{\sigma}}.$
(5.20)

These will be called the general $($half$)$ bubble tower solutions.

We also have the most basic definition of this section. We observe that
although in the definition the solution is indexed by
$(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})$,
one should recall that the configuration map from Definition 5.7 relates them
and by the perturbation map from Definition 5.7, namely
$(\boldsymbol{x},\boldsymbol{L})=(\boldsymbol{a}_{j}(\boldsymbol{q},\boldsymbol{a}_{0},\boldsymbol{R}),\boldsymbol{\lambda}_{j}(\boldsymbol{q},\boldsymbol{a}_{0},\boldsymbol{R}))$
and
$(\boldsymbol{q},\boldsymbol{a}_{0},\boldsymbol{R})=(\boldsymbol{q}(\boldsymbol{x},\boldsymbol{L}),\boldsymbol{a}_{0}(\boldsymbol{x},\boldsymbol{L}),\boldsymbol{R}(\boldsymbol{x},\boldsymbol{L}))$.

######  Definition 5.12.

Let $\sigma\in(1,+\infty)$, $n>2\sigma$, and $N\geqslant 2$. For any $L\gg 1$
sufficiently large and
$\boldsymbol{L}=(L^{1},\ldots,L^{N})\in\ell^{\infty}(\mathbb{R}^{N}_{+})$ and
$\boldsymbol{R}=(R^{1},\ldots,R^{N})\in\ell^{\infty}(\mathbb{R}^{N}_{+})$ such
that (5.10) holds let
$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\in\ell^{\infty}(\mathbb{R}^{(n+1)N})$
be its associated perturbation parameters. We define its associated solution
$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda})}\in\mathcal{C}^{\infty}(\mathbb{R}^{n}\setminus\Sigma)$
as

$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}(x)=\sum_{i=1}^{N}\bar{u}_{(x_{i},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}(x).$
(5.21)

Here

$\bar{u}_{(x_{i},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}(x)=:\widehat{U}^{+}_{(x_{i},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}(x)+\chi_{i}(x)\phi_{(x_{i},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}(x),$
(5.22)

where
$\widehat{U}^{+}_{(x_{i},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}\in\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\Sigma)$
is the generalized bubble tower solution given by (5.16) and

$\phi_{(x_{i},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}(x)=\phi_{(\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}(x-x_{i})\quad{\rm
and}\quad\chi_{i}(x)=\chi(x-x_{i})\quad{\rm for\ all}\quad
i\in\\{1,\dots,N\\}$ (5.23)

with $\phi_{(\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}$ the
error function from Lemma 4.10. We say that
$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}\in\mathcal{C}^{\infty}(\mathbb{R}^{n}\setminus\Sigma)$
is an approximate solution to ($\mathcal{Q}_{2\sigma,\Sigma}$), denote by
$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}\in{\rm
Apx}_{\sigma}(\Sigma)$, whenever
$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\in{\rm Adm}_{\sigma}(\Sigma)$.
We then define the so-called perturbation map $\Upsilon_{\rm
sol}:\ell^{\infty}_{\tau}(\mathbb{R}^{(n+1)N})\rightarrow\mathcal{C}^{\infty}(\mathbb{R}^{n}\setminus\Sigma)$
such that it associates balanced configurations with sequences of admissible
perturbations.

###  5.5. Normalized approximate kernels

In this subsection, we will use the aforementioned parameters to define a
family of projections on the (normalized) approximate kernels. At least for
low Fourier eigenmodes, this family is entirely constructed by varying the
parameters in the approximate solution.

######  Definition 5.13.

Let $\sigma\in(1,+\infty)$, $n>2\sigma$ and $N\geqslant 2$. Assume that
$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\in{\rm Adm}_{\sigma}(\Sigma)$
is an admissible configuration as in Definition 5.5 with
$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}\in{\rm
Apx}_{\sigma}(\Sigma)$ their associated approximate solution as in Definition
5.12.

  1. (a)

Let us introduce some notation of normalized approximate kernels.

     * (i)

If $\ell=0$, we set

$Z_{j,0}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})=\partial_{r_{j}^{i}}U_{(x_{i},L^{i}_{j},\lambda_{j}^{i},a_{j}^{i})};$

for the zero-frequency Fourier eigenmodes.

     * (ii)

If $\ell\in\\{1,\dots,n\\}$, we set

$Z_{j,\ell}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})=\lambda_{j}^{i}\partial_{a_{j,\ell}^{i}}U_{(x_{i},L^{i}_{j},\lambda_{j}^{i},a_{j}^{i})}=-\lambda_{j}^{i}\partial_{x_{\ell}}U_{(x_{i},L^{i}_{j},\lambda_{j}^{i},a_{j}^{i})}.$

for the low-frequency Fourier eigenmodes.

We denote by
$\\{Z_{j,\ell}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\\}_{(i,j,\ell)\in\mathcal{I}_{\infty}}\subset\mathcal{C}^{0}(\mathbb{R}^{n}\setminus\Sigma)$
the family of normalized approximate kernels.

  2. (b)

Let us introduce some notation of normalized approximate cokernels.

     * (i)

If $\ell=0$, we set

$\overline{Z}_{j,0}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})=f^{\prime}_{\sigma}(U_{(x_{i},L^{i}_{j},\lambda_{j}^{i},a_{j}^{i})})Z_{j,0}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j});$

for the zero-frequency Fourier eigenmodes.

     * (ii)

If $\ell\in\\{1,\dots,n\\}$, we set

$\overline{Z}_{j,\ell}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})=f^{\prime}_{\sigma}(U_{(x_{i},L^{i}_{j},\lambda_{j}^{i},a_{j}^{i})})Z_{j,\ell}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})$

for the low-frequency Fourier eigenmodes.

We denote by
$\\{\overline{Z}_{j,\ell}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\\}_{(i,j,\ell)\in\mathcal{I}_{\infty}}\subset\mathcal{C}^{0}(\mathbb{R}^{n}\setminus\Sigma)$
the family of normalized approximate cokernels.

These normalized kernels satisfy some orthogonality conditions, which will be
important in applying a finite-dimensional reduction.

######  Lemma 5.14.

Let $\sigma\in(1,+\infty)$, $n>2\sigma$ and $N\geqslant 2$. Assume that
$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\in{\rm Adm}_{\sigma}(\Sigma)$
is an admissible configuration as in Definition 5.5 with
$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}\in{\rm
Apx}_{\sigma}(\Sigma)$ their associated approximate solution as in Definition
5.12. Then, one has

  * (i)

If $\ell\in\\{1,\dots,n\\}$, then

$\displaystyle\int_{\mathbb{R}^{n}}\overline{Z}_{j,\ell}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\overline{Z}_{j^{\prime},\ell^{\prime}}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\mathrm{d}x=\frac{4(n-2\sigma)^{2}}{n}\left(\delta_{\ell,\ell^{\prime}}+\mathrm{o}(1)\right)e^{-(\gamma_{\sigma}+1)|t_{j}^{i}-t_{j^{\prime}}^{i}|},$
(5.24)

where $\delta_{\ell,\ell^{\prime}}$ is Kronecker’s delta;

  * (ii)

If $\ell=0$, then

$\displaystyle\int_{\mathbb{R}^{n}}\overline{Z}_{j,0}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})Z_{j^{\prime},0}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\mathrm{d}x=C_{0}(1+\mathrm{o}(1))e^{-\gamma_{\sigma}|t_{j}^{i}-t_{j^{\prime}}^{i}|}$
(5.25)

for some $C_{0}>0$.

###### Proof.

Initially, let us observe that by Lemma 4.6, the set of bounded solutions to

$\phi-(-\Delta)^{-\sigma}(f^{\prime}_{\sigma}(U_{(x_{i},L^{i}_{j},\lambda_{j}^{i},a_{j}^{i})})\phi)=0\quad{\rm
in}\quad\mathbb{R}^{n}$

is spanned by
$\\{\overline{Z}_{j,0}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j}),\dots,\overline{Z}_{j,n}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\\}$
for any $i\in\\{1,\dots,N\\}$ and $j\in\mathbb{N}$.

Without loss of generality, assume in the following that $x_{i}=0$. For
$\ell=0$, we will repeatedly use the following estimates

$\left|Z_{j,0}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})(x)\right|\lesssim\begin{cases}|x|^{-\gamma_{\sigma}}V_{(x_{i},L_{i},a_{j}^{i},\lambda^{i}_{j})}(-\ln|x|),&{\rm
if}\ |x|\leqslant 1,\\\
|x|^{-2\gamma_{\sigma}}(\lambda_{j}^{i})^{\gamma_{\sigma}},&{\rm if}\
|x|\geqslant 1.\end{cases}$ (5.26)

In addition, we have also have

$Z_{j,\ell}^{i}(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})(x)=2\gamma_{\sigma}V_{a_{j}^{i},L_{i},t_{i}}(-\ln|x|)^{\frac{\gamma_{\sigma}+1}{\gamma_{\sigma}}}|x-a_{j}^{i}-x_{i}|^{-\gamma_{\sigma}-1}\left(x-a_{j}^{i}-x_{i}\right)_{\ell},$

Then, after recentering at $x_{i}=0$, it is easy to see that the following
orthogonality conditions (5.24) are in force. Similar estimates also hold true
for $\ell=0$; the orthogonality condition in (5.25) is also satisfied.

The lemma is then proved. ∎

###  5.6. Weighted functional spaces

It is convenient to define the suitable function spaces on which we will run
our perturbation technique.

######  Definition 5.15.

Let $\alpha\in(0,1)$ and $\zeta_{1},\zeta_{2}\in\mathbb{R}$ such that
$\zeta_{1}<0$ and $\zeta_{2}>0$. We set the weighted norm

$\|u\|_{\mathcal{C}^{\alpha}_{\zeta_{1},\zeta_{2}}(\mathbb{R}^{n}\setminus\Sigma)}=\|\operatorname{dist}(x,\Sigma)^{-\zeta_{1}}u\|_{\mathcal{C}^{\alpha}(B_{1}(\Sigma))}+\||x|^{-\zeta_{2}}u\|_{\mathcal{C}^{\alpha}(\mathbb{R}^{n}\setminus
B_{1}(\Sigma))}.$

In other words, one that
$u\in\mathcal{C}^{\alpha}_{\zeta_{1},\zeta_{2}}(\mathbb{R}^{n}\setminus\Sigma)$
if and only if

  * (i)

$($Near the singular set$)$ it is bounded by a constant times
$|x-x_{i}|^{\zeta_{1}}$ and has its $\ell$-th order partial derivatives
bounded by a constant times $|x-x_{i}|^{\zeta_{1}-\ell}$ for
$\ell\leqslant\alpha$ near each singular point $x_{i}\in\Sigma$.

  * (ii)

$($Away from the singular set$)$ it is bounded by $|x|^{\zeta_{2}}$ and has
its $\ell$-th order partial derivatives bounded by a constant times
$|x|^{\zeta_{2}-\ell}$ for $\ell\leqslant\alpha$.

Note that we are implicitly assuming that $0\in\Sigma$, in order to simplify
the notation.

######  Definition 5.16.

Let $\sigma\in(1,+\infty)$, $n>2\sigma$ and $N\geqslant 2$. We define the
following weighted norms

$\|u\|_{\mathcal{C}_{*,\tau}(\mathbb{R}^{n}\setminus\Sigma)}=\|u\|_{\mathcal{C}_{\min\\{\zeta_{1},-\gamma_{\sigma}+\tau\\},-n-2\sigma}(\mathbb{R})}$
(5.27)

and

$\|h\|_{\mathcal{C}_{**,\tau}(\mathbb{R}^{n}\setminus\Sigma)}=\|h\|_{\mathcal{C}_{n+\tau,-n+2\sigma}(\mathbb{R})},$
(5.28)

where

$-\gamma_{\sigma}<\zeta_{1}<\min\left\\{-\gamma_{\sigma}+2\sigma,0\right\\}.$
(5.29)

Here $0<\tau\ll 1$ small enough is given in the definition of the perturbation
parameters (5.13) and (5.14). In this fashion, we denote by
$\mathcal{C}_{*,\tau}(\mathbb{R}^{n}\setminus\Sigma)$ and
$\mathcal{C}_{**,\tau}(\mathbb{R}^{n}\setminus\Sigma)$ the corresponding
weighted Hölder spaces.

Let us make some observations regarding the last definition.

######  Remark 5.17.

We emphasize that to simplify the notation, many times we will ignore the
small perturbation and just the weight near the singular set as dist
$(x,\Sigma)^{-\zeta_{1}},\operatorname{dist}(x,\Sigma)^{2\sigma-\zeta_{1}}$,
respectively. The weights in Definition 5.16 are suitably chosen to guarantee
the invertibility and Fredholmness of the linearized operator around
approximate solutions on weighted Hölder spaces; this will be clear in the
reduction method we apply in the remaining subsections.

###  5.7. Perturbation of the approximate solution

This subsection is devoted to performing a perturbation method based on the
approximated solution, which requires linearizing
(${\mathcal{Q}}^{\prime}_{2\sigma,\Sigma}$) around the approximate solution
(5.21) and estimating both the weighted norm of the right-inverse for the
linearized operator given by and the associated remainder error. We emphasize
that the balancing formulas and the orthogonality conditions for the
normalized kernels discussed above will be building blocks of our
construction.

Let us explain our strategy in more detail. First, we consider the nonlinear
operator defined
$\mathscr{N}_{\sigma}:\mathcal{C}^{0}(\mathbb{R}^{n}\setminus\Sigma)\rightarrow\mathcal{C}^{2\sigma}(\mathbb{R})$
given by

$\mathscr{N}_{\sigma}({u})={u}-(-\Delta)^{-\sigma}(f_{\sigma}\circ{u}).$
(5.30)

Notice that (${\mathcal{Q}}^{\prime}_{2\sigma,\Sigma}$) can be reformulated as

$\mathscr{N}_{\sigma}({u})=0\quad{\rm in}\quad\mathbb{R}^{n}\setminus\Sigma.$

Next, by linearizing this operator around the approximate solution, we find a
linear operator
$\mathscr{L}_{\sigma}[\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}]:\mathcal{C}^{0}(\mathbb{R}^{n}\setminus\Sigma)\rightarrow\mathcal{C}^{2\sigma}(\mathbb{R}^{n}\setminus\Sigma)$
given by

$\mathscr{L}_{\sigma}[\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}](\phi)=\phi-(-\Delta)^{-\sigma}(f^{\prime}_{\sigma}\circ\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})})\phi.$
(5.31)

For the sake of simplicity, let us denote

$\mathscr{L}_{\sigma}[\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}]:=\mathscr{L}_{\sigma}(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j}).$

####  5.7.1. Quantitative estimates

Our first estimate concerns the nonlinear operator defined as (5.30) applied
to the approximate solution
$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}\in\mathcal{C}^{\infty}(\mathbb{R}^{n}\setminus\Sigma)$
given by (5.21), namely

$\displaystyle\mathcal{N}_{\sigma}(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})$
$\displaystyle:=\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})})$
$\displaystyle=\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}-(-\Delta)^{-\sigma}(f_{\sigma}\circ\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}).$

We emphasize that we must suitably choose the weighted norm in (5.28) so that
our following estimates have the correct decay.

######  Lemma 5.18.

Let $\sigma\in(1,+\infty)$, $n>2\sigma$ and $N\geqslant 2$. Assume that
$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\in{\rm Adm}_{\sigma}(\Sigma)$
is an admissible configuration as in Definition 5.5 with
$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})}\in{\rm
Apx}_{\sigma}(\Sigma)$ their associated approximate solution as in Definition
5.12. Then, there exists a weight $\zeta_{1}<0$ satisfying (5.29) such that

$\|\mathcal{N}_{\sigma}(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\|_{\mathcal{C}_{**,\tau}(\mathbb{R}^{n}\setminus\Sigma)}\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}$ (5.32)

for some $\xi>0$ uniformly on $L\gg 1$ large.

###### Proof.

For the sake of simplicity, we shall prove the estimate in (5.32) for the
$L^{\infty}-$norm. Namely, we need to quantitatively estimate the term
$|\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})|$
and then a applying a classical perturbation technique.

The rest of the proof will be divided into two cases.

Case 1:
$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})=(\boldsymbol{0},\boldsymbol{1})$
for all ${j\in\mathbb{N}}$.

In this case, the approximate solution
$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}\in\mathcal{C}^{\infty}(\mathbb{R}^{n}\setminus\Sigma)$
is given by

$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}(x)=\sum_{i=1}^{N}\left[\sum_{j\in\mathbb{N}}U_{(x_{i},L_{j}^{i},\boldsymbol{0},\lambda_{j}^{i})}(x)+\chi_{i}(x)\phi_{i}(x-x_{i})\right],$

where

$U_{(x_{i},L_{j}^{i},\boldsymbol{0},\lambda_{j}^{i})}(x):=\left(\frac{\lambda_{j}^{i}}{|\lambda_{j}^{i}|^{2}+|x-x_{i}|^{2}}\right)^{\gamma_{\sigma}}.$

Without loss of generality, assume $x_{1}=0$. Before we prove the estimate of
$|\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})|$,
we first prove the following claim:

Claim 1: The following estimate holds

$\displaystyle|\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})|\lesssim\begin{cases}|x-x_{i}|^{\zeta_{1}-2\sigma}e^{-\gamma_{\sigma}L(1+\xi)},\
&{\rm if}\ 0<\mathrm{d}(x,\Sigma)<\frac{1}{2},\\\
e^{-\gamma_{\sigma}L(1+\xi)},\ &{\rm if}\
\frac{1}{2}\leqslant\mathrm{d}(x,\Sigma)<1,\\\
|x|^{-(n+2\sigma)}e^{-\gamma_{\sigma}L(1+\xi)},\ &{\rm if}\
\mathrm{d}(x,\Sigma)\geqslant 1,\end{cases}$

where

$\displaystyle\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}):=\sum_{i=1}^{N}\sum_{j\in\mathbb{N}}f_{\sigma}(U_{(x_{i},L_{j}^{i},\boldsymbol{0},\lambda_{j}^{i})})-f_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}).$
(5.33)

As a matter of fact, we proceed by a direct estimate in terms in the
asymptotic behavior of the bubble tower solution. Without loss of generality,
assume $x_{1}=0$. The proof will be divided into three steps: the exterior,
transition, and interior region, respectively.

Step 1: If $\mathrm{d}(x,\Sigma)\geqslant 1$, then
$|\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})|\lesssim|x|^{-(n+2\sigma)}e^{-\gamma_{\sigma}L(1+\xi)}$.

In this region, we notice $\chi_{i}(x)=0$ for all $i\in\\{2,\dots,N\\}$ when
$\mathrm{d}(x,\Sigma)\geqslant 1$. Next, using that

$|U_{(x_{i},L_{j}^{i},\boldsymbol{0},\lambda_{j}^{i})}(x)|\sim\left(\lambda_{j}^{i}\right)^{\gamma_{\sigma}}|x|^{-(n-2\sigma)}\quad{\rm
as}\quad|x|\rightarrow+\infty$

and recalling the relation in (5.10), we have

$\displaystyle\left|\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})\right|$
$\displaystyle=c_{n,\sigma}\left|\left(\sum_{i=1}^{N}\sum_{j\in\mathbb{N}}U_{(x_{i},L_{j}^{i},\boldsymbol{0},\lambda_{j}^{i})}\right)^{\frac{n+2\sigma}{n-2\sigma}}-\sum_{i=1}^{N}\sum_{j\in\mathbb{N}}\left(U_{(x_{i},L_{j}^{i},\boldsymbol{0},\lambda_{j}^{i})}\right)^{\frac{n+2\sigma}{n-2\sigma}}\right|$
$\displaystyle\lesssim\left(e^{-\frac{(n-2\sigma)L}{2}}|x|^{-(n-2\sigma)}\right)^{\frac{n+2\sigma}{n-2\sigma}}$
$\displaystyle\lesssim e^{-\frac{(n+2\sigma)L}{2}}|x|^{-(n+2\sigma)},$

which finishes the proof of the first step.

Step 2: If $\frac{1}{2}\leqslant|x|\leqslant 1$, then
$|\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})|\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}$ for some $\xi>0$.

In this case, it is easy to verify the estimate

$|\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})|\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}$

for some $\xi>0$.

Step 3: If $0<|x|\leqslant\frac{1}{2}$, then
$|\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})|\lesssim|x-x_{i}|^{\zeta_{1}-2\sigma}e^{-\gamma_{\sigma}L(1+\xi)}$
for some $\xi>0$.

Notice that $\chi_{1}(x)\equiv 1$ and $\chi_{i}(x)0$ for
$i\in\\{2,\dots,N\\}$. By definition, it follows that

$\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}={\widehat{U}}_{(x_{1},L_{j}^{1},0,\lambda_{j}^{1})}-\left(1-\chi_{1}\right)\phi_{1}+\sum_{i=2}^{N}\left(\sum_{j\in\mathbb{N}}U_{(x_{i},L_{j}^{i},\boldsymbol{0},\lambda_{j}^{i})}+\chi_{i}\phi_{i}\right)-\sum_{j\in\mathbb{Z}\setminus\mathbb{N}}U_{(x_{1},L_{j}^{1},0,\lambda_{j}^{1})}.$

Hence, by an easy computation, we obtain

$\displaystyle\left|\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})\right|$
$\displaystyle\lesssim\left|\sum_{j\in\mathbb{N}}(U_{(x_{1},L_{j}^{1},0,\lambda_{j}^{1})})^{\frac{n+2\sigma}{n-2\sigma}}-\left(\sum_{j\in\mathbb{N}}U_{(x_{1},L_{j}^{1},0,\lambda_{j}^{1})}+\mathcal{O}(e^{-\gamma_{\sigma}(1+\xi)})\right)^{\frac{n+2\sigma}{n-2\sigma}}\right|+\mathcal{O}(e^{-\gamma_{\sigma}(1+\xi)})$
$\displaystyle\lesssim\sum_{j\in\mathbb{N}}(U_{(x_{1},L_{j}^{1},0,\lambda_{j}^{1})})^{\frac{4\sigma}{n-2\sigma}}e^{-\gamma_{\sigma}}+\mathcal{O}(e^{-\gamma_{\sigma}(1+\xi)}))$
(5.34)
$\displaystyle\lesssim|x|^{-2\sigma}\sum_{j\in\mathbb{N}}V_{(x_{1},L_{j}^{1},0,\lambda_{j}^{1})}^{\frac{4\sigma}{n-2\sigma}}e^{-\gamma_{\sigma}L}+\mathcal{O}(e^{-\gamma_{\sigma}(1+\xi)}),$

where

$V_{(x_{1},L_{j}^{1},0,\lambda_{j}^{1})}(-\ln|x|):=V_{\rm
sph}(-\ln|x|-L_{1}-2jL_{1})$

and we recall that the spherical solution $v_{\rm sph}$ is defined as (4.13).

Furthermore, it is straightforward to see that when
$0<|x|\leqslant\frac{1}{2}$ there exists $\xi>0$ and $\zeta_{1}<0$ satisfying

$|x|^{-\zeta_{1}}\left(\sum_{j\in\mathbb{Z}}V_{(x_{1},L_{j}^{1},0,\lambda_{j}^{1})}(-\ln|x|)\right)^{\frac{4\sigma}{n-2\sigma}}\lesssim
e^{-\xi L_{1}}.$ (5.35)

Indeed, if $-\infty<t<{L_{1}}$ there exists $C_{1}>0$ such that $|x|\leqslant
C_{1}$ and

$\sum_{j\in\mathbb{Z}}V_{(x_{1},L_{j}^{1},0,\lambda_{j}^{1})}(-\ln|x|)\leqslant
C_{1}e^{-\gamma_{\sigma}L_{1}}.$ (5.36)

Also, if ${L_{1}}\leqslant t<+\infty$, there exists $C_{2}>0$ such that
$|x|\leqslant C_{2}e^{-{L_{1}}/{2}}$ and

$\sum_{j\in\mathbb{Z}}V_{(x_{1},L_{j}^{1},0,\lambda_{j}^{1})}(-\ln|x|)\leqslant
C_{2}.$ (5.37)

Finally, combining (5.7.1) and (5.35) implies

$\displaystyle\left|\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})\right|\lesssim|x|^{\zeta_{1}-2\sigma}e^{-\gamma_{\sigma}(1+\xi)},$

which finishes the proof of the first claim.

We now proceed to the proof of our preliminary estimate.

Claim 2: The following estimates holds

$\displaystyle|\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})})|\lesssim\begin{cases}|x-x_{i}|^{\min\\{\zeta_{1}-\tau,-\gamma_{\sigma}+\tau\\}}e^{-\gamma_{\sigma}L(1+\xi)},\
&{\rm if}\ 0<\mathrm{d}(x,\Sigma)<\frac{1}{2},\\\
e^{-\gamma_{\sigma}L(1+\xi)},\ &{\rm if}\
\frac{1}{2}\leqslant\mathrm{d}(x,\Sigma)<1,\\\
|x|^{2\sigma-n}e^{-\gamma_{\sigma}L(1+\xi)},\ &{\rm if}\
\mathrm{d}(x,\Sigma)\geqslant 1.\end{cases}$

As before, the proof will be divided into three steps as follows:

Step 1: If $\mathrm{d}(x,\Sigma)\geqslant 1$, then
$|\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})|\lesssim|x|^{2\sigma-n}e^{-\gamma_{\sigma}L(1+\xi)}$.

Notice that $\chi_{i}(x)\equiv 0$ for all $i\in\\{1,\dots,N\\}$ and
$x\in\mathbb{R}^{n}\setminus\Sigma$ such that $\mathrm{d}(x,\Sigma)\geqslant
1$. From this, we get

$\displaystyle\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})$
$\displaystyle=\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}-(-\Delta)^{-\sigma}(f_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}))$
$\displaystyle=\sum_{i=1}^{N}\sum_{j\in\mathbb{N}}U_{(x_{i},L_{j}^{i},\boldsymbol{0},\lambda_{j}^{i})}(x)-(-\Delta)^{-\sigma}f_{\sigma}\left(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}\right)$
$\displaystyle=\int_{\mathbb{R}^{n}}\left(\sum_{i=1}^{N}\sum_{j\in\mathbb{N}}f_{\sigma}(U_{(x_{i},L_{j}^{i},\boldsymbol{0},\lambda_{j}^{i})}(y))-f_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}(y))\right)\mathcal{R}_{\sigma}(x-y)\mathrm{d}y$
$\displaystyle=\int_{\mathbb{R}^{n}}\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}(y))\mathcal{R}_{\sigma}(x-y)\mathrm{d}y$
$\displaystyle=\left(\int_{|y|\leqslant
1}+\int_{1\leqslant|y|\leqslant|x|}+\int_{|y|\geqslant|x|}\right)\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}(y))\mathcal{R}_{\sigma}(x-y)\mathrm{d}y$
$\displaystyle=:I_{11}+I_{12}+I_{13}.$

where we recall that
$\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})$
is given by (5.33).

Applying Step 1 of Claim 1, we have

$|I_{11}|\lesssim e^{-\gamma_{\sigma}(1+\xi)}\int_{|y|\leqslant
1}{|x-y|^{2\sigma-n}}{|y|^{\zeta_{1}-2\sigma}}\mathrm{d}y\lesssim|x|^{2\sigma-n}e^{-\gamma_{\sigma}(1+\xi)},$
$|I_{12}|\lesssim
e^{-\gamma_{\sigma}(1+\xi)}\int_{1\leqslant|y|\leqslant|x|}{|x-y|^{n-2\sigma}}{|y|^{-(n+2\sigma)}}\mathrm{d}y\lesssim|x|^{2\sigma-n}e^{-\gamma_{\sigma}L(1+\xi)},$

and

$|I_{13}|\lesssim
e^{-\gamma_{\sigma}(1+\xi)}\int_{|y|\geqslant|x|}{|x-y|^{2\sigma-n}}{|y|^{-(n+2\sigma)}}\mathrm{d}y\lesssim|x|^{-n}e^{-\gamma_{\sigma}(1+\xi)}.$

Combining the above estimates, we finish the proof of Step 1.

Step 2: If $\frac{1}{2}\leqslant|x|\leqslant 1$, then
$|\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})})|\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}$ for some $\xi>0$.

In this case, it holds

$\displaystyle\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})$
$\displaystyle=\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}-(-\Delta)^{-\sigma}(f_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}))$
$\displaystyle=\sum_{i=1}^{N}\sum_{j\in\mathbb{N}}U_{(x_{i},L_{j}^{i},\boldsymbol{0},\lambda_{j}^{i})}+\chi_{1}\phi_{1}-(-\Delta)^{-\sigma}f_{\sigma}\left(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}\right)$
$\displaystyle=\int_{\mathbb{R}^{n}}\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})(y)\mathcal{R}_{\sigma}(x-y)\mathrm{d}y+e^{-\gamma_{\sigma}L(1+\xi)}$
$\displaystyle=\left(\int_{|y|\leqslant\frac{|x|}{2}}+\int_{\frac{|x|}{2}\leqslant|y|\leqslant
1}+\int_{1\leqslant|y|\leqslant 2|x|}+\int_{|y|\geqslant
2|x|}\right)\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})(y)\mathcal{R}_{\sigma}(x-y)\mathrm{d}y+\mathcal{O}(e^{-\gamma_{\sigma}L(1+\xi)})$
$\displaystyle=:I_{21}+I_{22}+I_{23}+I_{24}+\mathcal{O}(e^{-\gamma_{\sigma}L(1+\xi)}).$

Applying Step 2 of Claim 1, we get

$\displaystyle|I_{21}|\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}\int_{|y|\leqslant\frac{|x|}{2}}{|x-y|^{2\sigma-n}}{|y|^{\zeta_{1}-2\sigma}}\mathrm{d}y\lesssim|x|^{\zeta_{1}}e^{-\gamma_{\sigma}L(1+\xi)}\lesssim
e^{-\gamma_{\sigma}L(1+\xi)},$ $\displaystyle|I_{22}|\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}\int_{\frac{|x|}{2}\leqslant|y|\leqslant
1}{|x-y|^{2\sigma-n}}{|y|^{\zeta_{1}-2\sigma}}\mathrm{d}y\lesssim|x|^{\zeta_{1}-2\sigma}e^{-\gamma_{\sigma}L(1+\xi)}\int_{|y-x|\leqslant\frac{3}{2}}{|x-y|^{2\sigma-n}}\mathrm{d}y\lesssim
e^{-\gamma_{\sigma}L(1+\xi)},$ $\displaystyle|I_{23}|\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}\int_{1\leqslant|y|\leqslant
2|x|}{|x-y|^{2\sigma-n}}{|y|^{-(n+2\sigma)}}\mathrm{d}y\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}\int_{|y-x|\leqslant
3|x|}{|x-y|^{2\sigma-n}}\mathrm{d}y\lesssim e^{-\gamma_{\sigma}L(1+\xi)},$

and

$\displaystyle|I_{24}|\lesssim e^{-\gamma_{\sigma}L(1+\xi)}\int_{|y|\geqslant
2|x|}{|x-y|^{2\sigma-n}}{|y|^{-(n+2\sigma)}}\mathrm{d}y\lesssim|x|^{2\sigma-n}e^{-\gamma_{\sigma}L(1+\xi)}\int_{|y|\geqslant
1}{|y|^{-(n+2\sigma)}}\mathrm{d}y\lesssim e^{-\gamma_{\sigma}L(1+\xi)}.$

Consequently, the proof of Step 2 follows.

Step 3: If $0<|x|\leqslant\frac{1}{2}$, then
$|\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})|\lesssim|x|^{\zeta_{1}-\tau}e^{-\gamma_{\sigma}L(1+\xi)}$.

Similarly to the previous steps, we obtain

$\displaystyle\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})$
$\displaystyle=\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}-(-\Delta)^{-\sigma}(f_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}))$
$\displaystyle=\sum_{i=1}^{N}\sum_{j\in\mathbb{N}}U_{(x_{i},L_{j}^{i},\boldsymbol{0},\lambda_{j}^{i})}+\phi_{1}-(-\Delta)^{-\sigma}f_{\sigma}\left(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}\right)$
$\displaystyle=\int_{\mathbb{R}^{n}}\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}(y))\mathcal{R}_{\sigma}(x-y)\mathrm{d}y+e^{-\gamma_{\sigma}L(1+\xi)}$
$\displaystyle=\left(\int_{|y|\leqslant\frac{|x|}{2}}+\int_{\frac{|x|}{2}\leqslant|y|\leqslant
2|x|}+\int_{2|x|\leqslant|y|\leqslant 1}+\int_{|y|\geqslant
1}\right)\mathcal{D}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}(y))\mathcal{R}_{\sigma}(x-y)\mathrm{d}y+\mathcal{O}(e^{-\gamma_{\sigma}L(1+\xi)})$
$\displaystyle=:I_{31}+I_{32}+I_{33}+I_{34}+\mathcal{O}(e^{-\gamma_{\sigma}L(1+\xi)}).$

Applying Step 3 of Claim 1, we get

$\displaystyle|I_{31}|\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}\int_{|y|\leqslant\frac{|x|}{2}}{|x-y|^{2\sigma-n}}{|y|^{\zeta_{1}-2\sigma}}\mathrm{d}y\lesssim|x|^{\zeta_{1}}e^{-\gamma_{\sigma}L(1+\xi)},$
$\displaystyle|I_{32}|\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}\int_{\frac{|x|}{2}\leqslant|y|\leqslant
2|x|}{|x-y|^{2\sigma-n}}{|y|^{\zeta_{1}-2\sigma}}\mathrm{d}y\lesssim|x|^{\zeta_{1}}e^{-\gamma_{\sigma}L(1+\xi)},$
$\displaystyle|I_{33}|\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}\int_{2|x|\leqslant|y|\leqslant
1}{|x-y|^{2\sigma-n}}{|y|^{\zeta_{1}-2\sigma}}\mathrm{d}y\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}|x|^{\zeta_{1}}\int_{2|x|\leqslant|y|\leqslant
1}{|y|^{-n}}\mathrm{d}y\lesssim|x|^{\zeta_{1}-\tau}e^{-\gamma_{\sigma}L(1+\xi)},$

and

$\displaystyle|I_{34}|\lesssim e^{-\gamma_{\sigma}L(1+\xi)}\int_{|y|\geqslant
1}{|x-y|^{2\sigma-n}}{|y|^{-(n+2\sigma)}}\mathrm{d}y\lesssim
e^{-\gamma_{\sigma}L(1+\xi)}\int_{|y|\geqslant
1}{|y|^{-(n+2\sigma)}}\mathrm{d}y\lesssim e^{-\gamma_{\sigma}L(1+\xi)}.$

Therefore, for $|x|\leqslant\frac{1}{2}$, we conclude

$|\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})|\lesssim|x|^{\zeta_{1}-\tau}e^{-\gamma_{\sigma}L(1+\xi)},$

which gives us the desired estimate in Step 3.

By combining the last three steps, the proof of the first case is concluded.

Now we consider the case of a general configuration. We will use a
perturbation technique based on the last case in this situation.

Case 2:
$(\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})\neq(\boldsymbol{0},\boldsymbol{1})$
for some ${j\in\mathbb{N}}$.

Initially, we will prove the following decomposition.

Claim 3: It holds that

$\displaystyle|\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{a}_{j},\boldsymbol{\lambda}_{j})})-\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})|$
$\displaystyle\leqslant\sum_{i=1}^{N}\sum_{j\in\mathbb{N}}(-\Delta)^{-\sigma}\left[|f^{\prime}_{\sigma}(U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})})-f^{\prime}_{\sigma}(U_{(x_{i},L_{j}^{i},\boldsymbol{0},\lambda_{j}^{i})})|\left(|\partial_{r_{j}^{i}}U_{(x_{\ell},L_{\ell}^{i},a_{\ell}^{i},\lambda_{\ell}^{i})}||r_{j}^{i}|+\sum_{\ell=1}^{n}|\partial_{a_{j,\ell}^{i}}U_{(x_{\ell},L_{\ell}^{i},a_{\ell}^{i},\lambda_{\ell}^{i})}||a_{j,\ell}^{i}|\right)\right]$
$\displaystyle=\sum_{i=1}^{N}\sum_{j\in\mathbb{N}}(-\Delta)^{-\sigma}\left[|\partial_{r_{j}^{i}}\widetilde{\mathcal{D}}^{\prime}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})||r_{j}^{i}|+\sum_{\ell=1}^{n}|\partial_{a_{j,\ell}^{i}}\tilde{\mathcal{D}}^{\prime}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})||a_{j,\ell}^{i}|\right],$

where

$\displaystyle\widetilde{\mathcal{D}}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}):=\sum_{i=1}^{N}\sum_{j\in\mathbb{N}}f_{\sigma}(U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})})-f_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}).$
(5.38)

To prove this fact, we will differentiate
$\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})$
with respect to the parameters $r_{j}^{i},a^{i}_{j,\ell}$. Since the variation
is linear in the displacements of the parameters, we vary the parameter of one
point at one time. First, with respect to $r_{j}^{i}$, we have

$\displaystyle{\partial_{r_{j}^{i}}}\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})$
$\displaystyle=\partial_{r_{j}^{i}}U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})}-(-\Delta)^{-\sigma}(f_{\sigma}^{\prime}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}){\partial_{r_{j}^{i}}}U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})})$
$\displaystyle=(-\Delta)^{-\sigma}[(f^{\prime}_{\sigma}(U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})})-f^{\prime}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}){\partial_{r_{j}^{i}}}U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})}].$

Second, with respect to $a_{j,\ell}^{i}$, we obtain

$\displaystyle{\partial_{a_{j,\ell}^{i}}}\mathscr{N}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})$
$\displaystyle=(-\Delta)^{-\sigma}\left[(f^{\prime}_{\sigma}(U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})})-f^{\prime}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})}))\sum_{\ell=1}^{n}{\partial_{a_{j,\ell}^{i}}}U_{(x_{\ell},L_{\ell}^{i},a_{\ell}^{i},\lambda_{\ell}^{i})}\right].$

This fact concludes the proof of Claim 3.

Next, we shall obtain $L^{\infty}$-estimate in the sense below. We first
consider the case of the parameters $r_{j}^{i}$.

Claim 4: The following estimate holds

$\displaystyle\left|\partial_{r_{j}^{i}}\widetilde{\mathcal{D}}^{\prime}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})\right|\lesssim\begin{cases}\mathrm{d}(x,\Sigma)^{\min\\{\zeta_{1},-\gamma_{\sigma}+\tau\\}-2\sigma}e^{-\gamma_{\sigma}L(1+\xi)},\
&{\rm if}\ 0<\mathrm{d}(x,\Sigma)<1,\\\
|x|^{-(n+2\sigma)}e^{-\gamma_{\sigma}L(1+\xi)},\ &{\rm if}\
\mathrm{d}(x,\Sigma)\geqslant 1.\end{cases}$

As before, we consider two cases separately.

Step 1: If $\mathrm{d}(x,\Sigma)\geqslant 1$, then

$[f^{\prime}_{\sigma}(U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})})-f^{\prime}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})]{\partial_{r_{j}^{i}}}U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})}\lesssim|x|^{-(n+2\sigma)}e^{-\gamma_{\sigma}L(1+\xi)}e^{-\nu
t_{j}^{i}},$

for a suitable choice of $\nu>0$.

As a matter of fact, we have

$\displaystyle[f^{\prime}_{\sigma}(U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})})-f^{\prime}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})]{\partial_{r_{j}^{i}}}U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})}$
$\displaystyle\lesssim\left(e^{-\gamma_{\sigma}L}|x|^{-(n-2\sigma)}\right)^{\frac{4n}{n-2\sigma}}e^{-\gamma_{\sigma}L(2j+1)}|x|^{-(n-2\sigma)}$
$\displaystyle\lesssim|x|^{-(n+2\sigma)}e^{-\gamma_{\sigma}L(1+\xi)}e^{-\nu
t_{j}^{i}},$

which by (5.13) and (5.14) concludes the proof of this step.

Step 2: If $0<\mathrm{d}(x,\Sigma)<1$, then

$\displaystyle[f^{\prime}_{\sigma}(U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})})-f^{\prime}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})]{\partial_{r_{j}^{i}}}U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})}\lesssim\mathrm{d}(x,\Sigma)^{\min\\{-\gamma_{\sigma}+\tau,\zeta_{1}\\}-2\sigma}e^{-\gamma_{\sigma}L(1+\xi)}$

for some $-\gamma_{\sigma}<\zeta_{1}<\min\\{0,-\gamma_{\sigma}+2\sigma\\}$ and
$0<\tau\ll 1$ small enough.

In this situation, we may assume without loss of generality that
$|x-x_{i}|\leqslant 1$ for $i\in\\{2,\dots,N\\}$. Hence, we proceed similarly
to the proof of the estimates (5.36) and (5.37) to find

$\displaystyle[f^{\prime}_{\sigma}(U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})})-f^{\prime}_{\sigma}(\bar{u}_{(\boldsymbol{x},\boldsymbol{L},\boldsymbol{0},\boldsymbol{1})})]{\partial_{r_{j}^{i}}}U_{(x_{i},L_{j}^{i},a_{j}^{i},\lambda_{j}^{i})}$
$\displaystyle\lesssim|x-x_{i}|^{-2\sigma}\left(\sum_{j\in\mathbb{N}}V_{(x_{i^{\prime}},L_{j^{\prime}}^{i^{\prime}}\lambda_{j^{\prime}}^{i^{\prime}},a_{j^{\prime}}^{i^{\prime}})}(-\ln|x-x_{i}|)\right)^{\frac{n+2\sigma}{n-2\sigma}}e^{-\gamma_{\sigma}L(2j+1)}$
$\displaystyle\lesssim|x-x_{i}|^{\zeta_{1}-2\sigma}e^{-\gamma_{\sigma}L(1+\xi)}e^{-\nu
t_{j}^{i}}.$