\section{Introduction}
Let $N\geq 2$, $s\in (0, 1)$ and $p\in (1, \frac{N}{s})$.
Let us consider the following nonlinear problem
\begin{equation}\label{P}
\left\{
\begin{array}{ll}
(-\Delta)^{s}_{p} u= g(u) \, \mbox{ in } \mathbb{R}^{N}, \\
u\in W^{s, p}(\mathbb{R}^{N}),
\end{array}
\right.
\end{equation}
where $(-\Delta)^{s}_{p}$ is the fractional $p$-Laplacian operator (see \cite{CL, DTGCV, IMS, P}) defined by
$$
(-\Delta)^{s}_{p}u(x)=C_{N, s, p} P.V. \int_{\mathbb{R}^{N}} \frac{|u(x)-u(y)|^{p-2}(u(x)-u(y))}{|x-y|^{N+sp}}dy, \quad
C_{N, s, p}=\frac{\frac{sp}{2}(1-s)2^{2s-1}}{\pi^{\frac{N-1}{2}}} \frac{\Gamma(\frac{N+sp}{2})}{\Gamma(\frac{p+1}{2})}\Gamma(2-s),
$$
for all $u\in C^{1, 1}_{loc}(\mathbb{R}^{N})\cap L^{p-1}_{sp}$, with
$$
L^{p-1}_{sp}=\Bigl\{u\in L^{p-1}_{loc}(\mathbb{R}^{N}): \int_{\mathbb{R}^{N}}\frac{|u(x)|^{p-1}}{1+|x|^{N+sp}}dx<\infty\Bigr\},
$$
and
$g:\mathbb{R}\rightarrow \mathbb{R}$ is a continuous nonlinearity such that $g(0)=0$, $g$ is odd and satisfies the following Berestycki-Lions type conditions \cite{BL}:
\begin{compactenum}[$(g1)$]
\item $\displaystyle{-\infty<\liminf_{t\rightarrow 0} \frac{g(t)}{|t|^{p-2}t}\leq \limsup_{t\rightarrow 0} \frac{g(t)}{|t|^{p-2}t}=-m}$, with $m>0$,
\item $\displaystyle{-\infty\leq \limsup_{t\rightarrow 0} \frac{g(t)}{|t|^{p^{*}_{s}-1}}\leq 0}$, where $p^{*}_{s}=\frac{Np}{N-sp}$ is the critical fractional exponent,
\item there exists $\zeta>0$ such that $G(\zeta)>0$, where $G(t)=\int_{0}^{t}g(\tau)d\tau$.
\end{compactenum}
Our first result deals with the existence of a weak solution to \eqref{P} satisfying a Pohozaev type identity. More precisely, we prove the following result.
\begin{thm}\label{thm1}
Let $N\geq 2$, $s\in (0, 1)$ and $p\in (1, \frac{N}{s})$. Assume that $(g1)$-$(g3)$ hold. Then there exists a nontrivial weak solution $u\in W^{s, p}(\mathbb{R}^{N})$ to \eqref{P} fulfilling the following Pohozaev identity:
\begin{align}\label{POHOZAEV}
\frac{C_{N, s, p}(N-sp)}{2p}\iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}}dxdy-N\int_{\mathbb{R}^{N}} G(u)dx=0.
\end{align}
\end{thm}
We recall that the Pohozaev identity for the $p$-Laplacian operator $\Delta_{p}u=\dive(|\nabla u|^{p-2}\nabla u)$ has been extensively investigated in literature; see for instance \cite{GV, PS}.
On the other hand, in the fractional setting, if $s\in (0, 1)$ and $p=2$, then $(-\Delta)^{s}_{p}$ becomes the fractional Laplacian operator $(-\Delta)^{s}$ and the corresponding Pohozaev identity for weak solutions to \eqref{P} has been established in \cite[Proposition 4.1]{CW} (see also \cite[Theorem 3.5.1]{A} for more details and \cite{RS} for the case of bounded domains). More precisely, in \cite{CW} the authors required $g\in C^{1}(\mathbb{R})$ and obtained the Pohozaev identity for $(-\Delta)^{s}$ by combining the extension method \cite{CS}, some regularity results for $(-\Delta)^{s}$,
an integration by parts and some ideas inspired by the proof given in \cite[Proposition 1]{BL}.
When $p\neq 2$, we are not able to adapt this approach. In fact, even if the extension method for $(-\Delta)^{s}_{p}$ has been recently studied in \cite{DTGCV}, we do not have enough informations about the regularity of the extension $U(x, y)$ of $u(x)$ in $\mathbb{R}^{N+1}_{+}=\{(x, y)\in \mathbb{R}^{N+1}: y>0 \}$ whenever $p\neq 2$, and contrary to the linear case $p=2$ we do not know if the following identity is valid:
$$
C\iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}}dxdy=\iint_{\mathbb{R}^{N+1}_{+}} y^{-1+p(1-s)} |\nabla U|^{p}dxdy,
$$
where $C>0$ is an appropriate constant.
Therefore, we need to follow a different strategy to arrive at \eqref{POHOZAEV}. We stress that the main difficulty in reaching \eqref{POHOZAEV}, it consists in using $x\cdot \nabla u$, where $u$ is a weak solution to \eqref{P}, as test function in the weak formulation of \eqref{P} and then apply an integration by parts. Indeed, due to the fact that for $p\neq 2$ we have to handle the nonlinearity of the operator $(-\Delta)^{s}_{p}$ and its nonlocal character, it seems to be an awkward task to verify that $x\cdot \nabla u\in W^{s, p}(\mathbb{R}^{N})$. Moreover, an integration by parts formula for $(-\Delta)^{s}_{p}$ seems to be hard to accomplish.
In this paper, by means of suitable variational methods for potential operators with covariance condition (see \cite{B, Chab}), we prove the existence of a weak solution to \eqref{P} and that satisfies \eqref{POHOZAEV}.
We recall that $u\in W^{s, p}(\mathbb{R}^{N})$ is a weak solution to \eqref{P} if  $\langle I'(u),\phi\rangle=0$ for all $\phi\in W^{s, p}(\mathbb{R}^{N})$, where $I:W^{s, p}(\mathbb{R}^{N})\rightarrow \mathbb{R}$ is  the energy functional associated with \eqref{P}, namely
$$
I(u)=\frac{C_{N, s, p}}{2p}\iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}}dxdy-\int_{\mathbb{R}^{N}} G(u)dx.
$$
Note that the restriction $N\geq 2$ ensures that \eqref{P} admits a radial weak solution (see Theorem \ref{THMBL}).
We stress that Theorem \ref{thm1} only guarantees the existence of a weak solution to \eqref{P} fulfilling \eqref{POHOZAEV} and not that, as in the case $p=2$, any weak solution to \eqref{P} satisfies \eqref{POHOZAEV}.
However, this last statement is true for $C^{1,1}$ functions satisfying equation \eqref{P} in the pointwise sense.
More precisely, our second main result can be stated as follows.
\begin{thm}\label{thm2}
Let $N\geq 2$, $s\in (0, 1)$ and $p\in (1, \frac{N}{s})$.
When $p\in (1, 2)$, we assume that $s<\frac{2(p-1)}{p}$. Assume that $(g1)$-$(g3)$ hold.
Let $u\in W^{s, p}(\mathbb{R}^{N})\cap C^{1, 1}(\mathbb{R}^{N})$ be such that $(-\Delta)^{s}_{p}u=g(u)$ in $\mathbb{R}^{N}$. Then \eqref{POHOZAEV} is valid.
\end{thm}
The proof of Theorem \ref{thm2} goes as follows. We first prove an integral representation formula for $(-\Delta)^{s}_{p}v$, with $v\in C^{1, 1}(\mathbb{R}^{N})$, in the spirit of \cite[Lemma 3.2]{DPV} (see Lemma \ref{RIF}). Thus, motivated by \cite[Lemma 4.2]{DFW1}, we establish an integration by parts formula for $C^{1, 1}$ functions and vector fields of class $C^{0, 1}_{c}$ (see Lemma \ref{LEMMA}).
Therefore, multiplying \eqref{P} by $\varphi(\lambda x) x\cdot \nabla u$,
where $\varphi\in C^{1}_{c}(\mathbb{R}^{N})$ is such that $\varphi=1$ in a neighborhood of $0$ and $\lambda>0$,  integrating over $\mathbb{R}^N$ and using our integration by parts formula, we can pass to the limit as $\lambda\rightarrow 0$ in the relation that we find to infer the desired result.

As far as we know, this is the first time that a Pohozaev identity for the fractional $p$-Laplacian operator $(-\Delta)^{s}_{p}$ in $\mathbb{R}^{N}$ is obtained in literature. \\
The paper is organized as follows. In section $2$ we give the proof of Theorem \ref{thm1}. Section $3$ is devoted to the proof of Theorem \ref{thm2}.

\subsection*{Notations:}
Henceforth, $\|\cdot\|_{L^{p}(\mathbb{R}^{N})}$ denotes the $L^{p}(\mathbb{R}^N)$-norm with $p\in [1, \infty]$, for $s\in (0, 1)$ and $p\in [1, \infty)$ the fractional Sobolev space
$W^{s, p}(\mathbb{R}^{N})$ is the set of functions $u\in L^{p}(\mathbb{R}^{N})$ such that $\iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}}dxdy<\infty$,
$D^{\beta} u:=\frac{\partial^{|\beta|}u}{{\partial x_{1}^{\beta_{1}}}\dots {\partial x_{N}^{\beta_{N}}}}$ stands for the distributional derivative of order $\beta=(\beta_{1}, \dots, \beta_{N})\in (\mathbb{N}\cup\{0\})^{N}$ ($\beta$ is called a multi-index), where $|\beta|:=\beta_{1}+\dots+\beta_{N}$ is the order of $\beta$. With $C^{k, \gamma}(\mathbb{R}^{N})$, where $k\in \mathbb{N}\cup \{0\}$ and $\gamma\in (0, 1]$, we denote the set of functions $u\in C^{k}(\mathbb{R}^N)$ whose derivatives up to order $k$ are bounded and such that
$$
 \|u\|_{C^{k, \gamma}(\mathbb{R}^{N})}=\sum_{|\beta|\leq k}\|D^{\beta}u\|_{L^{\infty}(\mathbb{R}^{N})}+\sum_{|\beta|=k}\sup_{x\neq y} \frac{|D^{\beta}u(x)-D^{\beta}u(y)|}{|x-y|^{\gamma}}<\infty.
$$

\section{Proof of Theorem \ref{thm1}}
In this section, we investigate the existence of a weak solution to \eqref{P} by means of two abstract results found in \cite[Chapter 4]{Chab} (see also \cite{B}). For the reader's convenience, we state them below. We first introduce some definitions.

Let $X$ be a reflexive Banach space equipped with norm $\|\cdot\|$ and denote by $X^{*}$ its dual.
A map $A:X\rightarrow X^{*}$ is said to be a potential operator with a potential $a:X\rightarrow \mathbb{R}$, if $a$ is Gateaux differentiable and
$$
\lim_{t\rightarrow 0} \frac{a(u+tv)-a(u)}{t}=\langle A(u), v\rangle \quad \mbox{ for all } u, v\in X.
$$
For a potential, we alway assume that $a(0)=0$.
Let $A:X\rightarrow X^{*}$ and $B:X\rightarrow X^{*}$ two potential operators with potentials $a:X\rightarrow \mathbb{R}$ and $b:X\rightarrow \mathbb{R}$, respectively.
Suppose that $b$ is defined on a set $D(b)\subset X$.
For all $\sigma>0$, let $T_{\sigma}: X\rightarrow X$ be a linear map such that $T_{\sigma_1}\circ T_{\sigma_{2}}=T_{\sigma_{1}\sigma_{2}}$ for all $\sigma_{1}, \sigma_{2}>0$, and $T_{1}=id$. Let $E$ be a subspace of $X$ such that $T_{\sigma}(E)\subset E$ for each $\sigma>0$. We assume the following covariance conditions on the potentials $a$ and $b$:
\begin{compactenum}[$(C1)$]
\item $a: X\rightarrow \mathbb{R}$ and there exists $q_{1}\in \mathbb{R}$ such that $a(T_{\sigma}u)=\sigma^{q_{1}}a(u)$ for each $\sigma>0$ and $u\in D(b)$.
\item $b: D(b)\rightarrow \mathbb{R}$ with $E\subset D(b)$, $D(b)+E\subset D(b)$, $T_{\sigma}(D(b))\subset D(b)$, and there exists $q_{2}\in \mathbb{R}$ such that $b(T_{\sigma}u)=\sigma^{q_{2}}b(u)$ for each $\sigma>0$ and $u\in D(b)$.
\end{compactenum}
Assuming $\{u\in D(b): b(u)=1\}\neq \emptyset$, we consider
\begin{align}\label{J1}
J=\inf\{a(u): u\in D(b), \, b(u)=1 \}.
\end{align}
\begin{thm}\label{THM4.1.1}\cite[Theorem 4.1.1-(i)]{Chab}
Suppose that $A$ is a potential operator with a potential $a$. Moreover, suppose that the problem \eqref{J1} has a solution $u$ and that $b$ has a linear continuous Gateaux derivative $\langle B(u), v\rangle$ at all directions $v\in E$. Then
$$
\langle A(u), v\rangle=J \frac{q_{1}}{q_{2}}\langle B(u), v\rangle \quad \mbox{ for all } v\in E.
$$
The functionals $a$ and $b$ have linear continuous Gateaux derivatives at all directions $T_{\sigma}u$, $\sigma>0$, satisfying
$$
\langle A(T_{\sigma}u), v\rangle=\frac{q_{1}}{q_{2}}J \sigma^{q_{1}-q_{2}}\langle B(T_{\sigma}u), v\rangle \quad \mbox{ for all } v\in E \mbox{ and } \sigma>0.
$$
Furthermore, if $\frac{q_{1}}{q_{2}}J>0$ and $q_{1}-q_{2}\neq 0$, then the scaled minimizer $\bar{u}=T_{\bar{\sigma}}u$, with $\bar{\sigma}=\left(\frac{q_{1}}{q_{2}}J\right)^{\frac{1}{q_{2}-q_{1}}}$, satisfies
$$
\langle A(\bar{u}), v\rangle=\langle B(\bar{u}), v\rangle \quad \mbox{ for all } v\in E.
$$
\end{thm}
Next we assume that $a$ and $b$ satisfy $(C1)$ and $(C2)$ with $D(b)=X$ and $q_{1}\neq q_{2}$.
\begin{thm}\label{THM4.1.2}\cite[Theorem 4.1.2]{Chab}
Suppose that there exists a Banach subspace $X_{1}\subset X$ equipped with a norm from $X$ and a mapping $T:X\rightarrow X_{1}$ such that $a(Tu)\leq a(u)$ and $b(Tu)=b(u)$ for each $u\in X$ and let $b(u_{0})>0$ for some $u_{0}\in X$. Moreover, we assume that $a$ is weakly lower semicontinuous and that $b$ admits a decomposition $b(u)=b_{1}(u)-b_{2}(u)$ on $X$, with $b_{1}$ and $b_{2}$ nonnegative on $X$, where $b_{1}$ restricted on $X_{1}$ is weakly continuous, $b_{2}$ is weakly lower semicontinuous and
\begin{compactenum}[$(i)$]
\item there exist constants $\e\in (0, 1)$, $C>0$ and $\alpha>0$ such that
$$
b_{1}(u)\leq \e b_{2}(u)+Ca(u)^{\alpha} \quad \mbox{ for all } u\in X,
$$
\item there exist constants $\beta>0$, $C'>0$ and $\gamma>0$ such that
$$
\|u\|^{\beta}\leq C'(b_{2}(u)+a(u)^{\gamma}) \quad \mbox{ for all } u\in X.
$$
\end{compactenum}
Then problem \eqref{J1} has a nontrivial solution.
\end{thm}

\noindent
Now we are able to prove an existence result for \eqref{P}.
\begin{thm}\label{THMBL}
Let $N\geq 2$, $s\in (0, 1)$ and $p\in (1, \frac{N}{s})$.
Then, \eqref{P} admits a nontrivial weak solution.
\end{thm}
\begin{proof}
Set $X=W^{s,p}(\mathbb{R}^{N})$ and $X_{1}=W^{s, p}_{rad}(\mathbb{R}^{N})=\{u\in W^{s, p}(\mathbb{R}^{N}): u \mbox{ is radially symmetric } \}$. It is easy to check that $X$ is a reflexive Banach space. In fact, let $Y=L^{p}(\mathbb{R}^{N})\times L^{p}(\mathbb{R}^{2N})$ and  introduce the linear isometry $\mathcal{L} :W^{s, p}(\mathbb{R}^{N})\rightarrow Y$ defined as
$$
\mathcal{L}(u)=\left(u, \frac{u(x)-u(y)}{|x-y|^{\frac{N+sp}{p}}} \right).
$$
Clearly, the product space $Y$ is reflexive.
Since $W^{s, p}(\mathbb{R}^{N})$ is a Banach space, $\mathcal{L}(W^{s, p}(\mathbb{R}^{N}))$ is a closed subspace of $Y$. It follows that $\mathcal{L}(W^{s, p}(\mathbb{R}^{N}))$ is reflexive. As a result, $W^{s, p}(\mathbb{R}^{N})$ is a reflexive space for all $s\in (0, 1)$ and $p\in (1, \infty)$.

Set
$$
\|u\|=\left(\frac{C_{N, s, p}}{2p} \iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}}dxdy+\|u\|^{p}_{L^{p}(\mathbb{R}^{N})}\right)^{\frac{1}{p}} \quad \mbox{ for all } u\in X.
$$
Let $T_{\sigma}: X\rightarrow X$ defined as $T_{\sigma}u(x)=u(\frac{x}{\sigma})$, $\sigma>0$, and $T:X\rightarrow X_{1}$ given by $Tu(x)=u^{*}(|x|)$ (Schwarz symmetrization).
Modifying $g$ as in \cite{BL}, we may assume that $g$ satisfies the stronger condition
\begin{align}\label{G2'}
\lim_{|t|\rightarrow \infty} \frac{|g(t)|}{|t|^{p^{*}_{s}-1}}=0,
\end{align}
instead of $(g2)$.
Put
$$
a(u)=\frac{C_{N, s, p}}{2p}\iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}}dxdy \,\,\mbox{ and } \,\, b(u)=\int_{\mathbb{R}^{N}} G(u)dx.
$$
Note that $a, b\in C^{1}(X, \mathbb{R})$, $a(T_{\sigma}u)=\sigma^{N-sp} a(u)$ and $b(T_{\sigma}u)=\sigma^{N}b(u)$ for each $u\in X$ and $\sigma>0$. We also have
$a(Tu)\leq a(u)$ (see \cite[Theorem A.1]{FS}) and $b(Tu)=b(u)$ (see \cite[Appendix]{BL}) for all $u\in X$. Clearly, $a$ is weakly lower semicontinuous on $X$.

Let us now show that the set
$\{u\in X: b(u)=1\}$ is nonempty.
For $R>0$, define
\begin{equation*}
w_{R}(x)=\left\{
\begin{array}{ll}
\zeta & \mbox{ for } |x|\leq R, \\
\zeta(R+1-|x|) & \mbox{ for } R\leq |x|\leq R+1, \\
0 & \mbox{ for } |x|\geq R+1,
\end{array}
\right.
\end{equation*}
where $\zeta>0$ is given in $(g3)$.
It is clear that $w_{R}\in X$ and
$$
b(w_{R})\geq G(\zeta)|B_{R}(0)|-|B_{R+1}(0)\setminus B_{R}(0)| \max_{t\in [0, \zeta]} |G(t)|.
$$
Then there exist two constants $C_{1}, C_{2}>0$ such that
$$
b(w_{R})\geq C_{1}R^{N}-C_{2}R^{N-1},
$$
and so $b(w_{R})>0$ for $R>0$ large enough.

Next we consider the following constrained minimization problem:
\begin{align}\label{J}
J=\inf\{a(u): u\in X, \, b(u)=1\}.
\end{align}
Let $g_{1}(t)=(g(t)+mt^{p-1})^{+}$ and $g_{2}(t)=g_{1}(t)-g(t)$ for all $t\geq 0$. We extend $g_{1}$ and $g_{2}$ as odd functions for $t\leq 0$. Set $b_{1}(u)=\int_{\mathbb{R}^{N}} G_{1}(u)dx$ and $b_{2}(u)=\int_{\mathbb{R}^{N}}G_{2}(u)dx$, where $G_{i}(t)=\int_{0}^{t}g_{i}(\tau)d\tau$ for $i=1, 2$. Then, $b(u)=b_{1}(u)-b_{2}(u)$ for all $u\in X$. Clearly, $b_{1}(u), b_{2}(u)\geq 0$ for all $u\in X$, and $b_{2}$ is weakly lower semicontinuous on $X$. Let us prove that $b_{1}$ is weakly continuous on $X_{1}$. Let $(u_{n})\subset X_{1}$ be such that $u_{n}\rightharpoonup u$ in $X_{1}$ and we show that $b_{1}(u_{n})\rightarrow b_{1}(u)$. By \cite[Theorem II.1]{Lions}, we know that $X_{1}$ is compactly embedded into $L^{r}(\mathbb{R}^{N})$ for all $r\in (p, p^{*}_{s})$. We note that $G_{1}\in C(\mathbb{R})$, $\lim_{|t|\rightarrow 0} \frac{G_{1}(t)}{|t|^{p}}=0=\lim_{|t|\rightarrow \infty} \frac{G_{1}(t)}{|t|^{p^{*}_{s}}}$ (by definition of $g_{1}$ and assumptions $(g1)$ and \eqref{G2'}), $\sup_{n\in \mathbb{N}}\|u_{n}\|\leq C$ and $G(u_{n})\rightarrow G(u)$ a.e in $\mathbb{R}^{N}$. Applying \cite[Lemma 2.4]{CW}, we deduce that $\|G(u_{n})-G(u)\|_{L^{1}(\mathbb{R}^{N})}\rightarrow 0$ and thus $b_{1}(u_{n})\rightarrow b_{1}(u)$ as desired.

Hereafter, we verify that assumptions $(i)$ and $(ii)$ of Theorem \ref{THM4.1.2} hold.
Using $(g1)$ and $(g2)$, we see that for all $\e>0$ there exists a constant $C_{\e}>0$ such that
$$
G_{1}(t)\leq \e G_{2}(t)+C_{\e}|t|^{p^{*}_{s}} \quad \mbox{ for all } t\in \mathbb{R},
$$
and thanks to the Sobolev inequality (see \cite[Theorem 6.5]{DPV}) we arrive at
$$
b_{1}(u)\leq \e b_{2}(u)+C_{\e} C_{*} a(u)^{\frac{p^{*}_{s}}{p}} \quad \mbox{ for all } u\in X,
$$
for some constant $C_{*}>0$.
On the other hand, because $g_{2}(t)\geq m t^{p-1}$ for all $t\geq 0$, we have
$$
b_{2}(u)\geq \frac{m}{p}\|u\|^{p}_{L^{p}(\mathbb{R}^{N})} \quad \mbox{ for all } u\in X,
$$
whence
$$
\|u\|^{p}\leq \frac{p}{m} b_{2}(u)+a(u) \quad \mbox{ for all } u\in X.
$$
Then, invoking Theorem \ref{THM4.1.2}, we infer that problem \eqref{J} has a nontrivial solution $u$. Observing that $\frac{(N-sp)}{N}J>0$ and $(N-sp)-N=-sp\neq 0$, we can invoke Theorem \ref{THM4.1.1} to deduce that $\bar{u}=T_{\bar{\sigma}}u$, where $\bar{\sigma}=(\frac{N-sp}{N}J)^{\frac{1}{sp}}$, satisfies
\begin{align*}
\langle a'(\bar{u}), v\rangle=\langle b'(\bar{u}), v\rangle \quad \mbox{ for all } v\in X,
\end{align*}
that is $\bar{u}$ is a weak solution to \eqref{P}.
\end{proof}

As byproduct of Theorem \ref{THMBL}, we obtain Theorem \ref{thm1}.
\begin{proof}[Proof of Theorem \ref{thm1}]
Let $\bar{u}$ be the solution found in Theorem \ref{THMBL}.
From $a(\bar{u})=\bar{\sigma}^{N-sp}J$, $J>0$, $b(\bar{u})=\bar{\sigma}^{N}$ and $\bar{\sigma}=(\frac{N-sp}{N}J)^{\frac{1}{sp}}$, we can see that
\begin{align*}
a(\bar{u})=\frac{N}{N-sp}b(\bar{u}),
\end{align*}
that is, \eqref{POHOZAEV} holds.
\end{proof}

\begin{remark}
The proof of Theorem \ref{thm1} still works in the case $s=1$.
\end{remark}

\begin{remark}
Reasoning as in \cite[Lemma 3.18]{AIDCDS}, we can show that if $u$ is a nonnegative weak solution to \eqref{P}, then $u\in C^{0, \alpha}(\mathbb{R}^{N})\cap L^{r}(\mathbb{R}^{N})$ for all $r\in [p, \infty]$ and for some $\alpha\in (0, 1)$. Since $(-\Delta)^{s}_{p}u+m u^{p-1}\leq g_{1}(u)$ in $\mathbb{R}^{N}$, where $g_{1}(t)=o(t^{p-1})$ as $t\rightarrow 0$, we can proceed as in  \cite[Theorem 1.5]{DPQ} to see that $|u(x)|\leq C|x|^{-\frac{N+sp}{p-1}}$ for all $|x|$ large enough. If in addition $u\in C^{1, 1}_{loc}(\mathbb{R}^{N})$ and $g\in C^{1}(\mathbb{R})$, then $u$ must be radially symmetric and monotone decreasing about some point in $\mathbb{R}^N$ in view of  \cite[Theorem 5]{CL}.
\end{remark}

\section{Proof of Theorem \ref{thm2}}

We start by giving a useful integral representation formula for $(-\Delta)^{s}_{p}$ along functions of class $C^{1, 1}(\mathbb{R}^{N})$.
This formula is motivated by \cite[Lemma 3.2]{DPV} in which the authors dealt with $p=2$ and functions in the Schwartz space of rapidly decaying functions.
\begin{lem}\label{RIF}
Let $N\in \mathbb{N}$, $s\in (0, 1)$ and $p\in (1, \infty)$. When $p\in (1, 2)$, we assume that $s<\frac{2(p-1)}{p}$. Then, for all $u\in C^{1, 1}(\mathbb{R}^{N})$, it holds, for all $x\in \mathbb{R}^N$,
\begin{align*}
(-\Delta)^{s}_{p}u(x)=\frac{C_{N, s, p}}{2} \int_{\mathbb{R}^{N}} \frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))+|u(x)-u(x-z)|^{p-2}(u(x)-u(x-z))}{|z|^{N+sp}} dz.
\end{align*}
\end{lem}
\begin{proof}
Pick $u\in C^{1, 1}(\mathbb{R}^{N})$ and fix $x\in \mathbb{R}^{N}$. Using the change of variables theorem and the symmetry of the kernel $\frac{1}{|z|^{N+sp}}$, we have
\begin{align*}
&P.V.\int_{\mathbb{R}^{N}} \frac{|u(x)-u(y)|^{p-2}(u(x)-u(y))}{|x-y|^{N+sp}} dy\\
&\quad =P.V.\int_{\mathbb{R}^{N}} \frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))}{|z|^{N+sp}} dz \\
&\quad =\frac{1}{2} P.V.\int_{\mathbb{R}^{N}} \frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))+|u(x)-u(x-z)|^{p-2}(u(x)-u(x-z))}{|z|^{N+sp}} dz.
\end{align*}
Let us now show that the $P.V.$ in the above formula can be removed.
Fix $\e\in (0, 1)$. Then we can write
\begin{align*}
&\int_{\mathbb{R}^{N}\setminus B_{\e}(0)} \frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))+|u(x)-u(x-z)|^{p-2}(u(x)-u(x-z))}{|z|^{N+sp}} dz\\
&\quad =\int_{\mathbb{R}^{N}\setminus B_{1}(0)} \frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))+|u(x)-u(x-z)|^{p-2}(u(x)-u(x-z))}{|z|^{N+sp}} dz\\
&\quad \quad +\int_{B_{1}(0)\setminus B_{\e}(0)} \frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))+|u(x)-u(x-z)|^{p-2}(u(x)-u(x-z))}{|z|^{N+sp}} dz.
\end{align*}
We first prove that
$$
\int_{\mathbb{R}^{N}\setminus B_{1}(0)} \frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))+|u(x)-u(x-z)|^{p-2}(u(x)-u(x-z))}{|z|^{N+sp}}dz<\infty.
$$
Take $|z|\geq 1$. Since $u\in L^{\infty}(\mathbb{R}^{N})$, we see that
\begin{align*}
\left|\frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))+|u(x)-u(x-z)|^{p-2}(u(x)-u(x-z))}{|z|^{N+sp}}\right|\leq \frac{2^{p}\|u\|^{p-1}_{L^{\infty}(\mathbb{R}^{N})}}{|z|^{N+sp}}.
\end{align*}
Because $\frac{1}{|z|^{N+sp}}$ is integrable at infinity, the assertion is proved.
Now we claim that
$$
\int_{B_{1}(0)\setminus B_{\e}(0)} \frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))+|u(x)-u(x-z)|^{p-2}(u(x)-u(x-z))}{|z|^{N+sp}}dz<\infty.
$$
Let $\e\leq |z|<1$. We first assume $p\geq 2$.
Recalling that (see \cite[pag. 255]{Simon})
$$
||a|^{p-2}a-|b|^{p-2}b|\leq C_{p}(|a|+|b|)^{p-2}|a-b| \quad \mbox{ for all } a, b\in \mathbb{R},
$$
and using the mean value theorem, we obtain
\begin{align*}
&\left|\frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))+|u(x)-u(x-z)|^{p-2}(u(x)-u(x-z))}{|z|^{N+sp}}\right|\\
&\quad =\left|\frac{|\nabla u(x+\theta_{1}z)\cdot z|^{p-2}\nabla u(x+\theta_{1}z)\cdot z-|\nabla u(x-\theta_{2}z)\cdot z|^{p-2}\nabla u(x-\theta_{2}z)\cdot z}{|z|^{N+sp}}\right|\\
&\quad \leq 2^{p-2}C_{p} \|\nabla u\|^{p-2}_{L^{\infty}(\mathbb{R}^{N})}|z|^{p-2} \frac{|[\nabla u(x+\theta_{1}z)-\nabla u(x-\theta_{2}z)]\cdot z|}{|z|^{N+sp}} \\
&\quad \leq 2^{p-2}C_{p} \|\nabla u\|^{p-2}_{L^{\infty}(\mathbb{R}^{N})} |z|^{p-2} \frac{[\nabla u]_{C^{0, 1}(\mathbb{R}^{N})}|z|^{2}}{|z|^{N+sp}} \\
&\quad \leq 2^{p-2} C_{p} \|u\|_{C^{1, 1}(\mathbb{R}^{N})}^{p-1}\frac{1}{|z|^{N+sp-p}},
\end{align*}
where $\theta_{1}, \theta_{2}\in (0, 1)$. Since
$$
\frac{1}{|z|^{N+sp-p}}
$$
is integrable near the origin (thanks to $s\in (0, 1)$), we reach the assertion. Now we suppose $p\in (1, 2)$. Then, observing that (see \cite[pag. 255]{Simon})
$$
||a|^{p-2}a-|b|^{p-2}b|\leq C'_{p}|a-b|^{p-1} \quad \mbox{ for all } a, b\in \mathbb{R},
$$
we get
\begin{align*}
&\left|\frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))+|u(x)-u(x-z)|^{p-2}(u(x)-u(x-z))}{|z|^{N+sp}}\right|\\
&\quad \leq \frac{|\nabla u(x+\theta_{1}z)\cdot z-\nabla u(x-\theta_{2}z)\cdot z|^{p-1}}{|z|^{N+sp}}\\
&\quad \leq C'_{p} \|\nabla u\|^{p-1}_{C^{0, 1}(\mathbb{R}^{N})} |z|^{p-1} \frac{|z|^{p-1}}{|z|^{N+sp}} \\
&\quad \leq C'_{p} \|u\|_{C^{1, 1}(\mathbb{R}^{N})}^{p-1}\frac{1}{|z|^{N+sp-2p+2}},
\end{align*}
where $\theta_{1}, \theta_{2}\in (0, 1)$. Due to the fact that
$$
\frac{1}{|z|^{N+sp-2p+2}}
$$
is integrable near the origin (because of $s<\frac{2(p-1)}{p}$), we can conclude the proof.
\end{proof}

\begin{remark}
From the proof of Lemma \ref{RIF}, we derive that $(-\Delta)^{s}_{p}u\in L^{\infty}(\mathbb{R}^{N})\cap C^{0}(\mathbb{R}^{N})$ for all $u\in C^{1, 1}(\mathbb{R}^{N})$, with $s<\frac{2(p-1)}{p}$ if $p\in (1, 2)$.
\end{remark}

Next we prove a useful integration by parts formula for $W^{s, p}(\mathbb{R}^{N})\cap C^{1, 1}(\mathbb{R}^{N})$ functions and vector fields of class $C^{0, 1}_{c}(\mathbb{R}^N, \mathbb{R}^N)$.
\begin{lem}\label{LEMMA}
Let $N\in \mathbb{N}$, $s\in (0, 1)$ and $p\in (1, \infty)$. When $p\in (1, 2)$, we assume that $s<\frac{2(p-1)}{p}$. Let $u\in W^{s, p}(\mathbb{R}^{N})\cap C^{1, 1}(\mathbb{R}^{N})$ and $\mathcal{X}\in C^{0, 1}_{c}(\mathbb{R}^{N}, \mathbb{R}^{N})$.
Then we have
\begin{align*}
&\frac{C_{N, s, p}}{2} \iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}} \left[\dive(\mathcal{X}(x))+\dive(\mathcal{X}(y))-(N+sp) \frac{(\mathcal{X}(x)-\mathcal{X}(y))\cdot (x-y)}{|x-y|^{2}} \right]  dxdy\\
&=-p \int_{\mathbb{R}^{N}} \mathcal{X}(x)\cdot \nabla u(x) \, (-\Delta)^{s}_{p}u(x) \,dx.
\end{align*}
\end{lem}
\begin{proof}
We modify in a suitable way some arguments found in \cite[Lemma 4.2]{DFW1}.
Exploiting the symmetry of the kernel $\frac{1}{|x|^{N+sp}}$ and Fubini's theorem, we see that
\begin{align*}
&\frac{C_{N, s, p}}{2} \iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}} \left[\dive(\mathcal{X}(x))+\dive(\mathcal{X}(y))-(N+sp) \frac{(\mathcal{X}(x)-\mathcal{X}(y))\cdot (x-y)}{|x-y|^{2}} \right]  dxdy \\
&\quad =-\frac{C_{N, s, p}(N+sp)}{2} \iint_{\mathbb{R}^{2N}} |u(x)-u(y)|^{p} \frac{(x-y)\cdot (\mathcal{X}(x)-\mathcal{X}(y))}{|x-y|^{N+sp+2}} dxdy\\
&\quad \quad+\frac{C_{N, s, p}}{2} \iint_{\mathbb{R}^{2N}}  |u(x)-u(y)|^{p} \frac{\dive \mathcal{X}(x)+\dive \mathcal{X}(y)}{|x-y|^{N+sp}}dxdy\\
&\quad =-\frac{C_{N, s, p}(N+sp)}{2} \lim_{\mu\rightarrow 0}\iint_{|x-y|>\mu} |u(x)-u(y)|^{p} \frac{(x-y)\cdot (\mathcal{X}(x)-\mathcal{X}(y))}{|x-y|^{N+sp+2}} dxdy\\
&\quad \quad+C_{N, s, p} \iint_{\mathbb{R}^{2N}}  |u(x)-u(y)|^{p} \frac{\dive \mathcal{X}(x)}{|x-y|^{N+sp}}dxdy\\
&\quad =-C_{N, s, p}(N+sp) \lim_{\mu\rightarrow 0} \int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}\setminus \overline{B_{\mu}(y)}} |u(x)-u(y)|^{p}  \frac{(x-y)\cdot \mathcal{X}(x)}{|x-y|^{N+sp+2}}dxdy\\
&\quad \quad+C_{N, s, p} \iint_{\mathbb{R}^{2N}}  |u(x)-u(y)|^{p} \frac{\dive \mathcal{X}(x)}{|x-y|^{N+sp}}dxdy.
\end{align*}
Now, applying the divergence theorem in the domain $\mathbb{R}^{N}\setminus \overline{B_{\mu}(y)}$, with $y\in \mathbb{R}^{N}$ and $\mu>0$ fixed, and using the fact that
$$
\nabla_{x}(|x-y|^{-(N+sp)})=-(N+sp)\frac{x-y}{|x-y|^{N+sp+2}},
$$
we have
\begin{align}\label{IBPF1}
&\frac{C_{N, s, p}}{2} \iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}} \left[\dive(\mathcal{X}(x))+\dive(\mathcal{X}(y))-(N+sp) \frac{(\mathcal{X}(x)-\mathcal{X}(y))\cdot (x-y)}{|x-y|^{2}} \right]  dxdy \nonumber\\
&\quad=C_{N, s, p} \lim_{\mu\rightarrow 0} \int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}\setminus \overline{B_{\mu}(y)}} |u(x)-u(y)|^{p} \, \nabla_{x}(|x-y|^{-(N+sp)}) \cdot \mathcal{X}(x) \, dxdy \nonumber\\
&\quad\quad+C_{N, s, p} \iint_{\mathbb{R}^{2N}}  |u(x)-u(y)|^{p} \frac{\dive \mathcal{X}(x)}{|x-y|^{N+sp}}dxdy \nonumber\\
&\quad=-C_{N, s, p} \lim_{\mu\rightarrow 0}\int_{\mathbb{R}^{N}} \int_{\mathbb{R}^{N}\setminus \overline{B_{\mu}(y)}} |u(x)-u(y)|^{p} \frac{\dive \mathcal{X}(x)}{|x-y|^{N+sp}}dxdy \nonumber\\
&\quad\quad -p C_{N, s, p} \lim_{\mu\rightarrow 0}\int_{\mathbb{R}^{N}} \int_{\mathbb{R}^{N}\setminus \overline{B_{\mu}(y)}} |u(x)-u(y)|^{p-2} (u(x)-u(y))  \frac{\nabla u(x)\cdot \mathcal{X}(x)}{|x-y|^{N+sp}}dxdy \nonumber\\
&\quad\quad +C_{N, s, p}\lim_{\mu\rightarrow 0}  \int_{\mathbb{R}^{N}} \int_{\partial B_{\mu}(y)} |u(x)-u(y)|^{p} \frac{(y-x)\cdot \mathcal{X}(x)}{|x-y|^{N+sp+1}} \, d\sigma(y) dx \nonumber\\
&\quad\quad+C_{N, s, p} \iint_{\mathbb{R}^{2N}}  |u(x)-u(y)|^{p} \frac{\dive \mathcal{X}(x)}{|x-y|^{N+sp}}dxdy \nonumber\\
&\quad=-pC_{N, s, p} \lim_{\mu\rightarrow 0} \iint_{|x-y|>\mu} |u(x)-u(y)|^{p-2}(u(x)-u(y)) \frac{\nabla u(x)\cdot \mathcal{X}(x)}{|x-y|^{N+sp}} dxdy\nonumber\\
&\quad\quad+C_{N, s, p} \lim_{\mu\rightarrow 0}  \mu^{-N-1-sp} \iint_{|x-y|=\mu} |u(x)-u(y)|^{p} (y-x)\cdot \mathcal{X}(x) d\sigma(x, y) \nonumber\\
&\quad=-\frac{pC_{N, s, p}}{2} \lim_{\mu\rightarrow 0} \int_{\mathbb{R}^{N}} \nabla u(x)\cdot \mathcal{X}(x) \int_{\mathbb{R}^{N}\setminus \overline{B_{\mu}(y)}} \left( \frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))}{|z|^{N+sp}} \right. \nonumber\\
&\quad\quad \left.+\frac{|u(x)-u(x-z)|^{p-2}(u(x)-u(x-z))}{|z|^{N+sp}}\right) dzdx\nonumber\\
&\quad\quad+\frac{C_{N, s, p}}{2} \lim_{\mu\rightarrow 0}  \mu^{-N-1-sp} \iint_{|x-y|=\mu} |u(x)-u(y)|^{p} (y-x)\cdot (\mathcal{X}(x)-\mathcal{X}(y)) d\sigma(x, y)\nonumber\\
&\quad=-p Y_{1}+Y_{2}.
\end{align}
Since $u\in C^{1, 1}(\mathbb{R}^{N})$, it follows from Lemma \ref{RIF} that
\begin{align}\label{IBPF2}
Y_{1}&=\frac{C_{N, s, p}}{2}  \int_{\mathbb{R}^{N}}  \nabla u(x)\cdot \mathcal{X}(x) \int_{\mathbb{R}^{N}} \left(\frac{|u(x)-u(x+z)|^{p-2}(u(x)-u(x+z))}{|z|^{N+sp}}\right. \nonumber\\
&\quad \left. +\frac{|u(x)-u(x-z)|^{p-2}(u(x)-u(x-z))}{|z|^{N+sp}} \right) dzdx \nonumber\\
&=\int_{\mathbb{R}^{N}} \nabla u(x)\cdot \mathcal{X}(x)  \, (-\Delta)^{s}_{p}u(x) \,dx.
\end{align}
On the other hand, since $\mathcal{X}$ has compact support, we can find $R>0$ large enough such that $(\mathcal{X}(x)-\mathcal{X}(y))=0$ for all $x, y\in B_{R}(0)$ with $|x-y|<1$. For all $\mu\in (0, 1)$, we put
$$
\mathcal{N}_{\mu}=\Bigl\{(x, y)\in B_{R}(0)\times B_{R}(0): |x-y|=\mu \Bigr\}.
$$
Then, using $u\in C^{0, 1}(\mathbb{R}^{N})$, $\mathcal{X}\in C^{0, 1}(\mathbb{R}^{N}, \mathbb{R}^{N})$, and that the $2N-1$-dimensional measure of the set $\mathcal{N}_{\mu}$ is of order $O(N-1)$ as $\mu\rightarrow 0$, we have
\begin{align}\begin{split}\label{IBPF3}
&\mu^{-N-1-sp} \iint_{|x-y|=\mu} |u(x)-u(y)|^{p}(y-x)\cdot (\mathcal{X}(x)-\mathcal{X}(y)) d\sigma(x, y)\\
&\quad=\mu^{-N-1-sp} \iint_{\mathcal{N}_{\mu}} |u(x)-u(y)|^{p}(y-x)\cdot (\mathcal{X}(x)-\mathcal{X}(y)) d\sigma(x, y)=O(\mu^{p(1-s)})\rightarrow 0 \quad \mbox{ as } \mu\rightarrow 0,
\end{split}\end{align}
that is $Y_{2}=0$. Combining \eqref{IBPF1}, \eqref{IBPF2} and \eqref{IBPF3}, we obtain the desired formula.
\end{proof}
Now, we are ready to give the proof of Theorem \ref{thm2}.
\begin{proof}[Proof of Theorem \ref{thm2}]
Let $u\in W^{s, p}(\mathbb{R}^{N})\cap C^{1, 1}(\mathbb{R}^{N})$ be such that $(-\Delta)^{s}_{p}u=g(u)$ in $\mathbb{R}^{N}$. Pick $\varphi\in C^{1}_{c}(\mathbb{R}^{N})$ such that $0\leq \varphi\leq 1$ in $\mathbb{R}^{N}$, $\varphi(x)=1$ for $|x|\leq 1$ and $\varphi(x)=0$ for $|x|\geq 2$. Put $\varphi_{\lambda}(x)=\varphi(\lambda x)$ for all $x\in \mathbb{R}^{N}$ and $\lambda>0$.
Note that, for all $x\in \mathbb{R}^{N}$ and $\lambda>0$,
\begin{align}\label{WS}
0\leq\varphi_{\lambda}(x)\leq 1 \,\mbox{ and } \, |x||\nabla \varphi_{\lambda}(x)|\leq C_{1},
\end{align}
 where the constant $C_{1}>0$ is independent of $\lambda$.
Since
$\mathcal{X}(x)=\varphi_{\lambda}(x)x\in C^{0, 1}_{c}(\mathbb{R}^{N}, \mathbb{R}^{N})$ (note that \eqref{WS} implies that $|x\varphi_{\lambda}(x)-y\varphi_{\lambda}(y)|\leq C_{2}|x-y|$ for some constant $C_{2}>0$ independent of $\lambda$) and $u\in W^{s, p}(\mathbb{R}^{N})\cap C^{1, 1}(\mathbb{R}^{N})$, we can apply
Lemma \ref{LEMMA} to deduce that
\begin{align*}
&-\frac{C_{N, s, p}}{2p} \iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}} \left[\dive(\mathcal{X}(x))+\dive(\mathcal{X}(y))-(N+sp) \frac{(\mathcal{X}(x)-\mathcal{X}(y))\cdot (x-y)}{|x-y|^{2}} \right]  dxdy\\
&=\int_{\mathbb{R}^{N}} \mathcal{X}(x)\cdot \nabla u(x) \, (-\Delta)^{s}_{p}u(x) \,dx\\
&=\int_{\mathbb{R}^{N}} \mathcal{X}(x)\cdot \nabla u(x) \, g(u(x)) \,dx.
\end{align*}
An integration by parts ensures that
\begin{align*}
-\int_{\mathbb{R}^{N}} g(u(x)) \varphi_{\lambda}(x)x\cdot \nabla u(x)dx=N \int_{\mathbb{R}^{N}} \varphi_{\lambda}(x)G(u(x))dx+\lambda \int_{\mathbb{R}^{N}} G(u(x)) x\cdot \nabla \varphi_{\lambda}(x)dx,
\end{align*}
and so
\begin{align}\label{BRUN1}
&-\frac{C_{N, s, p}}{2p} \iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}} \left[\dive(\mathcal{X}(x))+\dive(\mathcal{X}(y))-(N+sp) \frac{(\mathcal{X}(x)-\mathcal{X}(y))\cdot (x-y)}{|x-y|^{2}} \right]  dxdy \nonumber\\
&=N \int_{\mathbb{R}^{N}} \varphi_{\lambda}(x)G(u(x))dx+\lambda \int_{\mathbb{R}^{N}} G(u(x)) x\cdot \nabla \varphi_{\lambda}(x)dx.
\end{align}
From \eqref{WS} we derive that
$$
\left| \dive(\mathcal{X}(x))+\dive(\mathcal{X}(y))-(N+sp) \frac{(\mathcal{X}(x)-\mathcal{X}(y))\cdot (x-y)}{|x-y|^{2}}\right|\leq C_{3},
$$
for some constant $C_{3}>0$ independent of $\lambda$. In view of this fact, $u\in W^{s, p}(\mathbb{R}^{N})$, \eqref{WS}, the pointwise convergences $\varphi_{\lambda}(x)\rightarrow 1$, $\mathcal{X}(x)\rightarrow x$ and $\dive\mathcal{X}(x)\rightarrow N$ for a.e. $x\in \mathbb{R}^N$, as $\lambda\rightarrow 0$, and $G(u)\in L^{1}(\mathbb{R}^{N})$, we can apply the dominated convergence theorem to infer that, as $\lambda\rightarrow 0$,
\begin{align}\label{BRUN2}
&-\frac{C_{N, s, p}}{2p} \iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}} \left[\dive(\mathcal{X}(x))+\dive(\mathcal{X}(y))-(N+sp) \frac{(\mathcal{X}(x)-\mathcal{X}(y))\cdot (x-y)}{|x-y|^{2}} \right]  dxdy \nonumber\\
&\rightarrow -\frac{C_{N, s, p}(N-sp)}{2p} \iint_{\mathbb{R}^{2N}} \frac{|u(x)-u(y)|^{p}}{|x-y|^{N+sp}} dxdy,
\end{align}
and
\begin{align}\label{BRUN3}
N \int_{\mathbb{R}^{N}} \varphi_{\lambda}(x)G(u(x))dx+\lambda \int_{\mathbb{R}^{N}} G(u(x)) x\cdot \nabla \varphi_{\lambda}(x)dx\rightarrow N \int_{\mathbb{R}^{N}} G(u(x))dx.
\end{align}
On account of \eqref{BRUN1}, \eqref{BRUN2} and \eqref{BRUN3}, we arrive at \eqref{POHOZAEV}.
\end{proof}


\section{Introduction} \label{s:intro}

Magnetism is found to be present in about 10 to 20\% of white dwarfs \citep{Kawka2007, Giammichele2012, Bagnulo2021}.
According to the compilation of \citet{Ferrario2015a} and from the Montreal White Dwarf Database \citep[MWDD,][]{MWDD}, the vast majority of magnetic white dwarfs (MWDs) identified in the literature are found to have a hydrogen-rich atmosphere.
Such white dwarfs have been examined in detail in the first paper of this series \citep[][hereafter Paper I]{Hardy2022}.
On the other hand, magnetic white dwarfs with elements other than hydrogen as their main constituent appear to be much rarer.
Magnetism has been detected in a handful of cool helium-rich objects such as DQs \citep{Schmidt1995, Schmidt1999, Vornanen2010, Vornanen2013} and DZs \citep{Reid2001, Dufour2006, Hollands2015}.
Magnetism is also detected through circular spectropolarimetry in a number of featureless cool DC white dwarfs, some of which are also probably helium-dominated \citep[][]{Putney1997}.
On the hotter side, a large fraction of the carbon-oxygen dominated atmosphere white dwarfs, the so-called "Hot DQ", also exhibit traces of magnetism \citep{Dufour2008, Dufour2010}.
In the DB temperature range (roughly 12,000 K to 30,000 K), however, the incidence of magnetism appear to be quite rare.
While a very weak magnetic field (1.3 kG) has been detected through asteroseismology in GD358 \citep[][]{Winget1994}, unambiguous detection of helium in a magnetic white dwarf was only first confirmed with detailed modeling by \citet{Achilleos1992} for the star Feige 7 ($\sim$ 35 MG).
A weaker magnetic field ($\sim$ 100 kG) has also been detected from circular polarization in the DBA white dwarf LB 8827 \citep[][]{Wesemael2001}. More recently, \citet{Richer2019} announced the discovery of a magnetic DB star in the M47 cluster.
This star shows magnetic splitting of helium lines coherent with a field between 2 and 3 MG, but also features that indicate it may be a binary system (MWD + late-type companion). Four magnetic DBs have also been identified in the course of the Hamburgh ESO survey by \citet{Reimers1998}, although no formal analysis has ever been attempted to our knowledge.
Finally, from comparison of observed features with Zeeman splitting calculations of helium lines, three more magnetic DB, GD 229, HE 1211$-$1707, and HE 1043$-$0502 (they are included in our sample under the names J2012+3113, J1214$-$1724 and J1046$-$0518 respectively), have been identified \citep[][]{Jordan1998, Wickramasinghe2002}.
Detailed modeling similar to that presented for magnetic DAs in Paper I, however, has never been attempted for all these objects.

In Paper I, we examined 651 objects that had previously been classified as magnetic hydrogen-rich white dwarfs in the literature.
During this process, we found 38 objects that we could not fit within our assumed theoretical framework, suggesting that some of these could have a non-hydrogen-rich composition.
Furthermore, the MWDD lists numerous white dwarfs classified as DBH that, to our knowledge, have never been analyzed properly using magnetic model atmosphere.
We thus aim here to visit/revisit all those potentially helium-rich objects in light of the mass availability of precise trigonometric parallax measurements \citep[][]{Gaia2016,Gaia2018,Gaia2021} and broadband photometry \citep[][]{York2000, Chambers2016}.
However, since our sample consists mostly of spectra with low signal-to-noise ratio and medium resolution from the Sloan Digital Sky Survey (see below), we restrict ourselves to objects in the DB temperature range ($\sim$ 12,000-30,000 K) that should show detectable signs of helium-line splitting when highly magnetic ($B \ge 2 MG $).
This paper aims to model those objects with the assumption of a helium-rich composition and offset dipole geometry, leaving other compositions and geometries to future studies.
In \autoref{s:observations}, we present our sample, in \autoref{s:theory} we describe our theoretical framework and fitting method.
\autoref{s:results} present our detailed analysis and assessment for each star and we conclude in \autoref{s:conclu}.

\section{Observations} \label{s:observations}

\subsection{Sample Selection} \label{ss:sample}

We aim to analyze in a homogeneous fashion all known DB magnetic white dwarfs that have been identified in the literature over the years \citep[MWDD,][and references therein]{Ferrario2015a, Ferrario2020}.
We also include the many objects that have been flagged as DBH by visual inspection of thousands of SDSS spectra \citep[32 stars, mostly from][]{Kleinman2013,Kepler2015, Kepler2016} as well as all the objects that we were unable to fit using hydrogen-rich models in Paper I (38 stars).
This brings the total number of potentially magnetic DB white dwarfs in our sample to 83 stars when removing duplicates found in more than one of the above categories.
However, optical spectra were not available for 4 objects, leaving us with a total sample of 79 stars to examine (see Table 1).
Apart from those mentioned in the introduction, these white dwarfs have never, to our knowledge, been subject to a detailed analysis using magnetic synthetic spectra.

\input{names}

In order to have an analysis as homogeneous as possible, we use in priority SDSS \citep{York2000} \emph{ugriz} point-spread function magnitudes when available. If \emph{ugriz} data was not available for an object, we relied on Pan-STARRS \citep{Chambers2016} \emph{grizy} or previously published \emph{JHK} and \emph{BVRI} photometry, in that order of priority.
Of our 79 objects, 13 ($\sim 16$\%) had no \emph{ugriz} data.
For those situations, we had to use \emph{grizy} data for 10 stars, \emph{BVRI} photometry for 2 stars and \emph{B + JHK} for a single star.

The majority (67 out of 79; 84\%) of the spectra used for this work are collected from SDSS.
The spectra of the 12 objects without SDSS spectroscopy were taken from
\citet{Kilic2020, Bergeron1997, Bergeron2001, Giammichele2012, Limoges2015}
or archival data obtained by the Montreal group over the years.
Spectra for the stars presented in \citet{Wickramasinghe2002} (GD 229, HE 1211$-$1707, and HE 1043$-$0502) were kindly made available to us by L. Ferrario (private communication).

\section{Theoretical Framework} \label{s:theory}

The theoretical framework and fitting method used in this study are very similar to that described in Paper I, except that a pure helium composition is assumed. Below we briefly review the methodology, emphasizing the differences with Paper I.

The first detailed investigation of magnetic white dwarfs with helium absorption features, that of \citet{Achilleos1992}, relied on line data calculation from \citet{Kemic1974}.
Those first order calculations, however, are only valid up to about 20 MG and tables quickly need to be extrapolated outside their range of validity for most practical purposes.
Subsequently, more detailed calculation performed in Heidelberg provided precise data for a large number of energy states in the 300-700 MG range that allowed the identification of stationary neutral helium lines \citep[][]{Jordan1998}, but the full parameter space has yet to be explored systematically using those data.

In this work we use the helium line splitting data described in \citet{Jordan2001}, \citet{Becken1999} and \citet{Becken2000}.
The data contained transitional energies, wavelengths, and various quantum quantities needed to compute oscillator strengths, for a wide range of magnetic field intensities from 0 to over 200 GG ($200 \times 10^3$ MG), for each of those components.
We implemented an interpolating scheme developed by \citet{Steffen1990} to obtain the component data for arbitrary field intensities, with a linear extrapolation when required.
This extrapolation occurs at over 200 GG for most components, although some weaker components only have calculations up to 500 MG.
As most stars in our sample have a dipolar magnetic field $B_p$ of under 100 MG with only a very limited number of stars with $B_p$ over 500 MG, this extrapolation is not expected to be frequent.
A wide subset of the line data is illustrated in \autoref{f:helium}, where the line position in relation to the applied magnetic field can be seen.
For magnetic field strength above about 5 MG, large non-linearities become very important with many components shifting rapidly by hundreds of Angstrom as the field strength increases.
As a result, it becomes difficult to recognize the helium signatures in the most magnetic objects.
Note that for this exploratory work, we simply assume that the broadening treatment for each line component is the same as for non-magnetic white dwarfs as no broadening calculations in the presence of strong magnetic field have been performed to our knowledge.
Strongly magnetic white dwarf spectra may thus not be accurately be represented within our theoretical framework.
\autoref{f:repr_synth} shows typical helium-rich synthetic spectra assuming a dipolar geometry for various field strength, inclination and dipole offset.
For some geometries at the highest field intensities, the spread of magnetic field strength over the surface leads to absorption features spread all over large wavelength range, making lines almost unrecognizable.
In some cases, the density of lines can almost mimic a featureless DC spectrum while it could easily be interpreted as noise in others.
As we will see below, some of these configurations may explain hot stars with spectra that are globally featureless but should normally show strong hydrogen/helium lines if non-magnetic, assuming the effective temperature inferred by their photometric energy distribution.
We believe these "Hot DC" are probably highly magnetic white dwarfs with line features that are so spread out that almost no spectral features can be recognized anymore.

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/MAGLINHE}
    \caption{\label{f:helium} Position of helium lines between 3000 and 8000 \AA{} with increasing magnetic field. We show only lines with zero-field wavelengths below 12 000 \AA{} for clarity.}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.4\linewidth]{figures/representative_synth_spectra}
    \caption{\label{f:repr_synth} Representative reconstructed synthetic spectra of magnetic DB white dwarfs at various field strength, inclination and offsets for an effective temperature of 20,000 K and a surface gravity of $\log{g}=8.0$. The spectra are normalized by the flux value of the continuum at 3500 \AA{} and offset by an arbitrary factor for clarity.}
\end{figure}

The methodology to obtain the atmospheric parameters is the same as in Paper I.
Briefly, photometry is first fit using multiple model grids with varying magnetic field intensities, and the resulting atmospheric parameters are then used to determine the dipolar geometry that best reproduces the spectroscopic observations.
The magnetic field geometry is reconstructed from pre-calculated grid (13 GB of data) covering 5 parameters: $\log{g}=7.0$ to $9.5$ (in steps of 0.5 dex), $T_{\rm{eff}}=$12,000 K to 30,000 K (in steps of 3,000 K), field intensity (0 to 500 MG in 54 log steps), the angle between the line of sight and the normal of the surface ($\psi$ from 0 to $\pi$ in 6 steps), and the limb darkening ($\mu = \cos{\theta}$ from 0 to 1 in 4 steps).
While this exploratory work only considers offset dipole geometries, our technique allows us to rapidly construct any desired geometry.

\section{Results} \label{s:results}

Using the fitting technique described in detail in Paper I, we fitted all 79 objects in our sample in a homogeneous way assuming single DB synthetic spectra with a dipolar geometry.
Each solution was then carefully inspected by eye to judge the quality of the fit.
From this exercise, we found that, similarly to many DAH from Paper I, the DBH classification for 18 objects did not hold up against detailed comparison with magnetic synthetic spectra.
The alleged splitting observed in some lines was probably, in hindsight, simply noise.
While better spectroscopic observation may reveal that some of these stars are magnetic after all, we believe classifying these objects as magnetic based on these low signal-to-noise ratio spectra was perhaps a bit premature.
We thus reject these objects from our analysis and reclassify them as DB until proof to the contrary is obtained.
Five other objects were also rejected since the spectra were too noisy to assess their magnetic nature.
Furthermore, fit of the energy distribution for 20 objects indicated an effective temperature below $\sim$12,000K, too cool to show helium features.
Again, noise in the spectra was probably interpreted as a sign of magnetism for most of these. Some of these objects may turn to genuinely be magnetic white dwarfs when spectropolarimetric observations become available, but no such assessment could be made based on the spectra at hand.
The list of white dwarfs previously classified as magnetic that we reject for our analysis is presented in \autoref{t:rejected}.

\input{rejected_table}

From our fit assuming offset dipole geometry, we could only firmly confirm magnetism in 8 DB stars in our sample (Figures \ref{f:good_fit_0} to \ref{f:good_fit_7}).
With the exception of J0043-1000 \citep[Feige 7][]{Achilleos1992} these are, to our knowledge, the very first magnetic DB to be modeled with detailed helium line splitting calculations.
Dipole solutions appear to be an excellent representation for those stars except for J0043-1000 and J0856+1611.
In the first case, the small discrepancies that are observed can be explained by the fact that our grid does not include small traces of hydrogen that are known to be present (hints of H$\alpha$ can also be observed in J0305+3747 and J2323-0046).
Furthermore, \citet{Achilleos1992} also showed that time resolved spectroscopy revealed composition changes across the surface as the star rotates, something that is not captured by the single spectra used in this study.
In the case of J2323-0046, the predicted depth of the helium lines does not agree with the observations.
This is reminiscent of the shallow features found for numerous magnetic DA in Paper I, although the differences are more subtle here. It is also possible magnetic spots on the surface produce temperature variations that may explain why the lines are not as deep as expected from its photometric atmospheric parameters.
\autoref{t:results} present the atmospheric parameters for the confirmed magnetic DB stars in our sample.
We note that two objects have surface gravities lower than $\log{g}$=7, the lowest value in our grid.
Although this is a very small number of individuals, the average mass is around \input{average_mass}.
It thus appears that magnetic helium rich white dwarfs can be found at low, normal and high mass range, offering little indication of preferred formation channel.

\input{figures/mwd_db.tex}

\input{results_table}

The rest of our sample comprises stars that we are unable to model using our magnetic helium models.
In fact, with the exception of the confirmed magnetic DB with a relatively small field strength J0212+0644 (\autoref{f:good_fit_4}), none of the stars that we were unable to fit (or that we tentatively labeled as Magnetic non DA) in Paper I could successfully be reproduced using the current theoretical framework.
Likewise, we were also unsuccessful in modeling the few objects that are known to exhibit spectroscopic signatures (stationary components) indicative of neutral helium in very high magnetic field \citep[][]{Jordan1998, Wickramasinghe2002}, GD 229, HE 1211$-$1707, and HE 1043$-$0502.
These objects possibly harbor extremely strong magnetic fields that cannot be simply reproduced assuming a simple offset dipole and other geometry will need to be explored in future studies.
Finally, we find several additional objects (see \autoref{t:nonDB}) showing a number of features that could also not be reproduced with either the hydrogen or helium grid that we computed.
While our work ruled out a simple dipole geometry for both hydrogen and helium-rich composition, it is possible that a much more complex field structure is at play.
Another possibility is that some of the absorption features are from another element than hydrogen and helium, carbon and oxygen being the prime suspects.
\citet{Dufour2008, Dufour2010} has shown that a large fraction of the so-called Hot DQ (see \autoref{f:hotdq} for an example spectrum) are magnetic at some level (a few hundred kG to 2-3MG).
Similar objects but with much stronger fields could be very difficult to identify in the absence of detailed calculation for the position of the various transitions as a function of magnetic field strength.
It is interesting to note, for instance, that many objects show a significant absorption feature around 4200 \AA{} (for example \autoref{f:4200A}).
It is tempting to attribute this to carbon as Hot DQ white dwarf tend to show a similar strong feature but at a slightly longer wavelength (see \autoref{f:hotdq}).
This hypothesis would necessitate that all other features would be smeared out.
However, in the absence of calculation of carbon/oxygen line position in very large magnetic fields, we refrain from drawing any conclusion at this point.

\input{magnetic_non_DB_table}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{figures/fits/SDSS_J140750.65+301130.2}
    \caption{\label{f:4200A} Star J1407+3011 (SDSS J140750.65+301130.2) where the spectrum is mostly featureless, with the exception of one unidentified (and not reproduced) feature near 4200 \AA{}.}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{figures/hot_dq_example}
    \caption{\label{f:hotdq} An example of carbon/oxygen dominated atmosphere (Hot DQ) with a strong feature near 4200\AA{}, J1153+0056 \citep[WD 1150+012,][]{Dufour2008}.}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{figures/fits/US_736}
    \caption{\label{f:hotdc} Star J0935+4429 (US 736) where the spectrum is almost featureless, and the offset dipole solution leads to minimal spectral features.}
\end{figure}

Finally, we note the presence of several objects that show almost continuous spectra.
This is very unusual, as both DA and DB white dwarfs should show strong lines in the optical at their inferred effective temperature (see \autoref{t:hotdc}).
It is, of course, difficult if not impossible to infer their chemical composition from spectrometry in these cases. We believe, however, that these "Hot DC" have geometries that produce surface field strength variations that result in spectral components almost completely smeared across the spectra.
\autoref{f:hotdc} shows an example of an offset dipole solution that leads to practically no spectral features.
It is important to note that, while the magnetic field geometries do reproduce the observations to an extent, no physical information about the field geometry should be extracted from those solutions.

\input{hotdc_table}

\section{Conclusion} \label{s:conclu}

In this study, we performed the first homogeneous analysis of all known magnetic white dwarfs showing splitting in neutral helium lines reported in the literature.
After a careful examination of all candidates, we could not confirm the magnetic nature of 18 objects and reclassify them as simple DB.
The fact that the average mass of those stars is \Msol{M_*  = 0.61} lends support to the idea that these white dwarfs were misclassified as magnetic.
Assuming an offset dipole geometry for the magnetic field, we could well reproduce the spectra of 8 magnetic DB white dwarfs.
While we were able to reproduce the observed spectra of these few stars to a very high level of detail, there are still some where we cannot reach such precision, even though the presence of stationary components seems to indicate the presence of neutral helium lines.
It is likely that these stars have a much more complex magnetic field structure.
It is also possible that some stars have a different chemical composition, possibly carbon and oxygen, but the lack of detailed calculation of the line strength and position under very large magnetic field prevent a formal attribution.
We also identify featureless hot objects that are likely highly magnetic white dwarfs with lines spread beyond recognition.
Interestingly, these stars are also among some of the most massive white dwarfs of our sample (\Msol{M_*  > 1.2}).
Although the number of individuals is rather low, magnetic DB white dwarfs also appear to be more massive, on average, than typical DB white dwarfs \citep[][]{Bergeron2011, GenestBeaulieu2019}, suggesting, as has already been proposed by several authors for their hydrogen-rich counterpart, that they could be the result of mergers.

Future work on such objects should consider higher multipole expansion and other chemical composition when the appropriate data become available.
Also, given that rapid rotation is often found among highly magnetic white dwarfs, time resolved spectroscopic observations should be obtained in order to get clues about the exact field geometry.
Finally, it is still unknown how exactly the atmospheric structure and line broadening are affected in the presence of a very large magnetic field.
While neglecting these effects does not prevent us from achieving relatively good fit at moderate fields, those approximations may be responsible in part for our failure to model objects with much higher field strength.
Theoretical calculations exploring those aspects are highly awaited.

\section*{Data Availability}

The data underlying this article are available on the Montreal White Dwarf Database\footnote{\url{https://www.montrealwhitedwarfdatabase.org/}} \citep{MWDD}.

\bibliographystyle{apj}