11institutetext: The Hong Kong University of Science and Technology, Hong
Kong, China
11email<EMAIL_ADDRESS>
22institutetext: HKUST Shenzhen Research Institute, Shenzhen, China
33institutetext: Southern University of Science and Technology, Shenzhen,
China

# DiffCMR: Fast Cardiac MRI Reconstruction with Diffusion Probabilistic Models

Tianqi Xiang 1222These authors contributed equally.1222These authors
contributed equally. Wenjun Yue 112222These authors contributed
equally.2222These authors contributed equally. Yiqun Lin 11 Jiewen Yang 11
Zhenkun Wang 33 Xiaomeng Li 1111Corresponding author.1111Corresponding author.

###### Abstract

Performing magnetic resonance imaging (MRI) reconstruction from under-sampled
k-space data can accelerate the procedure to acquire MRI scans and reduce
patients’ discomfort. The reconstruction problem is usually formulated as a
denoising task that removes the noise in under-sampled MRI image slices.
Although previous GAN-based methods have achieved good performance in image
denoising, they are difficult to train and require careful tuning of
hyperparameters. In this paper, we propose a novel MRI denoising framework
DiffCMR by leveraging conditional denoising diffusion probabilistic models.
Specifically, DiffCMR perceives conditioning signals from the under-sampled
MRI image slice and generates its corresponding fully-sampled MRI image slice.
During inference, we adopt a multi-round ensembling strategy to stabilize the
performance. We validate DiffCMR with cine reconstruction and T1/T2 mapping
tasks on MICCAI 2023 Cardiac MRI Reconstruction Challenge (CMRxRecon) dataset.
Results show that our method achieves state-of-the-art performance, exceeding
previous methods by a significant margin. Code is available at
https://github.com/xmed-lab/DiffCMR.

###### Keywords:

Under-sampled MRI Cardiac MRI MRI reconstruction Denoising diffusion
probabilistic models

##  1 Introduction

Magnetic resonance imaging (MRI) is an important non-invasive imaging
technique that visualizes internal anatomical structures without radiation
doses. However, the acquisition time for MRI is significantly longer than
X-ray-based imaging since a series of data points should be collected in the
k-space. Particularly, cardiac magnetic resonance imaging (CMR) requires more
time for acquisition as the heart beats uncontrollably and the data
acquisition period should cover several heartbeat cycles. Long scanning time
for MRI usually brings discomfort and stress to patients, which may induce
artifacts in the MRI reconstruction process. In this work, we study the under-
sampled MRI reconstruction that sparsely samples k-space data for image
reconstruction, which is one of the ways to accelerate MRI acquisition.

The reconstruction from under-sampled k-space data is an ill-posed problem,
which is challenging and has received much attention. In the past, compressive
sensing [9] is used to solve an optimization problem that seeks to find the
most compressible representation that is consistent with the under-sampled
k-space data. With the development of deep learning, the reconstruction
problem is usually formulated as a denoising task that removes the artifacts
in under-sampled MRI images. For example, encoder-decoder-based methods [15,
13, 10] are proposed to learn a mapping from under-sampled images to fully-
sampled images. GAN-based methods [3] introduce a discriminator network that
learns to identify the differences between the generated and real images,
which can further improve image quality. However, GAN-based methods are
difficult to train and require careful tuning of hyperparameters.

Recently, denoising diffusion probabilistic models [5, 14] (DDPMs) are
proposed to use a series of transformations to increase the complexity of the
generated output iteratively. Compared with GAN-based methods, DDPMs have been
shown to be more stable during training and demonstrate outstanding
performance on a variety of computer vision tasks [4], including image
synthesis [5, 14, 2], inpainting [19, 8], segmentation [1, 17], and denoising
[18, 21, 20]. To this end, we leverage conditional DDPMs into under-sampled
MRI reconstruction. Specifically, we employ a conditional DDPM that generates
fully-sampled MRI slices with a conditioning signal from the input under-
sampled MRI slices and further adopt multi-round inference ensembling to
stabilize the denoising process.

To summarize, the main contributions of this work include 1.) we propose
DiffCMR for fast MRI reconstruction from under-sampled k-space data by
leveraging conditional diffusion models; 2.) extensive experiments are
conducted on MICCAI 2023 CMRxRecon dataset, showing DiffCMR’s state-of-the-art
performance that outperforms previous methods by a large margin.

##  2 Methodology

###  2.1 Problem Definition

In this work, we formulate the reconstruction of high-quality MRI from under-
sampled k-space data as a denoising task. Specifically, inverse Fast Fourier
Transform (iFFT) is performed to transform fully-sampled and under-sampled
k-space data into 2D image slices, which are referred to as fully-sampled and
under-sampled image slices in subsequent sections. We denote the set of
generated image slices as
$D^{(k)}=\\{(I_{i,k}^{u},I_{i,k}^{f})_{i=1}^{N_{k}}\\}$, where $k$ refers to
three different acceleration factors; $N_{k}$ is the number of samples
generated from the raw data with acceleration factor $k$; $I_{i,k}^{u}$
indicates the $i$-th under-sampled image slice with acceleration factor $k$
and $I_{i,k}^{f}$ is the corresponding fully-sampled image slice. For each
acceleration factor $k$, we aim to train a denoising model
$\phi_{\theta}^{(k)}$ with the dataset $D^{(k)}$, which has the ability to
recover fully-sampled image slice $I_{i,k}^{u}$ from under-sampled image slice
$I_{i,k}^{f}$.

###  2.2 Data Preprocess

The original dataset provided by the challenge organizer is composed of
single-coil and multi-coil data. Each data is stored in .mat format,
accompanied by striped masks of different acceleration factors (e.g., 4, 8,
and 10). Under-sampled k-space data can be generated by covering the striped
mask on the fully-sampled data. Based on the problem formulation, we first
process fully-sampled and under-sampled k-space data to make data pairs (i.e.,
ground truth and input) for network training. Specifically, for single-coil
data, we directly apply iFFT to obtain 2D image slices; for mult-coil data, we
first apply RSS [6] to aggregate k-space data from multiple coils and then
apply iFFT to obtain 2D image slices. In other words, for a multi-coil input
in shape $[t,s,c,h,w]$ (time-frame, slice, coil, height, width) or a single-
coil input in shape $[t,s,h,w]$, the preprocessed data will be $t\times s$ 2D
image slices in shape $[h,w]$.

We observe that blur noise is introduced by missing frequency information
blocked by the striped mask. Hence, for padding these slices to a fixed shape,
we choose to add zero padding to the k-space instead of the image space to
keep the purify of the source of blur noise because padding zeros in the
k-space will not bring new information from the frequency perspective while
padding zeros in the image space will introduce unnecessary bias.

Figure 1: Forward and backward denoising diffusion processes.

###  2.3 Framework — DiffCMR

We briefly introduce the formulation of the diffusion model proposed in [5,
14]. Diffusion models are generative models based on parametrized Markov chain
and are comprised of forward and backward processes, as shown in Fig. 1. The
forward process iteratively transforms a clean image $x_{0}$ into a series of
noisier images $\\{x_{1},x_{2},...,x_{T}\\}$. The following formulation can
express the iteration of the forward process:

$q(x_{t}|x_{t-1})=\mathcal{N}(x_{t};\sqrt{1-\beta_{t}}x_{t-1},\beta_{t}I_{n\times
n}),$ (1)

where $\beta_{t}$ is a constant to define the ratio of adding Gaussian noise.
The reverse process is parametrized by $\theta$ and can be simplified as:

$p_{\theta}(x_{t-1}|x_{t})=\mathcal{N}(x_{t-1};\epsilon_{\theta}(x_{t},t),\sigma_{t}^{2}I_{n\times
n}),$ (2)

where $\sigma_{t}$ is a fixed variance, and $\epsilon_{\theta}$ is a trained
step estimation function.

Conditional generations with diffusion models [2, 11, 1] are formulated by
letting backward process $p_{\theta}$ simultaneously manipulate the current
step noisy image $x_{t}$ and the given condition, which allows the generation
of images based on extra conditions without any additional learning. Here we
merge the information of the under-sampled image and the current step
denoising result by adding the extracted features from their corresponding
encoders.

Our proposed DiffCMR, as illustrated in Fig. 2, employs a conditional
diffusion model that conditions its step estimation function
$\epsilon_{\theta}$ on the aggregated information from both the input under-
sampled image slice $I^{u}$ and the current step recovery $x_{t}$. In our
architecture, the estimation function $\epsilon_{\theta}$ is a modified U-Net
[12] and can be further expressed as:

$\epsilon_{\theta}(x_{t},I^{u},t)=E(F(x_{t})+G(I^{u}),H(t)),$ (3)

where $H$ encodes the current time step $t$ into timestamp embedding; $G$ and
$F$ encode the under-sampled input slice $I^{u}$ and the current step
denoising result $x_{t}$, respectively. $E$ is a modified U-Net encoder-
decoder structure that receives summed features from $F$ and $G$, and
estimates the noise for the current step. The current time step $t$ is
embedded using a learned look-up table $H$ and inserted into layers of both
the encoder and the decoder of network $E$ by summation.

Figure 2: The figure above illustrates the pipeline of our proposed DiffCMR.
$F$ and $G$ encode the noisy signal $x_{t}$ and the under-sampled image
$I^{u}$, respectively. $H$ encodes timestamp $t$ to obtain timestamp
embeddings. $E$ is a modified U-net that receives both summed features and
timestamp embeddings for denoising.

###  2.4 Training and Inference Procedure

The training process for DiffCMR is demonstrated in Alg. 1. For each step, we
sample a random data pair $(I^{u},I^{f})$ from the dataset $D_{k}$, a
timestamp $t\in[1,T]$ from a Uniform distribution, and a noise $\epsilon$ from
a standard Gaussian distribution. We then obtain the current step recovery
$x_{t}$ by reparametrizing Eq. 1:

$\alpha_{t}=1-\beta_{t},\;\overline{\alpha}_{t}=\prod_{s=0}^{t}\alpha_{s},\;x_{t}=\sqrt{\overline{\alpha}_{t}}I^{f}+\sqrt{1-\overline{\alpha}_{t}}\epsilon.$
(4)

Then according to our pipeline in Fig. 2, $x_{t}$, $t$, and $I^{u}$ are sent
through the networks $F$, $H$, $G$, and $E$ to obtain
$\epsilon_{\theta}(x_{t},I_{i}^{u},t)$. Our training target is to minimize the
term:

$\|\epsilon-\epsilon_{\theta}(\sqrt{\overline{\alpha}_{t}}I^{f}+\sqrt{1-\overline{\alpha}_{t}}\epsilon,I_{i}^{u},t)\|.$
(5)

Our inference process is described in Alg. 2. As the procedure to recover
$x_{t-1}$ includes an addition with a random sampled
$z\sim\mathcal{N}(\mathbf{0},\mathbf{I_{n\times n}})$, which yields discrete
random noise points in the denoising results. We adopt a multi-round
ensembling inference strategy to diminish the noise points and stabilize the
denoising results. The inference procedure on the same input is carried out
for multiple rounds, and the final result is ensembled by taking the average.
The effectiveness of this strategy is proved by the experiment results from
the ablation study, see Tab.4 (b).

Algorithm 1 Training Algorithm

0: total denoising steps $T$, under-sampled and fully-sampled image pair
dataset $D_{k}=\\{(I_{i}^{u},I_{i}^{f})\\}_{i=1}^{N_{k}}$

repeat

Sample $(I_{i}^{u},I_{i}^{f})\sim
D_{k},\;\epsilon\sim\mathcal{N}(\mathbf{0},\mathbf{I_{n\times
n}}),\;t\sim\mathbf{Uniform(\\{1,...,T\\})}$

$\beta_{t}=\frac{10^{-4}(T-t)+2\times 10^{-2}(t-1)}{T-1}$,
$\alpha_{t}=1-\beta_{t}$, $\overline{\alpha}_{t}=\prod_{s=0}^{t}\alpha_{s}$

$x_{t}=\sqrt{\overline{\alpha}_{t}}I_{i}^{f}+\sqrt{1-\overline{\alpha}_{t}}\epsilon$

Compute gradient
$\nabla_{\theta}\|\epsilon-\epsilon_{\theta}(x_{t},I_{i}^{u},t)\|$

until convergence

Algorithm 2 Inference Algorithm

0: total denoising steps $T$, under-sampled image $I^{u}$, ensemble rounds $R$

for $r=1,2,...,R$ do

Sample $x_{T,r}\sim\mathcal{N}(\mathbf{0},\mathbf{I_{n\times n}})$

for $t=T,T-1,...,1$ do

Sample $z\sim N(\mathbf{0},\mathbf{I_{n\times n}})$

$\beta_{t}=\frac{10^{-4}(T-t)+2\times 10^{-2}(t-1)}{T-1}$

$\alpha_{t}=1-\beta_{t}$, $\overline{\alpha}_{t}=\prod_{s=0}^{t}\alpha_{s}$,
$\widetilde{\beta_{t}}=\frac{1-\overline{\alpha}_{t-1}}{1-\overline{\alpha}_{t}}\beta_{t}$

$x_{t-1,r}=\alpha_{t}^{\frac{1}{2}}(x_{t}-\frac{1-\alpha_{t}}{\sqrt{1-\overline{\alpha}_{t}}}\epsilon_{\theta}(x_{t},I^{u},t))+\mathbbm{1}_{[t>1]}\widetilde{\beta}_{t}^{\frac{1}{2}}z$

return $\sum_{r=1}^{R}x_{0,r}/R$

##  3 Experiments

###  3.1 Dataset

The CMRxRecon dataset is released in the MICCAI 2023 CMRxRecon Challenge,
comprises cine reconstruction and T1/T2 mapping tasks. Both tasks have two
coil types, each with 120 cases for training and 60 cases for validation. The
training cases provide fully sampled k-space data and under-sampled k-space
data with acceleration factors (AccFactor) 4, 8, and 10. The validation cases
only provide under-sampled k-space data with acceleration factors 4, 8, and
10.

As described in Sec. 2.2, we first preprocess the raw data by zero-padding the
k-space to size 512$\times$512, and transforming it to 2D images with iFFT.
Finally, we resize the image slices to 128$\times$128 to speed up the
experiment process. Our local training-validation split is set by assigning
images extracted from case P001 to P110 to the training set and those from
case P111 to P120 to the validation set. We randomly shuffle the validation
set and select the first 240 samples for inference. The detailed numbers of
our preprocessed samples for training and validation are listed in Table 1.

Table 1: Training and Validation pairs for both tasks and acceleration factors (AccFactor) 4, 8 and 10. | AccFactor04 | AccFactor08 | AccFactor10
---|---|---|---
| Train | Valid | Train | Valid | Train | Valid
Task 1 | 14304 | 1272 | 14304 | 1272 | 14304 | 1272
Task 2 | 32904 | 2904 | 32904 | 2904 | 32904 | 2904

###  3.2 Implementation Details

As for model architectures, we follow the conventions in [1] to build our
model. The network $G$ has 10 Residual in Residual Dense Blocks [16] and a
depth of six. The number of channels was set to $[C,C,2C,2C,4C,4C]$ with
$C=128$. The augmentation schemes include horizontal and vertical flips with a
probability of 0.5. The training process took place with a batch size of 6
images at size 128 $\times$ 128\. We used AdamW [7] optimizer for all
experiments. We used 100 diffusion steps for training and 1000 for inference.
As there are two different tasks together with 3 acceleration factors, we
trained the network 6 times to obtain 6 different sets of weight. The whole
training and inference process is carried out on a single NVIDIA GeForce RTX
3090 GPU.

Table 2: Experiment results for Task 1 - Cine Reconstruction | AccFactor04 | AccFactor08 | AccFactor10
---|---|---|---
| PSNR$\uparrow$ | SSIM$\uparrow$ | NMSE$\downarrow$ | PSNR$\uparrow$ | SSIM$\uparrow$ | NMSE$\downarrow$ | PSNR$\uparrow$ | SSIM$\uparrow$ | NMSE$\downarrow$
RAW | 28.79 | 0.8150 | 0.2187 | 27.94 | 0.8082 | 0.2781 | 27.75 | 0.8113 | 0.2883
U-Net [12] | 33.11 | 0.9212 | 0.0673 | 33.31 | 0.9460 | 0.0646 | 32.71 | 0.9391 | 0.0740
cGAN [3] | 33.85 | 0.9435 | 0.0573 | 33.14 | 0.9449 | 0.0669 | 32.56 | 0.9209 | 0.0760
DiffCMR | 36.10 | 0.9277 | 0.0346 | 34.85 | 0.9061 | 0.0457 | 34.47 | 0.9016 | 0.0493
Table 3: Experiment results for Task 2 - T1/T2 Mapping | AccFactor04 | AccFactor08 | AccFactor10
---|---|---|---
| PSNR$\uparrow$ | SSIM$\uparrow$ | NMSE$\downarrow$ | PSNR$\uparrow$ | SSIM$\uparrow$ | NMSE$\downarrow$ | PSNR$\uparrow$ | SSIM$\uparrow$ | NMSE$\downarrow$
RAW | 28.17 | 0.8167 | 0.2003 | 27.15 | 0.8041 | 0.2578 | 27.17 | 0.8100 | 0.2711
U-Net [12] | 32.06 | 0.9340 | 0.0537 | 31.05 | 0.9286 | 0.0695 | 29.49 | 0.9106 | 0.0987
cGAN [3] | 32.67 | 0.9460 | 0.0488 | 30.65 | 0.8899 | 0.0805 | 31.38 | 0.9304 | 0.0655
DiffCMR | 34.60 | 0.9071 | 0.0372 | 33.17 | 0.8937 | 0.0537 | 33.04 | 0.8941 | 0.0536

###  3.3 Results

As the online validation platform limits the daily attempts to 3 trials per
task, we perform the validation and report the results on our local split
dataset to speed up the upgrading procedure of our method. For the fairness of
comparison, all the experiments are trained and validated with the same split,
input resolution, and the same data augmentation scheme. The evaluation
metrics for our experiments are peak signal-to-noise ratio (PSNR), structural
similarity index measure (SSIM), and normalized mean square error (NMSE).

We compare our proposed DiffCMR with U-Net [12] and cGAN [3]. Qualitative
visualization comparisons are shown in Fig. 3. Quantitative results for Task1
and Task2 are listed in Tab. 2 and Tab. 3, where RAW means the direct
comparison results between the input under-sampled image slices and the fully-
sampled image slices. As can be seen, our method outperforms both baseline
methods across most tasks and acceleration factors.

Figure 3: Qualitative comparisons between raw input, results of U-Net [12],
results of cGAN [3], results of DiffCMR with single-round ensembling inference
($R=1$), results of DiffCMR with multi-round ensembling inference ($R=4$), and
ground truth.

##  4 Ablation Study

We evaluate two alternatives of hyper-parameters in our DiffCMR method at the
inference stage. The first variant determines the number of inference
diffusion steps. The second variant determines the number of inference
ensembling rounds. These variant experiments were carried out on the data from
Task 1 with an acceleration factor equal to 4.

Table 4: Ablation study on Task 1-AccFactor04. (a) Results with different
inference diffusion steps. (b) Results with different inference ensembling
rounds.

#Steps $T$ | PSNR$\uparrow$ | SSIM$\uparrow$ | NMSE$\downarrow$
---|---|---|---
Raw | 28.79 | 0.8150 | 0.2187
20 | 25.52 | 0.3664 | 0.3830
100 | 29.86 | 0.6006 | 0.1393
500 | 35.39 | 0.9010 | 0.0399
1000 | 36.10 | 0.9277 | 0.0346
(a)

#Rounds $R$ | PSNR$\uparrow$ | SSIM$\uparrow$ | NMSE$\downarrow$
---|---|---|---
Raw | 28.79 | 0.8150 | 0.2187
1 | 33.89 | 0.8491 | 0.0561
2 | 35.25 | 0.8994 | 0.0417
4 | 36.10 | 0.9277 | 0.0346
8 | 36.68 | 0.9430 | 0.0305
(b)

Varying the number of inference diffusion steps $T$. In this part, we set the
ensembling rounds $R=4$ and explore the effect of inference diffusion steps on
the denoising quality. Quantitative results are shown in Tab. 4(a). As can be
observed, the denoising performance is positively correlated with $T$ and
starts to outperform the raw input at around $T=100$. We choose $T=1000$ in
all other experiments.

Varying the number of inference ensembling rounds $R$. In this part, we set
the diffusion steps $T=1000$ and explore the effect of inference ensembling
rounds on the denoising quality. Quantitative results are shown in Tab. 4(b),
and Qualitative comparisons between $R=1$ and $R=4$ are visualized in Fig. 3.
The results show a positive correlation between denoising effectiveness and
$R$ with a diminishing marginal effect when $R$ is large. Therefore, we set
$R=4$ in all other experiments for a reasonable performance to runtime
tradeoff.

##  5 Conclusion

In this paper, we present DiffCMR, a conditional DDPM-based approach for high-
quality MRI reconstruction from under-sampled k-space data. Our framework
receives conditioning signals from the under-sampled MRI image slice at each
denoising diffusion step and generates the corresponding fully-sampled MRI
image slice. In addition, we adopt the multi-round ensembling strategy during
inference which largely enhances the stableness of our approach. Experiment
results show our DiffCMR method outperforms the existing popular denoisers
qualitatively and quantitatively. In conclusion, our proposed DiffCMR offers a
novel perspective for handling fast MRI reconstruction problems and
demonstrates impressive robustness.

Acknowledgement. This work was supported by the Hong Kong Innovation and
Technology Fund under Projects PRP/041/22FX.

## References

  * [1] Amit, T., Shaharbany, T., Nachmani, E., Wolf, L.: Segdiff: Image segmentation with diffusion probabilistic models. arXiv preprint arXiv:2112.00390 (2021)
  * [2] Choi, J., Kim, S., Jeong, Y., Gwon, Y., Yoon, S.: Ilvr: Conditioning method for denoising diffusion probabilistic models. arXiv preprint arXiv:2108.02938 (2021)
  * [3] Defazio, A., Murrell, T., Recht, M.: Mri banding removal via adversarial training. In: Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., Lin, H. (eds.) Advances in Neural Information Processing Systems. vol. 33, pp. 7660–7670. Curran Associates, Inc. (2020)
  * [4] Dhariwal, P., Nichol, A.: Diffusion models beat gans on image synthesis. Advances in neural information processing systems 34, 8780–8794 (2021)
  * [5] Ho, J., Jain, A., Abbeel, P.: Denoising diffusion probabilistic models. Advances in neural information processing systems 33, 6840–6851 (2020)
  * [6] Larsson, E.G., Erdogmus, D., Yan, R., Principe, J.C., Fitzsimmons, J.R.: Snr-optimality of sum-of-squares reconstruction for phased-array magnetic resonance imaging. Journal of Magnetic Resonance 163(1), 121–123 (2003)
  * [7] Loshchilov, I., Hutter, F.: Decoupled weight decay regularization. arXiv preprint arXiv:1711.05101 (2017)
  * [8] Lugmayr, A., Danelljan, M., Romero, A., Yu, F., Timofte, R., Van Gool, L.: Repaint: Inpainting using denoising diffusion probabilistic models. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 11461–11471 (2022)
  * [9] Lustig, M., Donoho, D., Pauly, J.M.: Sparse mri: The application of compressed sensing for rapid mr imaging. Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 58(6), 1182–1195 (2007)
  * [10] Muckley, M.J., Ades-Aron, B., Papaioannou, A., Lemberskiy, G., Solomon, E., Lui, Y.W., Sodickson, D.K., Fieremans, E., Novikov, D.S., Knoll, F.: Training a neural network for gibbs and noise removal in diffusion mri. Magnetic resonance in medicine 85(1), 413–428 (2021)
  * [11] Nichol, A.Q., Dhariwal, P.: Improved denoising diffusion probabilistic models. In: International Conference on Machine Learning. pp. 8162–8171. PMLR (2021)
  * [12] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18. pp. 234–241. Springer (2015)
  * [13] Schlemper, J., Caballero, J., Hajnal, J.V., Price, A.N., Rueckert, D.: A deep cascade of convolutional neural networks for dynamic mr image reconstruction. IEEE transactions on Medical Imaging 37(2), 491–503 (2017)
  * [14] Song, J., Meng, C., Ermon, S.: Denoising diffusion implicit models. arXiv preprint arXiv:2010.02502 (2020)
  * [15] Sriram, A., Zbontar, J., Murrell, T., Defazio, A., Zitnick, C.L., Yakubova, N., Knoll, F., Johnson, P.: End-to-end variational networks for accelerated mri reconstruction. In: Martel, A.L., Abolmaesumi, P., Stoyanov, D., Mateus, D., Zuluaga, M.A., Zhou, S.K., Racoceanu, D., Joskowicz, L. (eds.) Medical Image Computing and Computer Assisted Intervention – MICCAI 2020. pp. 64–73. Springer International Publishing, Cham (2020)
  * [16] Wang, X., Yu, K., Wu, S., Gu, J., Liu, Y., Dong, C., Qiao, Y., Change Loy, C.: Esrgan: Enhanced super-resolution generative adversarial networks. In: Proceedings of the European conference on computer vision (ECCV) workshops. pp. 0–0 (2018)
  * [17] Wu, J., Fang, H., Zhang, Y., Yang, Y., Xu, Y.: Medsegdiff: Medical image segmentation with diffusion probabilistic model. arXiv preprint arXiv:2211.00611 (2022)
  * [18] Xiang, T., Yurt, M., Syed, A.B., Setsompop, K., Chaudhari, A.: DDM$^2$: Self-supervised diffusion MRI denoising with generative diffusion models (2023)
  * [19] Xie, S., Zhang, Z., Lin, Z., Hinz, T., Zhang, K.: Smartbrush: Text and shape guided object inpainting with diffusion model. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 22428–22437 (2023)
  * [20] Xie, Y., Yuan, M., Dong, B., Li, Q.: Diffusion model for generative image denoising. arXiv preprint arXiv:2302.02398 (2023)
  * [21] Zhu, Y., Zhang, K., Liang, J., Cao, J., Wen, B., Timofte, R., Van Gool, L.: Denoising diffusion models for plug-and-play image restoration. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 1219–1229 (2023)


#  On Isolated Singularities of Fractional Semi-Linear Elliptic
Equations††thanks: Supported by NSFC. E-mail addresses: hui-
<EMAIL_ADDRESS>(H. Yang<EMAIL_ADDRESS>(W. Zou)

Hui Yang1, Wenming Zou2
1Yau Mathematical Sciences Center, Tsinghua University, Beijing 100084, China
2Department of Mathematical Sciences, Tsinghua University, Beijing 100084,
China

Abstract

In this paper, we study the local behavior of nonnegative solutions of
fractional semi-linear equations $(-\Delta)^{\sigma}u=u^{p}$ with an isolated
singularity, where $\sigma\in(0,1)$ and
$\frac{n}{n-2\sigma}<p<\frac{n+2\sigma}{n-2\sigma}$. We first use blow up
method and a Liouville type theorem to derive an upper bound. Then we
establish a monotonicity formula and a sufficient condition for removable
singularity to give a classification of the isolated singularities. When
$\sigma=1$, this classification result has been proved by Gidas and Spruck
(Comm. Pure Appl. Math. 34: 525-598, 1981).

Mathematics Subject Classification (2010): 35B09; 35B40; 35J70; 35R11

##  1 Introduction and Main results

The purpose of this paper is to study the local behavior of nonnegative
solutions of

$(-\Delta)^{\sigma}u=u^{p}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textmd{in}~{}B_{1}\backslash\\{0\\}$
(1.1)

with an isolated singularity at the origin, where the punctured unit ball
$B_{1}\backslash\\{0\\}\subset\mathbb{R}^{n}$ with $n\geq 2$, $\sigma\in(0,1)$
and $(-\Delta)^{\sigma}$ is the fractional Laplacian.

When $\sigma=1$, the isolated singularity of nonnegative solutions for (1.1)
has been very well understand, see Lions [26] for $1<p<\frac{n}{n-2}$, Aviles
[4] for $p=\frac{n}{n-2}$, Gidas-Spruck [19] for
$\frac{n}{n-2}<p<\frac{n+2}{n-2}$, Caffarelli-Gidas-Spruck [7] for
$\frac{n}{n-2}\leq p\leq\frac{n+2}{n-2}$, Korevaar-Mazzeo-Pacard-Schoen [25]
for $p=\frac{n+2}{n-2}$, and Bidaut-Véron and Véron [5] for
$p>\frac{n+2}{n-2}$.

The semi-linear equation (1.1) involving the fractional Laplacian has
attracted a great deal of interest since they are of central importance in
many fields, such as see [1, 2, 3, 8, 10, 12, 13, 14, 15, 17, 20, 21, 22, 23,
24] and references therein. Recently, the existence of singular solutions of
equation (1.1) with prescribed isolated singularities for the critical
exponent $p=\frac{n+2\sigma}{n-2\sigma}$ were studied in [1, 2, 14, 15] and
for the subcritical regime $\frac{n}{n-2\sigma}<p<\frac{n+2\sigma}{n-2\sigma}$
were studied in [1, 3]. Solutions of (1.1) with an isolated singularity are
the simplest cases of those singular solutions. In a recent paper [8],
Caffarelli, Jin, Sire and Xiong study the local behavior of nonnegative
solution of (1.1) with $p=\frac{n+2\sigma}{n-2\sigma}$. More precisely, let
$u$ be a nonnegative solution of (1.1) with $p=\frac{n+2\sigma}{n-2\sigma}$
and suppose that the origin is not a removable singularity. Then, near the
origin

$c_{1}|x|^{-\frac{n-2\sigma}{2}}\leq u(x)\leq
c_{2}|x|^{-\frac{n-2\sigma}{2}},$ (1.2)

where $c_{1},c_{2}$ are positive constants.

In this paper, we are interested in the local behavior of nonnegative
solutions of (1.1) with $\frac{n}{n-2\sigma}<p<\frac{n+2\sigma}{n-2\sigma}$.
For the classical case $\sigma=1$, this has been proved in the pioneering
paper [19] by Gidas and Spruck.

We study the equation (1.1) via the well known extension theorem for the
fractional Laplacian $(-\Delta)^{\sigma}$ established by Caffarelli-Silvestre
[9]. We use capital letters, such as
$X=(x,t)\in\mathbb{R}^{n}\times\mathbb{R}_{+}$, to denote points in
$\mathbb{R}_{+}^{n+1}$. We also denote $\mathcal{B}_{R}$ as the ball in
$\mathbb{R}^{n+1}$ with radius $R$ and center at the origin,
$\mathcal{B}_{R}^{+}$ as the upper half-ball
$\mathcal{B}_{R}\cap\mathbb{R}_{+}^{n+1}$, and
$\partial^{0}\mathcal{B}_{R}^{+}$ as the flat part of
$\partial\mathcal{B}_{R}^{+}$ which is the ball $B_{R}$ in $\mathbb{R}^{n}$.
Then the problem (1.1) is equivalent to the following extension problem

$\begin{cases}-\textmd{div}(t^{1-2\sigma}\nabla
U)=0~{}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{in}~{}\mathcal{B}_{1}^{+},\\\
\frac{\partial
U}{\partial\nu^{\sigma}}(x,0)=U^{p}(x,0)~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{on}~{}\partial^{0}\mathcal{B}_{1}^{+}\backslash\\{0\\},\end{cases}$
(1.3)

where $\frac{\partial U}{\partial\nu^{\sigma}}(x,0):=-\lim_{t\rightarrow
0^{+}}t^{1-2\sigma}\partial_{t}U(x,t)$. By [9], we only need to analyze the
behavior of the traces

$u(x):=U(x,0)$

of the nonnegative solutions $U(x,t)$ of (1.3) near the origin, from which we
can get the behavior of solutions of (1.1) near the origin.

We say that $U$ is a nonnegative solution of (1.3) if $U$ is in the weighted
Sobolev space
$W^{1,2}(t^{1-2\sigma},\mathcal{B}^{+}_{1}\backslash\overline{\mathcal{B}}^{+}_{\epsilon})$
for every $\epsilon>0$, $U\geq 0$, and it satisies (1.3) in the sense of
distribution away from 0, i.e., for every nonnegative $\Phi\in
C_{c}^{\infty}\left((\mathcal{B}_{1}^{+}\cup\partial^{0}\mathcal{B}_{1}^{+})\backslash\\{0\\}\right)$,

$\int_{\mathcal{B}_{1}^{+}}t^{1-2\sigma}\nabla
U\nabla\Phi=\int_{\partial^{0}\mathcal{B}_{1}^{+}}U^{p}\Phi.$ (1.4)

See [23] for more details on this definition. Then it follows from the
regularity result in [23] that $U$ is locally Hölder continuous in
$\overline{\mathcal{B}}^{+}_{1}\backslash\\{0\\}$. We say that the origin 0 is
a removable singularity of solution $U$ of (1.3) if $U(x,0)$ can be extended
as a continuous function near the origin, otherwise we say that the origin 0
is a non-removable singularity. Our main result is the following

######  Theorem 1.1.

Let $U$ be a nonnegative solution of (1.3). Assume

$\frac{n}{n-2\sigma}<p<\frac{n+2\sigma}{n-2\sigma}.$

Then either the singularity near 0 is removable, or there exist two positive
constants $c_{1}$ and $c_{2}$ such that

$c_{1}|x|^{-\frac{2\sigma}{p-1}}\leq u(x)\leq
c_{2}|x|^{-\frac{2\sigma}{p-1}}.$ (1.5)

######  Remark 1.1.

We point out that, if (1.5) holds, then the Harnack inequality (3.2) implies
that

$C_{1}|X|^{-\frac{2\sigma}{p-1}}\leq U(X)\leq C_{2}|X|^{-\frac{2\sigma}{p-1}}$

holds as well, for some positive constants $C_{1}$ and $C_{2}$.

For the classical case $\sigma=1$, Theorem 1.1 were proved in [19] by Gidas
and Spruck. We may also see [7] for this classical case. The similar upper
bound in (1.5) obtained in [19] for the classical case is very complicated and
technical, here we use the blow up method and a Liouville type theorem to
prove the upper bound in (1.5). To obtain the lower bound, there are some
extra difficulties, one of which is that the Pohozaev identity is not
available. More precisely, for the critical case
$p=\frac{n+2\sigma}{n-2\sigma}$, the Pohozaev integral $P(U,R)$ is independent
of $R$ by the Pohozaev identity (see [8] for more detalis on $P(U,R)$). In
[8], the authors make use of this property of the Pohozaev integral to prove
the lower bound, however, this does not hold in the subcritical case. We will
establish a useful monotonicity formula to overcome this difficulty. The
others would be those extra techniques to get the estimates of $U$ from those
of its trace $u$.

The paper is organized as follow. In Section 2, we recall three propositions:
a Liouville theorem, a Harnack inequality and a Sobolev inequality. Section 3
is devoted to the proof of Theorem 1.1. We first derive an upper bound and a
special form of Harnack inequality. Then we establish a monotonicity formula
and a sufficient condition of removable singularity to prove Theorem 1.1.

##  2 Preliminaries

In this section, we introduce some notations and some propositions which will
be used in our arguments. We denote $\mathcal{B}_{R}$ as the ball in
$\mathbb{R}^{n+1}$ with radius $R$ and center 0, and $B_{R}$ as the ball in
$\mathbb{R}^{n}$ with radius $R$ and center 0. We also denote
$\mathcal{B}_{R}^{+}$ as the upper half-ball
$\mathcal{B}_{R}\cap\mathbb{R}_{+}^{n+1}$,
$\partial^{+}\mathcal{B}_{R}^{+}=\partial\mathcal{B}_{R}^{+}\cap\mathbb{R}_{+}^{n+1}$
as the positive part of $\partial\mathcal{B}_{R}^{+}$, and
$\partial^{0}\mathcal{B}_{R}^{+}=\partial\mathcal{B}_{R}^{+}\backslash\partial^{+}\mathcal{B}_{R}^{+}$
as the flat part of $\partial\mathcal{B}_{R}^{+}$ which is the ball $B_{R}$ in
$\mathbb{R}^{n}$.

We say $U\in W_{loc}^{1,2}(t^{1-2\sigma},\overline{\mathbb{R}^{n+1}_{+}})$ if
$U\in W^{1,2}(t^{1-2\sigma},\mathcal{B}_{R}^{+})$ for all $R>0$, and $U\in
W_{loc}^{1,2}(t^{1-2\sigma},\overline{\mathbb{R}^{n+1}_{+}}\backslash\\{0\\})$
if $U\in
W^{1,2}(t^{1-2\sigma},\mathcal{B}_{R}^{+}\backslash\overline{\mathcal{B}}_{\epsilon}^{+})$
for all $R>\epsilon>0$.

We next recall three propositions, which will be used frequently in our paper.
For convenience, we state them here. Their proofs can be found in [23]. The
first one is a Liouville type theorem.

######  Proposition 2.1.

Let $U\in W_{loc}^{1,2}(t^{1-2\sigma},\overline{\mathbb{R}^{n+1}_{+}})$ be a
nonnegative weak solution of

$\begin{cases}-\textmd{div}(t^{1-2\sigma}\nabla
U)=0~{}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{in}~{}\mathbb{R}^{n+1}_{+},\\\
\frac{\partial
U}{\partial\nu^{\sigma}}(x,0)=U^{p}(x,0)~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{on}~{}\mathbb{R}^{n},\end{cases}$
(2.1)

with

$1\leq p<\frac{n+2\sigma}{n-2\sigma}.$

Then

$U(x,t)\equiv 0.$

The second one is a Harnack inequality, see also [6].

######  Proposition 2.2.

Let $U\in W_{loc}^{1,2}(t^{1-2\sigma},\mathcal{B}^{+}_{1})$ be a nonnegative
weak solution of

$\begin{cases}-\textmd{div}(t^{1-2\sigma}\nabla
U)=0~{}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{in}~{}\mathcal{B}^{+}_{1},\\\
\frac{\partial
U}{\partial\nu^{\sigma}}(x,0)=a(x)U(x,0)~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{on}~{}\partial^{0}\mathcal{B}^{+}_{1},\end{cases}$
(2.2)

If $a\in L^{q}(B_{1})$ for some $q>\frac{n}{2\sigma}$, then we have

$\sup_{\mathcal{B}^{+}_{1/2}}U\leq C\inf_{\mathcal{B}^{+}_{1/2}}U,$ (2.3)

where $C$ depends only on $n,\sigma$ and $\|a\|_{L^{q}(B_{1})}$.

The last one is a Sobolev type inequality.

######  Proposition 2.3.

Let $D=\Omega\times(0,R)\subset\mathbb{R}^{n}\times\mathbb{R}_{+}$ with $R>0$
and $\partial\Omega$ Lipschitz. Then there exists $C_{n,\sigma}>0$ depending
only on $n$ and $\sigma$ such that

$\|U(\cdot,0)\|_{L^{2n/(n-2\sigma)}(\Omega)}\leq C_{n,\sigma}\|\nabla
U\|_{L^{2}(t^{1-2\sigma},D)}$

for all $U\in C_{c}^{\infty}(D\cup\partial^{0}D)$.

##  3 Classification of Isolated Singularities

In this section, we investigate the local behavior of nonnegative solutions of
(1.3) and classify their isolated singularities. We first prove an upper bound
and a special form of Harnack inequality for nonnegative solutions with a
possible isolated singularity. We remark that this result will be of basic
importance in classifying the isolated singularities.

######  Proposition 3.1.

Let $U$ be a nonnegative solution of (1.3), with
$1<p<\frac{n+2\sigma}{n-2\sigma}$. Then,

  * (1)

there exists a positive constant $c$ independent of $U$ such that

$u(x)\leq
c|x|^{-\frac{2\sigma}{p-1}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textmd{in}~{}B_{1/2};$
(3.1)

  * (2)

(Harnack inequality) for all $0<r<1/8$, we have

$\sup_{\mathcal{B}_{2r}^{+}\backslash\overline{\mathcal{B}}_{r/2}^{+}}U\leq
C\inf_{\mathcal{B}_{2r}^{+}\backslash\overline{\mathcal{B}}_{r/2}^{+}}U,$
(3.2)

where $C$ is a positive constant independent of $r$ and $U$. In particular,
for all $0<r<1/8$, we have

$\sup_{\partial^{+}\mathcal{B}_{r}^{+}}U\leq
C\inf_{\partial^{+}\mathcal{B}_{r}^{+}}U,$ (3.3)

where $C$ is a positive constant independent of $r$ and $U$.

###### Proof.

Suppose by contradiction that there exists a sequence of points
$\\{x_{k}\\}\subset B_{1/2}$ and a sequence of solutions $\\{U_{k}\\}$ of
(1.3), such that

$|x_{k}|^{\frac{2\sigma}{p-1}}u_{k}(x_{k})\rightarrow+\infty~{}~{}~{}~{}~{}~{}\textmd{as}~{}k\rightarrow\infty.$
(3.4)

As in [8], we define

$v_{k}(x):=\left(\frac{|x_{k}|}{2}-|x-x_{k}|\right)^{\frac{2\sigma}{p-1}}u_{k}(x),~{}~{}~{}~{}~{}|x-x_{k}|\leq\frac{|x_{k}|}{2}.$

Take $\bar{x}_{k}$ satisfy $|\bar{x}_{k}-x_{k}|<\frac{|x_{k}|}{2}$ and

$v_{k}(\bar{x}_{k})=\max_{|x-x_{k}|\leq\frac{|x_{k}|}{2}}v_{k}(x).$

Let

$2\mu_{k}:=\frac{|x_{k}|}{2}-|\bar{x}_{k}-x_{k}|.$

Then

$0<2\mu_{k}\leq\frac{|x_{k}|}{2}~{}~{}~{}\textmd{and}~{}~{}~{}\frac{|x_{k}|}{2}-|x-x_{k}|\geq\mu_{k}~{}~{}~{}~{}~{}\forall~{}~{}|x-\bar{x}_{k}|\leq\mu_{k}.$

By the definition of $v_{k}$, we have

$(2\mu_{k})^{\frac{2\sigma}{p-1}}u_{k}(\bar{x}_{k})=v_{k}(\bar{x}_{k})\geq
v_{k}(x)\geq(\mu_{k})^{\frac{2\sigma}{p-1}}u_{k}(x)~{}~{}~{}~{}~{}\forall~{}~{}|x-\bar{x}_{k}|\leq\mu_{k}.$

Hence, we obtain

$2^{\frac{2\sigma}{p-1}}u_{k}(\bar{x}_{k})\geq
u_{k}(x)~{}~{}~{}~{}~{}\forall~{}~{}|x-\bar{x}_{k}|\leq\mu_{k}.$ (3.5)

Moreover, by (3.4), we also have

$(2\mu_{k})^{\frac{2\sigma}{p-1}}u_{k}(\bar{x}_{k})=v_{k}(\bar{x}_{k})\geq
v_{k}(x_{k})=\left(\frac{|x_{k}|}{2}\right)^{\frac{2\sigma}{p-1}}u_{k}(x_{k})\rightarrow+\infty~{}~{}~{}~{}\textmd{as}~{}k\rightarrow\infty.$
(3.6)

Now, we define

$W_{k}(y,t):=\frac{1}{u_{k}(\bar{x}_{k})}U_{k}\left(\bar{x}_{k}+\frac{y}{u_{k}(\bar{x}_{k})^{\frac{p-1}{2\sigma}}},\frac{t}{u_{k}(\bar{x}_{k})^{\frac{p-1}{2\sigma}}}\right),~{}~{}~{}~{}~{}(y,t)\in\Omega_{k},$

where

$\Omega_{k}:=\left\\{(y,t)\in\mathbb{R}^{n+1}_{+}|\left(\bar{x}_{k}+\frac{y}{u_{k}(\bar{x}_{k})^{\frac{p-1}{2\sigma}}},\frac{t}{u_{k}(\bar{x}_{k})^{\frac{p-1}{2\sigma}}}\right)\in\mathcal{B}_{1}^{+}\backslash\\{0\\}\right\\}.$

Let $w_{k}(y):=W_{k}(y,0)$. Then $W_{k}$ satisfies $w_{k}(0)=1$ and

$\begin{cases}-\textmd{div}(t^{1-2\sigma}\nabla
W_{k})=0~{}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{in}~{}\Omega_{k},\\\
\frac{\partial
W_{k}}{\partial\nu^{\sigma}}(x,0)=w_{k}^{p}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{on}~{}\partial^{0}\Omega_{k}.\end{cases}$
(3.7)

Furthermore, by (3.5) and (3.6), we have

$w_{k}(y)\leq
2^{\frac{2\sigma}{p-1}}~{}~{}~{}~{}~{}~{}~{}\textmd{in}~{}B_{R_{k}}$

with

$R_{k}:=\mu_{k}u(\bar{x}_{k})^{\frac{p-1}{2\sigma}}\rightarrow+\infty~{}~{}~{}~{}\textmd{as}~{}k\rightarrow\infty.$

By Proposition 2.2, for any $T>0$, we have

$0\leq W_{k}\leq C(T)~{}~{}~{}~{}~{}\textmd{in}~{}B_{R_{k}/2}\times[0,T),$

where the constant $C(T)$ depends only on $n,\sigma$ and $T$. By Corollary
2.10 and Theorem 2.15 in [23] there exists $\alpha>0$ such that for every
$R>1$,

$\|W_{k}\|_{W^{1,2}(t^{1-2\sigma},\mathcal{B}_{R}^{+})}+\|W_{k}\|_{C^{\alpha}(\overline{\mathcal{B}}_{R}^{+})}+\|w_{k}\|_{C^{2,\alpha}(\overline{B}_{R})}\leq
C(R),$

where $C(R)$ is independent of $k$. Therefore, there is a subsequence of
$k\rightarrow\infty$, still denoted by itself, and a nonnegative function
$W\in W_{loc}^{1,2}(t^{1-2\sigma},\overline{\mathbb{R}^{n+1}_{+}})\cap
C_{loc}^{\alpha}(\overline{\mathbb{R}^{n+1}_{+}})$ such that, as
$k\rightarrow\infty$,

$\begin{cases}W_{k}\rightarrow
W~{}~{}~{}~{}~{}~{}~{}\textmd{weakly}~{}\textmd{in}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}&W_{loc}^{1,2}(t^{1-2\sigma},\overline{\mathbb{R}^{n+1}_{+}}),\\\
W_{k}\rightarrow
W~{}~{}~{}~{}~{}~{}~{}\textmd{in}~{}~{}~{}~{}~{}~{}~{}&C_{loc}^{\alpha/2}(\overline{\mathbb{R}^{n+1}_{+}}),\\\
w_{k}\rightarrow
w~{}~{}~{}~{}~{}~{}~{}~{}~{}\textmd{in}~{}~{}~{}~{}~{}~{}~{}&C_{loc}^{2}(\mathbb{R}^{n}),\end{cases}$

where $w(y)=W(y,0)$. Moreover, $W$ satisfies $w(0)=1$ and

$\begin{cases}-\textmd{div}(t^{1-2\sigma}\nabla
W)=0~{}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{in}~{}\mathbb{R}^{n+1}_{+},\\\
\frac{\partial
W}{\partial\nu^{\sigma}}(x,0)=w^{p}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{on}~{}\mathbb{R}^{n}.\end{cases}$
(3.8)

This contradicts Proposition 2.1 and proves part (1) of the proposition.

Now we prove the Harnack inequality, which is actually a consequence of the
upper bound (3.1). Let

$V_{r}(X)=U(rX)$

for each $r\in(0,\frac{1}{8}]$ and for $\frac{1}{4}\leq|X|\leq 4$. Obviously,
$V_{r}$ satisfies

$\begin{cases}-\textmd{div}(t^{1-2\sigma}\nabla
V_{r})=0~{}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{in}~{}\mathcal{B}_{4}\backslash\overline{\mathcal{B}}_{1/4},\\\
\frac{\partial
V_{r}}{\partial\nu^{\sigma}}(x,0)=a_{r}(x)v_{r}(x)~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{on}~{}B_{4}\backslash\overline{B}_{1/4},\\\
\end{cases}$ (3.9)

where $v_{r}(x)=V_{r}(x,0)$ and
$a_{r}(x)=r^{2\sigma}\left(u(rx)\right)^{p-1}$. It follows (3.1) that

$|a_{r}(x)|\leq
C~{}~{}~{}~{}~{}~{}~{}\textmd{for}~{}\text{all}~{}1/4\leq|x|\leq 4,$

where $C$ is a positive constant independent of $r$ and $U$. By the Harnack
inequality in Proposition 2.2 and the standard Harnack inequality for
uniformly elliptic equations, we have

$\sup_{\frac{1}{2}\leq|X|\leq 2}V_{r}(X)\leq C\inf_{\frac{1}{2}\leq|X|\leq
2}V_{r}(X),$

where $C$ is another positive constant independent of $r$ and $U$. Hence, we
get (3.2). ∎

In order to prove the lower bound in (1.5), we need to establish a
monotonicity formula for the nonnegative solutions $U$ of (1.3). More
precisely, take a nonnegative solutions $U$ of (1.3), let $0<r<1$ and define

$\displaystyle E(r;U):=$ $\displaystyle
r^{2\frac{(p+1)\sigma}{p-1}-n}\left[r\int_{\partial^{+}\mathcal{B}_{r}^{+}}t^{1-2\sigma}|\frac{\partial
U}{\partial\nu}|^{2}+\frac{2\sigma}{p-1}\int_{\partial^{+}\mathcal{B}_{r}^{+}}t^{1-2\sigma}\frac{\partial
U}{\partial\nu}U\right]$
$\displaystyle+\frac{1}{2}\frac{2\sigma}{p-1}\left(\frac{4\sigma}{p-1}-(n-2\sigma)\right)r^{2\frac{(p+1)\sigma}{p-1}-n-1}\int_{\partial^{+}\mathcal{B}_{r}^{+}}t^{1-2\sigma}U^{2}$
$\displaystyle-r^{2\frac{(p+1)\sigma}{p-1}-n+1}\left[\frac{1}{2}\int_{\partial^{+}\mathcal{B}_{r}^{+}}t^{1-2\sigma}|\nabla
U|^{2}-\frac{1}{p+1}\int_{\partial B_{r}}u^{p+1}\right].$

Then, we have the following monotonicity formula.

######  Proposition 3.2.

Let $U$ be a nonnegative solution of (1.3) with
$1<p<\frac{n+2\sigma}{n-2\sigma}$. Then, $E(r;U)$ is non-decreasing in
$r\in(0,1)$.

###### Proof.

Take standard polar coordinates in $\mathbb{R}^{n+1}_{+}$: $X=(x,t)=r\theta$,
where $r=|X|$ and $\theta=\frac{X}{|X|}$. Let $\theta_{1}=\frac{t}{|X|}$
denote the component of $\theta$ in the $t$ direction and

$\mathbb{S}^{n}_{+}=\\{X\in\mathbb{R}^{n+1}_{+}:r=1,\theta_{1}>0\\}$

denote the upper unit half-sphere.

We use the classical change of variable in Fowler [18],

$V(s,\theta)=r^{\frac{2\sigma}{p-1}}U(r,\theta),~{}~{}~{}~{}~{}s=\ln r.$

Direct calculations show that $V$ satisfies

$\begin{cases}V_{ss}-J_{1}V_{s}-J_{2}V+\theta_{1}^{2\sigma-1}\text{div}_{\theta}(\theta_{1}^{1-2\sigma}\nabla_{\theta}V)=0~{}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{in}~{}(-\infty,0)\times\mathbb{S}^{n}_{+},\\\
-\lim_{\theta_{1}\rightarrow
0^{+}}\theta_{1}^{1-2\sigma}\partial_{\theta_{1}}V=V^{p}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{on}~{}(-\infty,0)\times\partial\mathbb{S}^{n}_{+},\\\
\end{cases}$ (3.10)

where

$J_{1}=\frac{4\sigma}{p-1}-(N-2\sigma),~{}~{}~{}~{}~{}J_{2}=\frac{2\sigma}{p-1}\left(n-2\sigma-\frac{2\sigma}{p-1}\right).$

Multiplying (3.10) by $V_{s}$ and integrating, we have

$\displaystyle\int_{\mathbb{S}^{n}_{+}}\theta_{1}^{1-2\sigma}V_{ss}V_{s}-J_{2}\int_{\mathbb{S}^{n}_{+}}\theta_{1}^{1-2\sigma}VV_{s}-\int_{\mathbb{S}^{n}_{+}}\theta_{1}^{1-2\sigma}\nabla_{\theta}V\cdot\nabla_{\theta}V_{s}+\int_{\partial\mathbb{S}^{n}_{+}}V^{p}V_{s}$
(3.11)
$\displaystyle=J_{1}\int_{\mathbb{S}^{n}_{+}}\theta_{1}^{1-2\sigma}(V_{s})^{2}.$

For any $s\in(-\infty,0)$, we define

$\displaystyle\widetilde{E}(s):=$
$\displaystyle\frac{1}{2}\int_{\mathbb{S}^{n}_{+}}\theta_{1}^{1-2\sigma}(V_{s})^{2}-\frac{J_{2}}{2}\int_{\mathbb{S}^{n}_{+}}\theta_{1}^{1-2\sigma}V^{2}-\frac{1}{2}\int_{\mathbb{S}^{n}_{+}}\theta_{1}^{1-2\sigma}|\nabla_{\theta}V|^{2}$
$\displaystyle+\frac{1}{p+1}\int_{\partial\mathbb{S}^{n}_{+}}V^{p+1}.$

Then, by (3.11), we get

$\frac{d}{ds}\widetilde{E}(s)=J_{1}\int_{\mathbb{S}^{n}_{+}}\theta_{1}^{1-2\sigma}(V_{s})^{2}\geq
0.$ (3.12)

Here we have used the fact $J_{1}>0$ because
$1<p<\frac{n+2\sigma}{n-2\sigma}$. Hence, $\widetilde{E}(s)$ is non-decreasing
in $s\in(-\infty,0)$.

Now, rescaling back, we have

$\displaystyle\int_{\mathbb{S}^{n}_{+}}\theta_{1}^{1-2\sigma}(V_{s})^{2}$
$\displaystyle=\int_{\mathbb{S}^{n}_{+}}\theta_{1}^{1-2\sigma}\left(\frac{2\sigma}{p-1}r^{\frac{2\sigma}{p-1}-1}U+r^{\frac{2\sigma}{p-1}}U_{r}\right)^{2}r^{2}$
$\displaystyle=r^{2\frac{(p+1)\sigma}{p-1}-n}\int_{\partial^{+}\mathcal{B}_{r}^{+}}t^{1-2\sigma}\left(\frac{4\sigma^{2}}{(p-1)^{2}}r^{-1}U^{2}+\frac{4\sigma}{p-1}U\frac{\partial
U}{\partial\nu}+r|\frac{\partial U}{\partial\nu}|^{2}\right),$
$\int_{\mathbb{S}^{n}_{+}}\theta_{1}^{1-2\sigma}|\nabla_{\theta}V|^{2}=r^{2\frac{(p+1)\sigma}{p-1}-n+1}\int_{\partial^{+}\mathcal{B}_{r}^{+}}t^{1-2\sigma}\left(|\nabla
U|^{2}-|\frac{\partial U}{\partial\nu}|^{2}\right),$
$\int_{\mathbb{S}^{n}_{+}}\theta_{1}^{1-2\sigma}V^{2}=r^{2\frac{(p+1)\sigma}{p-1}-n-1}\int_{\partial^{+}\mathcal{B}_{r}^{+}}t^{1-2\sigma}U^{2},$
$\int_{\partial\mathbb{S}^{n}_{+}}V^{p+1}=r^{2\frac{(p+1)\sigma}{p-1}-n+1}\int_{\partial
B_{r}}u^{p+1}.$

Substituting these into (3.12) and noting that $s=\ln r$ is non-decreasing in
$r$, we easily obtain that $E(r;U)$ is also non-decreasing in $r\in(0,1)$. ∎

By the monotonicity of $E(r;U)$ we prove the following proposition, which will
play an essential role in deriving the lower bound in (1.5).

######  Proposition 3.3.

Let $U$ be a nonnegative solution of (1.3) with
$\frac{n}{n-2\sigma}<p<\frac{n+2\sigma}{n-2\sigma}$. If

$\liminf_{|x|\rightarrow 0}|x|^{\frac{2\sigma}{p-1}}u(x)=0,$

then

$\lim_{|x|\rightarrow 0}|x|^{\frac{2\sigma}{p-1}}u(x)=0.$

###### Proof.

Suppose by contradiction that

$\liminf_{|x|\rightarrow
0}|x|^{\frac{2\sigma}{p-1}}u(x)=0~{}~{}~{}~{}\textmd{and}~{}~{}~{}~{}\limsup_{|x|\rightarrow
0}|x|^{\frac{2\sigma}{p-1}}u(x)=C>0.$

Therefore, there exist two sequences of points $\\{x_{i}\\}$ and $\\{y_{i}\\}$
satisfying

$x_{i}\rightarrow 0,~{}~{}~{}y_{i}\rightarrow
0~{}~{}~{}~{}\textmd{as}~{}i\rightarrow\infty,$

such that

$|x_{i}|^{\frac{2\sigma}{p-1}}u(x_{i})\rightarrow
0~{}~{}~{}\textmd{and}~{}~{}~{}|y_{i}|^{\frac{2\sigma}{p-1}}u(y_{i})\rightarrow
C>0~{}~{}~{}\textmd{as}~{}i\rightarrow\infty.$

Let $g(r)=r^{\frac{2\sigma}{p-1}}\bar{u}(r)$, where
$\bar{u}(r)=\frac{1}{|\partial B_{r}|}\int_{\partial B_{r}}u$ denotes the
spherical average of $u$ over $\partial B_{r}$. Then, by the Harnack
inequality (3.3), we have

$\liminf_{r\rightarrow
0}g(r)=0~{}~{}~{}~{}\textmd{and}~{}~{}~{}~{}\limsup_{r\rightarrow 0}g(r)=C>0.$

Hence, there exists a sequence of local minimum points $r_{i}$ of $g(r)$ with

$\lim_{i\rightarrow\infty}r_{i}=0~{}~{}~{}~{}\textmd{and}~{}~{}~{}~{}\lim_{i\rightarrow\infty}g(r_{i})=0.$

Let

$V_{i}(X)=\frac{U(r_{i}X)}{U(r_{i}e_{1})},$

where $e_{1}=(1,0,\cdots,0)$. It follows from the Harnack inequality (3.2)
that $V_{i}$ is locally uniformly bounded away from the origin and satisfies

$\begin{cases}-\textmd{div}(t^{1-2\sigma}\nabla
V_{i})=0~{}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{in}~{}\mathbb{R}^{n+1}_{+},\\\
\frac{\partial
V_{i}}{\partial\nu^{\sigma}}(x,0)=\left(r_{i}^{\frac{2\sigma}{p-1}}U(r_{i}e_{1})\right)^{p-1}V_{i}^{p}(x,0)~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{on}~{}\mathbb{R}^{n}\backslash\\{0\\}.\end{cases}$
(3.13)

Note that by the Harnack inequality (3.3),
$r_{i}^{\frac{2\sigma}{p-1}}U(r_{i}e_{1})\rightarrow 0$ as
$i\rightarrow\infty$. Then by Corollary 2.10 and Theorem 2.15 in [23] that
there exists $\alpha>0$ such that for every $R>1>r>0$

$\|V_{i}\|_{W^{1,2}(t^{1-2\sigma},\mathcal{B}_{R}^{+}\backslash\overline{\mathcal{B}}_{r}^{+})}+\|V_{i}\|_{C^{\alpha}(\mathcal{B}_{R}^{+}\backslash\overline{\mathcal{B}}_{r}^{+})}+\|v_{i}\|_{C^{2,\alpha}(B_{R}\backslash
B_{r})}\leq C(R,r),$

where $v_{i}(x)=V_{i}(x,0)$ and $C(R,r)$ is independent of $i$. Then after
passing to a subsequence, $\\{V_{i}\\}$ converges to a nonnegative function
$V\in
W_{loc}^{1,2}(t^{1-2\sigma},\overline{\mathbb{R}^{n+1}_{+}}\backslash\\{0\\})\cap
C^{\alpha}_{loc}(\overline{\mathbb{R}^{n+1}_{+}}\backslash\\{0\\})$ satisfying

$\begin{cases}-\textmd{div}(t^{1-2\sigma}\nabla
V)=0~{}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{in}~{}\mathbb{R}^{n+1}_{+},\\\
\frac{\partial
V}{\partial\nu^{\sigma}}(x,0)=0~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{on}~{}\mathbb{R}^{n}\backslash\\{0\\}.\end{cases}$
(3.14)

By a Bôcher type theorem in [23], we have

$V(X)=\frac{a}{|X|^{n-2\sigma}}+b,$

where $a,b$ are nonnegative constants. Recall that $r_{i}$ are local minimum
of $g(r)$ for every $i$ and note that

$r^{\frac{2\sigma}{p-1}}\bar{v}_{i}(r)=r^{\frac{2\sigma}{p-1}}\frac{1}{|\partial
B_{r}|}\int_{\partial
B_{r}}v_{i}=\frac{1}{U(r_{i}e_{1})}r^{\frac{2\sigma}{p-1}}\bar{u}(r_{i}r)=\frac{1}{U(r_{i}e_{1})r_{i}^{\frac{2\sigma}{p-1}}}g(r_{i}r).$

Hence, for every $i$, we have

$\frac{d}{dr}\left[r^{\frac{2\sigma}{p-1}}\bar{v}_{i}(r)\right]\Bigg{|}_{r=1}=\frac{r_{i}}{U(r_{i}e_{1})r_{i}^{\frac{2\sigma}{p-1}}}g^{\prime}(r_{i})=0.$
(3.15)

Let $v(x)=V(x,0)$. Then we know that $v_{i}(x)\rightarrow v(x)$ in
$C_{loc}^{2}(\mathbb{R}^{n}\backslash\\{0\\})$. By (3.15), we obtain

$\frac{d}{dr}\left[r^{\frac{2\sigma}{p-1}}\bar{v}(r)\right]\Bigg{|}_{r=1}=0,$

which implies that

$a\left(\frac{2\sigma}{p-1}-(n-2\sigma)\right)+\frac{2\sigma b}{p-1}=0.$
(3.16)

On the other hand, by $V(e_{1})=1$, we have

$a+b=1.$ (3.17)

Combine (3.16) with (3.17), we get

$a=\frac{2\sigma}{(p-1)(n-2\sigma)}~{}~{}~{}~{}~{}\textmd{and}~{}~{}~{}~{}~{}b=1-\frac{2\sigma}{(p-1)(n-2\sigma)}.$

Since $\frac{n}{n-2\sigma}<p$, we have $0<a,b<1$. Now we compute $E(r;U)$.

It follows from Proposition 2.19 in [23] that $|\nabla_{x}V_{i}|$ and
$|t^{1-2\sigma}\partial_{t}V_{i}|$ are locally uniformly bounded in
$C_{loc}^{\beta}(\overline{\mathbb{R}^{n+1}_{+}}\backslash\\{0\\})$ for some
$\beta>0$. Hence, there exists a constant $C>0$ such that

$|\nabla_{x}U(X)|\leq
Cr_{i}^{-1}U(r_{i}e_{1})=o(1)r_{i}^{-\frac{2\sigma}{p-1}-1}~{}~{}~{}~{}~{}~{}\textmd{for}~{}\textmd{all}~{}|X|=r_{i}$

and

$|t^{1-2\sigma}\partial_{t}U(X)|\leq
Cr_{i}^{-2\sigma}U(r_{i}e_{1})=o(1)r_{i}^{-\frac{2\sigma}{p-1}-2\sigma}~{}~{}~{}~{}~{}~{}\textmd{for}~{}\textmd{all}~{}|X|=r_{i}.$

Thus, by a direct computation, we can get

$\lim_{i\rightarrow\infty}E(r_{i};U)=0.$

By the monotonicity of $E(r;U)$, we obtain

$E(r;U)\geq 0~{}~{}~{}~{}~{}\textmd{for}~{}\textmd{all}~{}r\in(0,1).$

On the other hand, by the scaling invariance of $E(r;U)$, for every $i$, we
have

$0\leq
E(r_{i};U)=E\left(1;r_{i}^{\frac{2\sigma}{p-1}}U(r_{i}X)\right)=E\left(1;r_{i}^{\frac{2\sigma}{p-1}}U(r_{i}e_{1})V_{i}\right).$

Hence, for every $i$, we have

$\displaystyle 0\leq$
$\displaystyle\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}|\frac{\partial
V_{i}}{\partial\nu}|^{2}+\frac{2\sigma}{p-1}\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}\frac{\partial
V_{i}}{\partial\nu}V_{i}$
$\displaystyle+\frac{1}{2}\frac{2\sigma}{p-1}\left(\frac{4\sigma}{p-1}-(n-2\sigma)\right)\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}V_{i}^{2}$
$\displaystyle-\frac{1}{2}\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}|\nabla
V_{i}|^{2}+\frac{1}{p+1}\int_{\partial
B_{1}}\left(r_{i}^{\frac{2\sigma}{p-1}}U(r_{i}e_{1})\right)^{p-1}V_{i}^{p+1}.$

Letting $i\rightarrow\infty$, we obtain

$\displaystyle 0$
$\displaystyle\leq\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}|\frac{\partial
V}{\partial\nu}|^{2}+\frac{2\sigma}{p-1}\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}\frac{\partial
V}{\partial\nu}V$
$\displaystyle~{}~{}~{}~{}+\frac{1}{2}\frac{2\sigma}{p-1}\left(\frac{4\sigma}{p-1}-(n-2\sigma)\right)\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}V^{2}-\frac{1}{2}\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}|\nabla
V|^{2}$
$\displaystyle=a^{2}(n-2\sigma)^{2}\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}-a(n-2\sigma)\frac{2\sigma}{p-1}\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}$
$\displaystyle~{}~{}~{}~{}+\frac{1}{2}\frac{2\sigma}{p-1}\left(\frac{4\sigma}{p-1}-(n-2\sigma)\right)\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}-\frac{1}{2}a^{2}(n-2\sigma)^{2}\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}$
$\displaystyle=\frac{\sigma}{p-1}\left(\frac{2\sigma}{p-1}-(n-2\sigma)\right)\int_{\partial^{+}\mathcal{B}_{1}^{+}}t^{1-2\sigma}<0.$

Note that in the last inequality we have used the fact
$\frac{2\sigma}{p-1}-(n-2\sigma)<0$ because $\frac{n}{n-2\sigma}<p$.
Obviously, we get a contradiction. ∎

To characterize the ”order” of an isolated singularity we establish the
following sufficient condition for removability of isolated singularities. For
its proof, we adapt the arguments from [19], but there are extra difficulties.
Such as, we need extra efforts to derive the estimates of $U$ from its trace
$u$.

######  Proposition 3.4.

Let $U$ be a nonnegative solution of (1.3) with
$\frac{n}{n-2\sigma}<p<\frac{n+2\sigma}{n-2\sigma}$. If

$\int_{\epsilon\leq|x|\leq 1}u^{\frac{(p-1)n}{2\sigma}}\leq c<+\infty$ (3.18)

with $c$ independent of $\epsilon$, then the singularity at the origin is
removable, i.e., u(x) can be extended to a continuous solution in the entire
ball $B_{1}$.

###### Proof.

Let

$p_{0}=\frac{n-2\sigma}{2\sigma}\left(p-\frac{n}{n-2\sigma}\right)~{}~{}~{}~{}\textmd{and}~{}~{}~{}~{}q_{0}=\frac{1}{2}(p_{0}+1)=\frac{1}{2}\frac{n-2\sigma}{2\sigma}(p-1).$

Following Serrin [28], we define, for $q\geq q_{0}$, $l>0$

$F(u)=\begin{cases}u^{q}~{}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{for}~{}0<u\leq
l,\\\
\frac{1}{q_{0}}\left[ql^{q-q_{0}}u^{q_{0}}+(q_{0}-q)l^{q}\right]~{}~{}~{}&\textmd{for}~{}l\leq
u,\end{cases}$

and

$G(u)=F(u)F^{\prime}(u)-q.$

Clearly, $F$ is a $C^{1}$ function of $u$ and $G$ is a piecewise smooth
function of $u$ with a corner at $u=l$. Moreover, since $p_{0}>0$ (implied by
$\frac{n}{n-2\sigma}<p$), $F,G$ satisfy

$F\leq\frac{q}{q_{0}}l^{q-q_{0}}u^{q_{0}},~{}~{}~{}~{}~{}~{}~{}~{}uF^{\prime}\leq
qF,$ (3.19) $|G|\leq FF^{\prime},$ (3.20)
$G^{\prime}\geq\begin{cases}\frac{1}{q}p{F^{\prime}}^{2}~{}~{}~{}~{}~{}&\textmd{for}~{}0<u<l,\\\
\frac{1}{q_{0}}p_{0}{F^{\prime}}^{2}~{}~{}~{}~{}~{}&\textmd{for}~{}l\leq
u.\end{cases}$ (3.21)

For any $0<R<1$, let $\eta$ and $\bar{\eta}$ be nonnegative $C^{\infty}$
function with $0\leq\eta,\bar{\eta}\leq 1$ in
$\mathcal{B}_{R}=\\{X\in\mathbb{R}^{n+1}:|X|<R\\}$, $\eta$ having compact
support in $\mathcal{B}_{R}$, and $\bar{\eta}$ vanishing in some neighborhood
of the origin. Take $(\eta\bar{\eta})^{2}G(U)$ as a test function into (1.4),
we have

$\displaystyle\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}(\eta\bar{\eta})^{2}G^{\prime}(U)|\nabla
U|^{2}+2\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}\eta\bar{\eta}G(U)\nabla
U\cdot\nabla(\eta\bar{\eta})$ (3.22)
$\displaystyle=\int_{B_{R}}(\eta\bar{\eta})^{2}G(u)u^{p}.$

Using (3.19) – (3.21) and simplifying we obtain from (3.22)

$\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}(\eta\bar{\eta})^{2}|\nabla(F(U))|^{2}\leq
C(q)\bigg{\\{}\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla(\eta\bar{\eta})|^{2}F^{2}+\int_{B_{R}}(\eta\bar{\eta})^{2}u^{p-1}F^{2}\bigg{\\}}.$
(3.23)

By Hölder and Proposition 2.3, we have

$\displaystyle\int_{B_{R}}(\eta\bar{\eta})^{2}u^{p-1}F^{2}$
$\displaystyle\leq\left(\int_{B_{R}}(u^{p-1})^{\frac{n}{2\sigma}}\right)^{\frac{2\sigma}{n}}\left(\int_{B_{R}}(\eta\bar{\eta}F)^{\frac{2n}{n-2\sigma}}\right)^{\frac{n-2\sigma}{n}}$
$\displaystyle\leq
C_{n,\sigma}^{2}\left(\int_{B_{R}}(u^{p-1})^{\frac{n}{2\sigma}}\right)^{\frac{2\sigma}{n}}\left(\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla(\eta\bar{\eta}F)|^{2}\right)$
$\displaystyle\leq
C_{n,\sigma}^{2}\|u^{p-1}\|_{L^{\frac{n}{2\sigma}}(B_{R})}\bigg{\\{}\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}(\eta\bar{\eta})^{2}|\nabla
F|^{2}$
$\displaystyle~{}~{}~{}~{}+\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla(\eta\bar{\eta})|^{2}F^{2}\bigg{\\}}.$

By the assumption (3.18), we can choose $R$ small enough (depending on $q,n$
and $\sigma$ ) such that

$\|u^{p-1}\|_{L^{\frac{n}{2\sigma}}(B_{R})}\leq\frac{1}{2}\frac{1}{C(q)C_{n,\sigma}^{2}},$

Hence, from (3.23), we obtain

$\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}(\eta\bar{\eta})^{2}|\nabla(F(U))|^{2}\leq
C(q)\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla(\eta\bar{\eta})|^{2}F^{2}$
(3.24)

with a new constant $C(q)$. Therefore, we have from (3.24)

$\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}(\eta\bar{\eta})^{2}|\nabla(F(U))|^{2}\leq
C(q)\bigg{\\{}\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla\eta|^{2}F^{2}+\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla\bar{\eta}|^{2}F^{2}\bigg{\\}}$
(3.25)

and

$\left(\int_{B_{R}}(\eta\bar{\eta}F)^{\frac{2n}{n-2\sigma}}\right)^{\frac{n-2\sigma}{n}}\leq
C(q)\bigg{\\{}\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla\eta|^{2}F^{2}+\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla\bar{\eta}|^{2}F^{2}\bigg{\\}}.$
(3.26)

For any $\epsilon>0$ small enough, we choose $\bar{\eta}_{\epsilon}$ satisfy

$\bar{\eta}_{\epsilon}(X)=\begin{cases}0~{}~{}~{}~{}~{}\textmd{for}~{}|X|\leq\epsilon,\\\
1~{}~{}~{}~{}~{}\textmd{for}~{}2\epsilon\leq|X|<R,\end{cases}$ (3.27)

and $|\nabla\bar{\eta}_{\epsilon}(X)|\leq\frac{c}{\epsilon}$ for all
$X\in\mathcal{B}_{R}$. By Hölder inequality

$\displaystyle\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla\bar{\eta}_{\epsilon}|^{2}F^{2}$
$\displaystyle\leq\left(\int_{\mathcal{B}_{2\epsilon}^{+}}t^{1-2\sigma}|\nabla\bar{\eta}_{\epsilon}|^{n+2-2\sigma}\right)^{\frac{2}{n+2-2\sigma}}\left(\int_{\mathcal{B}_{2\epsilon}^{+}}t^{1-2\sigma}F^{2\frac{n+2-2\sigma}{n-2\sigma}}\right)^{\frac{n-2\sigma}{n+2-2\sigma}}$
(3.28)
$\displaystyle\leq\frac{C}{\epsilon^{2}}\left(\int_{\mathcal{B}_{2\epsilon}^{+}}t^{1-2\sigma}\right)^{\frac{2}{n+2-2\sigma}}\left(\int_{\mathcal{B}_{2\epsilon}^{+}}t^{1-2\sigma}F^{2\frac{n+2-2\sigma}{n-2\sigma}}\right)^{\frac{n-2\sigma}{n+2-2\sigma}}$
$\displaystyle\leq
C\left(\int_{\mathcal{B}_{2\epsilon}^{+}}t^{1-2\sigma}F^{2\frac{n+2-2\sigma}{n-2\sigma}}\right)^{\frac{n-2\sigma}{n+2-2\sigma}},$

where $C$ is a positive constant independent of $\epsilon$. Since, by (3.19),
the Harnack inequality (3.3) and (3.1), we have

$\displaystyle\int_{\mathcal{B}_{2\epsilon}^{+}}t^{1-2\sigma}F^{2\frac{n+2-2\sigma}{n-2\sigma}}$
$\displaystyle\leq
C(l,q)\int_{\mathcal{B}_{2\epsilon}^{+}}t^{1-2\sigma}U^{2q_{0}\frac{n+2-2\sigma}{n-2\sigma}}$
(3.29) $\displaystyle\leq
C(l,q)\int_{0}^{2\epsilon}\left(\sup_{\partial^{+}\mathcal{B}_{s}^{+}}U\right)^{2q_{0}\frac{n+2-2\sigma}{n-2\sigma}}s^{n+1-2\sigma}ds$
$\displaystyle\leq
C(l,q)\int_{0}^{2\epsilon}\left(\inf_{\partial^{+}\mathcal{B}_{s}^{+}}U\right)^{2q_{0}\frac{n}{n-2\sigma}}s^{-2q_{0}\frac{4\sigma(1-\sigma)}{(p-1)(n-2\sigma)}}s^{n+1-2\sigma}ds$
$\displaystyle\leq C(l,q)\int_{0}^{2\epsilon}\left(\inf_{\partial
B_{s}}u\right)^{2q_{0}\frac{n}{n-2\sigma}}s^{n-1}ds$ $\displaystyle\leq
C(l,q)\int_{0}^{2\epsilon}\left(\frac{1}{|\partial B_{s}|}\int_{\partial
B_{s}}u^{2q_{0}\frac{n}{n-2\sigma}}\right)s^{n-1}ds$ $\displaystyle\leq
C(l,q)\int_{B_{2\epsilon}}u^{2q_{0}\frac{n}{n-2\sigma}}=C(l,q)\int_{B_{2\epsilon}}u^{\frac{(p-1)n}{2\sigma}}.$

Here we have used the fact
$2q_{0}\frac{n}{n-2\sigma}=\frac{(p-1)n}{2\sigma}>1$ because
$p-1>\frac{2\sigma}{n-2\sigma}$. Now, we have from (3.18), (3.28) and (3.29)

$\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla\bar{\eta}_{\epsilon}|^{2}F^{2}\rightarrow
0~{}~{}~{}~{}\textmd{as}~{}\epsilon\rightarrow 0.$

This together with (3.25) and (3.26), we obtain

$\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}\eta^{2}|\nabla(F(U))|^{2}\leq
C(q)\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla\eta|^{2}F^{2}$ (3.30)

and

$\left(\int_{B_{R}}(\eta
F)^{\frac{2n}{n-2\sigma}}\right)^{\frac{n-2\sigma}{n}}\leq
C(q)\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla\eta|^{2}F^{2}.$ (3.31)

Let $l\rightarrow+\infty$, since $F(u)\rightarrow u^{q}$, we have

$\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}\eta^{2}|\nabla U^{q}|^{2}\leq
C(q)\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla\eta|^{2}U^{2q}$ (3.32)

and

$\left(\int_{B_{R}}(\eta
u^{q})^{\frac{2n}{n-2\sigma}}\right)^{\frac{n-2\sigma}{n}}\leq
C(q)\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla\eta|^{2}U^{2q}.$ (3.33)

By a similar estimate as in (3.29), we can obtain

$\displaystyle\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla\eta|^{2}U^{2q}$
$\displaystyle\leq
C(q)\int_{0}^{R}\left(\inf_{\partial^{+}\mathcal{B}_{s}^{+}}U^{2q-2q_{0}\frac{2(1-\sigma)}{n-2\sigma}}\right)s^{n-1}ds$
(3.34) $\displaystyle\leq
C(q)\left(\int_{B_{R}}u^{2q}\right)^{1-\frac{q_{0}}{q}\frac{2(1-\sigma)}{n-2\sigma}}$

Inequality (3.32), (3.33) and (3.34) can be iterated a finite number of times
to show that

$U\in
W^{1,2}(t^{1-2\sigma},\mathcal{B}_{R}^{+})~{}~{}~{}~{}\textmd{and}~{}~{}~{}~{}u\in
L^{q}(B_{R})~{}~{}~{}~{}\textmd{for}~{}\textmd{all}~{}q>0.$

Furthermore, $U$ satisfies

$\begin{cases}-\textmd{div}(t^{1-2\sigma}\nabla
U)=0~{}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{in}~{}\mathcal{B}_{R}^{+},\\\
\frac{\partial
U}{\partial\nu^{\sigma}}(x,0)=U^{p}(x,0)~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{on}~{}\partial^{0}\mathcal{B}_{R}^{+}.\end{cases}$

Indeed, for $\epsilon>0$ small, let $\bar{\eta}_{\epsilon}$ be a smooth cut-
off function as in (3.27). Let $\psi\in
C_{c}^{\infty}(\mathcal{B}_{R}^{+}\cup\partial^{0}\mathcal{B}_{R}^{+})$. It
follows from (1.4) that

$\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}\nabla
U\nabla(\psi\bar{\eta}_{\epsilon})=\int_{\partial^{0}\mathcal{B}_{R}^{+}}U^{p}\psi\bar{\eta}_{\epsilon}.$
(3.35)

Since

$\left|\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}\nabla
U\nabla\bar{\eta}_{\epsilon}\psi\right|\leq
C\epsilon^{\frac{n-2\sigma}{2}}\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}|\nabla
U|^{2}\rightarrow 0~{}~{}~{}~{}\textmd{as}~{}\epsilon\rightarrow 0,$

by the dominated convergence theorem, let $\epsilon\rightarrow 0$ in (3.35),
we obtain

$\int_{\mathcal{B}_{R}^{+}}t^{1-2\sigma}\nabla
U\nabla\psi=\int_{\partial^{0}\mathcal{B}_{R}^{+}}U^{p}\psi.$

Since $u\in L^{q}(B_{R})$ for some $q>\frac{n}{2\sigma}$, it follows from
Proposition 2.10 in [23] that $U$ is Hölder continuous in
$\overline{\mathcal{B}_{R/2}^{+}}$. The proof of the proposition is completed.
∎

######  Corollary 3.1.

Let $U$ be a nonnegative solution of (1.3) with
$\frac{n}{n-2\sigma}<p<\frac{n+2\sigma}{n-2\sigma}$. Then either the origin 0
is a removable singularity or $\lim_{|X|\rightarrow 0}U(x,t)=+\infty$.

###### Proof.

By Proposition 3.4, if the origin is not a removable singularity, then there
exists a sequence of points $\\{x_{j}\\}$ such that

$r_{j}=|x_{j}|\rightarrow
0~{}~{}~{}~{}\textmd{and}~{}~{}~{}~{}U(x_{j},0)\rightarrow+\infty~{}~{}~{}\textmd{as}~{}j\rightarrow\infty.$

By the Harnack inequality (3.2), we have

$\inf_{|X|=r_{j}}U(X)\geq C^{-1}U(x_{j},0).$

By the maximum principle,

$U(X)\geq\inf_{|X|=r_{j},r_{j+1}}U(X)\geq
C^{-1}\min(U(x_{j},0),U(x_{j+1},0))~{}~{}~{}\textmd{in}~{}r_{j+1}\leq|X|\leq
r_{j}.$

Hence, we have $U(X)\rightarrow+\infty$ as $|X|\rightarrow 0$. ∎

Proof of Theorem 1.1. By Proposition 3.1,

$u(x)\leq c|x|^{-\frac{2\sigma}{p-1}}.$

If (1.5) does not hold, then

$\liminf_{x\rightarrow 0}|x|^{\frac{2\sigma}{p-1}}u(x)=0.$

It follows Proposition 3.3 that

$\lim_{x\rightarrow 0}|x|^{\frac{2\sigma}{p-1}}u(x)=0$ (3.36)

We will prove that the origin is a removable singularity. It suffices to
establish (3.18) by Proposition 3.4.

Let $\tau=\frac{n-2\sigma}{p-1}(p-\frac{n}{n-2\sigma})$ and $0<\delta<1$.
Define

$\Phi=|X|^{-\tau}-\delta t^{2\sigma}|X|^{-(\tau+2\sigma)}.$

Then we can choose $\delta$ small (depending only on $n$ , $\sigma$ and $p$)
such that

$\begin{cases}-\textmd{div}(t^{1-2\sigma}\nabla\Phi)\geq
0~{}~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{in}~{}\mathbb{R}^{n+1}_{+},\\\
\frac{\partial\Phi}{\partial\nu^{\sigma}}(x,0)=2\delta\sigma|x|^{-2\sigma}\phi(x)~{}~{}~{}~{}~{}~{}~{}~{}&\textmd{on}~{}\mathbb{R}^{n}\backslash\\{0\\},\end{cases}$
(3.37)

where $\phi(x)=\Phi(x,0)=|x|^{-\tau}$. Take $\xi_{1}(s)\in
C_{c}^{\infty}(\mathbb{R})$ satisfying $0\leq\xi_{1}\leq 1$ in $\mathbb{R}$
and

$\xi_{1}(s)=\begin{cases}0~{}~{}~{}~{}~{}~{}&\textmd{if}~{}|s|\leq 1,\\\
1~{}~{}~{}~{}~{}~{}&\textmd{if}~{}|s|\geq 2.\end{cases}$

Take $\xi_{2}(s)\in C_{c}^{\infty}(\mathbb{R})$ satisfying $0\leq\xi_{2}\leq
1$ in $\mathbb{R}$ and

$\xi_{2}(s)=\begin{cases}1~{}~{}~{}~{}~{}~{}&\textmd{if}~{}|s|\leq\frac{1}{2},\\\
0~{}~{}~{}~{}~{}~{}&\textmd{if}~{}|s|\geq\frac{3}{4}.\end{cases}$

For any $\epsilon>0$ small, we choose $\zeta(X)\in
C_{c}^{\infty}(\epsilon<|X|<1)$ as follow

$\zeta(X)=\begin{cases}\xi_{1}(\frac{|X|}{\epsilon})~{}~{}~{}~{}~{}~{}&\textmd{if}~{}|X|\leq\frac{1}{2},\\\
\xi_{2}(|X|)~{}~{}~{}~{}~{}~{}&\textmd{if}~{}|X|\geq\frac{1}{2}.\end{cases}$

Using $\zeta\Phi$ as a test function in (1.4), and the divergence theorem we
obtain

$\displaystyle\int_{B_{1}}u\zeta\phi\left(\frac{1}{\phi}\frac{\partial\Phi}{\partial\nu^{\sigma}}-u^{p-1}\right)=$
$\displaystyle\int_{\mathcal{B}_{1}^{+}}U\bigg{\\{}\textmd{div}(t^{1-2\sigma}\nabla\Phi)\zeta+2t^{1-2\sigma}\nabla\zeta\cdot\nabla\Phi$
(3.38)
$\displaystyle+\textmd{div}(t^{1-2\sigma}\nabla\zeta)\Phi\bigg{\\}}-\int_{B_{1}}u\phi\frac{\partial\zeta}{\partial\nu^{\sigma}}.$

Since $\zeta$ is radial, we have

$\frac{\partial\zeta}{\partial\nu^{\sigma}}(x,0)=-\lim_{t\rightarrow
0^{+}}t^{1-2\sigma}\partial_{t}\zeta=0~{}~{}~{}~{}~{}\textmd{in}~{}~{}B_{1}.$
(3.39)

By (3.37), we obtain

$\int_{\mathcal{B}_{1}^{+}}U\textmd{div}(t^{1-2\sigma}\nabla\Phi)\zeta\leq 0.$
(3.40)

Using (3.1) and the Harnack inequality (3.2), we have

$\displaystyle\int_{\mathcal{B}_{1}^{+}}U\bigg{\\{}2t^{1-2\sigma}\nabla\zeta\cdot\nabla\Phi+\textmd{div}(t^{1-2\sigma}\nabla\zeta)\Phi\bigg{\\}}$
(3.41) $\displaystyle\leq
C_{1}+\int_{\mathcal{B}_{2\epsilon}^{+}\backslash\mathcal{B}_{\epsilon}^{+}}U\bigg{\\{}2t^{1-2\sigma}|\nabla\zeta||\nabla\Phi|+|\textmd{div}(t^{1-2\sigma}\nabla\zeta)|\Phi\bigg{\\}}$
$\displaystyle\leq
C_{1}+C\epsilon^{-\frac{2\sigma}{p-1}-1}\int_{\mathcal{B}_{2\epsilon}^{+}\backslash\mathcal{B}_{\epsilon}^{+}}\left(t^{1-2\sigma}|X|^{-(\tau+1)}+|X|^{-(\tau+2\sigma)}\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+C\epsilon^{-\frac{2\sigma}{p-1}-\tau}\int_{\mathcal{B}_{2\epsilon}^{+}\backslash\mathcal{B}_{\epsilon}^{+}}t^{1-2\sigma}\left(|\Delta\zeta|+\epsilon^{-1}\left|\xi_{1}^{\prime}(\frac{|X|}{\epsilon})\right|\frac{1}{|X|}\right)$
$\displaystyle\leq C_{1}+C\epsilon^{n-(\tau+2\sigma)-\frac{2\sigma}{p-1}}\leq
C<+\infty$

with $C$ independent of $\epsilon$. Hence, we obtain

$\int_{B_{1}}u\zeta\phi\left(\frac{1}{\phi}\frac{\partial\Phi}{\partial\nu^{\sigma}}-u^{p-1}\right)\leq
C<+\infty.$

By (3.36), we have $u^{p-1}=o(|x|^{-2\sigma})$, and while

$\frac{1}{\phi}\frac{\partial\Phi}{\partial\nu^{\sigma}}=2\delta\sigma|x|^{-2\sigma}~{}~{}~{}~{}~{}~{}\textmd{in}~{}B_{1}\backslash\\{0\\}.$

Therefore

$\int_{B_{1}}u\zeta|x|^{-(n-\frac{2\sigma}{p-1})}=\int_{B_{1}}u\zeta|x|^{-(\tau+2\sigma)}\leq
C<+\infty$ (3.42)

with $C$ independent of $\epsilon$. Again by (3.1), we have

$\displaystyle\int_{2\epsilon\leq|x|\leq 1/2}u^{\frac{(p-1)n}{2\sigma}}$
$\displaystyle\leq C\int_{2\epsilon\leq|x|\leq
1/2}u|x|^{-\frac{2\sigma}{p-1}[\frac{(p-1)n}{2\sigma}-1]}$ $\displaystyle\leq
C\int_{B_{1}}u\zeta|x|^{-(n-\frac{2\sigma}{p-1})}\leq C<+\infty$

with $C$ independent of $\epsilon$. Thus, we establish (3.18) and complete the
proof of Theorem 1.1. $\square$

## References

  * [1] W. Ao, H. Chan, A. DelaTorre, M.A. Fontelos, M. González, J. Wei, On higher dimensional singularities for the fractional Yamabe problem: a non-local Mazzeo-Pacard program, arXiv:1802.07973.
  * [2] W. Ao, A. DelaTorre, M. González, J. Wei, A gluing approach for the fractional Yamabe problem with prescribed isolated singularities, Preprint, arXiv:1609. 08903.
  * [3] W. Ao, H. Chan, M. González, J. Wei, Existence of positive weak solutions for fractional Lane-Emden equations with prescribed singular set. Preprint.
  * [4] P. Aviles, Local behavior of solutions of some elliptic equations, Comm. Math. Phys., 108 (1987) 177-192.
  * [5] M.-F. Bidaut-Véron, L. Véron, Nonlinear elliptic equations on compact Riemannian manifolds and asymptotics of Emden equations, Invent. Math., 106 (3) (1991) 489-539.
  * [6] X. Cabré, Y. Sire, Nonlinear equations for fractional Laplacians, I: Regularity, maximum principles, and Hamiltonian estimates, Ann. Inst. H. Poincaré Anal. Non Linéaire, 31 (2014) 23-53.
  * [7] L. Caffarelli, B. Gidas, J. Spruck, Asymptotic symmetry and local behavior of semilinear elliptic equations with critical Sobolev growth, Comm. Pure Appl. Math. , 42 (1989) 271-297.
  * [8] L. Caffarelli, T. Jin, Y. Sire, J. Xiong, Local analysis of solutions of fractional semi-linear elliptic equations with isolated singularities, Arch. Ration. Mech. Anal., 213 (1) (2014) 245-268.
  * [9] L. Caffarelli, L. Silvestre, An extension problem related to the fractional Laplacian, Comm. Partial Differential Equations, 32 (8) (2007) 1245-1260.
  * [10] S.-Y. A. Chang, M. González, Fractional Laplacian in conformal geometry, Adv. Math.,  226 (2011) 1410-1432.
  * [11] C.-C. Chen, C.-S. Lin, On the asymptotic symmetry of singular solutions of the scalar curvature equations, Math. Ann., 313 (1999) 229-245.
  * [12] W. Chen, C. Li, Y. Li, A direct method of moving planes for the fractional Laplacian, Adv. Math., 308 (2017) 404-437.
  * [13] J. Dávila, L. Dupaigne, J. Wei, On the fractional Lane-Emden equation, Trans. Amer. Math. Soc., 369 (2017) 6087-6104.
  * [14] A. DelaTorre, M.d.M. González, Isolated singularities for a semilinear equation for the fractional Laplacian arising in conformal geometry, Preprint.  To appear in Rev. Mat. Iber.
  * [15] A. DelaTorre, M. del Pino, M.d.M. González, J. Wei, Delaunay-type singular solutions for the fractional Yamabe problem, Math. Ann., 369 (2017) 597-626.
  * [16] E. Di Nezza, G. Palatucci, E. Valdinoci. Hitchhiker’s guide to the fractional Sobolev spaces, Bull. Sci. Math., 136 (2012) 521-573.
  * [17] M. Fall, V. Felli, Unique continuation property and local asymptotics of solutions to fractional elliptic equations, Comm. Partial Differential Equations, 39 (2014) 354-397.
  * [18] R. H. Fowler, Further studies of Emden’s and similar differential equations, Q. J. Math., Oxf. Ser. 2 (1931) 259-288.
  * [19] B. Gidas, J. Spruck, Global and local behavior of positive solutions of nonlinear elliptic equations, Comm. Pure Appl. Math., 34 (1981) 525-598.
  * [20] M.d.M. González, R. Mazzeo, Y. Sire, Singular solutions of fractional order conformal Laplacians, J. Geom. Anal., 22 (2012) 845-863.
  * [21] M.d.M. González, J. Qing, Fractional conformal Laplacians and fractional Yamabe problems, Anal. PDE,  6 (2013) 1535-1576.
  * [22] T. Jin, O.S. de Queiroz, Y. Sire, J. Xiong, On local behavior of singular positive solutions to nonlocal elliptic equations, Calc. Var. Partial Differential Equations, 56 (1) (2017) Art. 9, 25 pp.
  * [23] T. Jin, Y.Y. Li, J. Xiong, On a fractional Nirenberg problem, part I: blow up analysis and compactness of solutions, J. Eur. Math. Soc., 16 (2014) 1111-1171.
  * [24] S. Kim, M. Musso, J. Wei, Existence theorems of the fractional Yamabe problem, arXiv:1603.06617.
  * [25] N. Korevaar, R. Mazzeo, F. Pacard, R. Schoen, Refined asymptotics for constant scalar curvature metrics with isolated singularities, Invent. Math., 135 (2) (1999) 233-272.
  * [26] P.-L. Lions, Isolated singularities in semilinear problems, J. Differential Equations, 38 (1980) 441-450.
  * [27] P. Poláčik, P. Quittner, Ph. Souplet, Singularity and decay estimates in superlinear problems via Liouville-type theorems, Duke Math. J., 139 (2007) 555-579.
  * [28] J, Serrin, Local behavior of solutions of quasi-linear equation, Acta. Math., 111 (1964) 247-302.
  * [29] J, Serrin, Isolated singularities of solutions of quasi-linear equation, Acta. Math., 113 (1965) 219-240.