\section{Introduction}\label{sect:intro}
Let us consider the model quasilinear elliptic problem
\begin{equation}\label{eq:p-laplacian}
\begin{cases}
-\Delta_{p} u  = \sigma u^{q}, \quad u\ge 0  \quad \text{in} \;\; \mathbb{R}^n, \\
\liminf\limits_{x\rightarrow \infty}u(x) = 0,
\end{cases}
\end{equation}
in the sub-natural growth case $0<q< p-1$.

Here $\Delta_{p}u = \text{div}( |\nabla u|^{p-2} \nabla u )$ is the $p$-Laplacian with
$1<p<\infty$, and $\sigma$ is a nonnegative locally integrable function, or
more generally, nonnegative Radon measure on $\mathbb{R}^{n}$ (in brief, $ \sigma \in \mathcal{M}^{+}(\mathbb{R}^n)$).

In this paper, we use  the notion of $p$-superharmonic solutions, or equivalently,
locally renormalized solutions. We refer to  \cite{HKM, KKT} for the corresponding definitions. Throughout, we tacitly assume that $u\in L^{q}_{loc}(\Omega, d\sigma)$, so that the right-hand side of \eqref{eq:p-laplacian} is a Radon measure.

We establish a sharp sufficient condition on $\sigma$, in terms of integrability of
nonlinear potentials, for the existence of a positive  solution $u \in L^{r}(\mathbb{R}^n)$, $0<r<\infty$,
to problem \eqref{eq:p-laplacian}.

The Wolff potential (or, more precisely, Havin-Maz'ya-Wolff potential, see \cite{HM, HW, KM})
${\bf{W}}_{\alpha, p}\sigma$ of a measure
$\sigma \in \mathcal{M}^{+}(\mathbb{R}^n)$ is defined, for $1<p<\infty$ and
$0<\alpha<\frac{n}{p}$, by
\[
{\bf{W}}_{\alpha, p}\sigma(x) = \int_{0}^{\infty} \left[ \frac{\sigma(B(x,r))}{r^{n-\alpha p}}\right]^{\frac{1}{p-1}}\; \frac{dr}{r}, \quad x \in \mathbb{R}^n
\]
where $B(x,r) = \{ y \in \mathbb{R}^n : |x-y|< r \}$ is a ball centered at $x \in \mathbb{R}^n$ of radius $r >0$.
Notice that ${\bf{W}}_{\alpha, p}\sigma=+\infty$ for $\alpha \ge \frac{n}{p}$ unless
$\sigma = 0$.

In the linear case $p=2$, it reduces (up to a normalization constant) to the Riesz potential
of order $2\alpha$, ${\bf{W}}_{\alpha, 2}\sigma = {\bf{I}}_{2\alpha}\sigma$, and in particular,
${\bf{W}}_{1, 2}\sigma ={\bf{I}}_{2}\sigma$ is the Newtonian potential.
The reader may wish to consult, for example, \cite{AH, HKM, KM, MZ, Maz} for a comprehensive discussion
of nonlinear  potentials and their applications to partial differential equations.

We observe that bilateral pointwise estimates of nontrivial (minimal) solutions $u$ to
\eqref{eq:p-laplacian} in the sub-natural growth case were obtained in \cite{CV2}:
\begin{equation} \label{two-sided}
c^{-1} [({\bf{W}}_{1, p} \sigma)^{\frac{p-1}{p-1-q}}+\mathbf{K}_{1, p, q}  \sigma]
\le u \le c [({\bf{W}}_{1, p} \sigma)^{\frac{p-1}{p-1-q}} +\mathbf{K}_{1, p, q}  \sigma],
\end{equation}
where $c>0$ is a constant which  depends only on $p$, $q$, and $n$.

Here $\mathbf{K}_{1, p, q}$ is the so-called \textit{intrinsic} Wolff potential
associated with \eqref{eq:p-laplacian} introduced in \cite{CV2}. Thus, a necessary
and sufficient condition for the existence of a solution $u \in L^r(\mathbb{R}^n)$
to \eqref{eq:p-laplacian} is given by:
\begin{equation} \label{cond-two-sided}
{\bf{W}}_{1, p}  \sigma \in L^{\frac{r(p-1)}{p-1-q}}(\mathbb{R}^n) \quad \textrm{and} \quad \mathbf{K}_{1, p, q}  \sigma \in L^r(\mathbb{R}^n).
\end{equation}

However, $\mathbf{K}_{1, p, q} \sigma$ is defined in terms of localized embedding
constants in certain weighted norm inequalities, which  makes the second condition  in \eqref{cond-two-sided} quite complicated, and difficult to verify.

In this paper, our goal is to give convenient sufficient conditions for the existence of
solutions $u \in L^r(\mathbb{R}^n)$ in terms of potentials ${\bf{W}}_{1, p} \sigma$. These conditions yield simpler, but more restrictive sufficient conditions
 of the type $\sigma \in L^s(\mathbb{R}^n)$, which are known to be sharp in the scale of Lebesgue spaces (see \cite{BO, CV3}).

Our approach is also applicable to  similar existence results for the fractional Laplace problem
\begin{equation} \label{eq:frac-laplacian}
\begin{cases}
\left(-\Delta \right)^{\alpha} u = \sigma u^{q}, \quad u \ge 0  \quad \text{in} \;\; \mathbb{R}^n, \\
\liminf\limits_{x \rightarrow \infty}u(x) = 0,
\end{cases}
\end{equation}
where $0<q<1$ and $\left(-\Delta \right)^{\alpha}$ is the fractional
Laplacian with $0< \alpha < \frac{n}{2}$.

In the classical case $\alpha = 1$ (or equivalently, $p=2$)
our methods are adapted to establish a sufficient condition for the existence of a positive
superharmonic solution $u \in L^{r}(\Omega)$,
$0<r<\infty$, to the problem
\begin{equation} \label{eq:laplacian}
\begin{cases}
-\Delta u = \sigma u^{q},  \quad u \ge 0 \quad \text{in} \;\; \Omega, \\
\liminf\limits_{x \rightarrow y} u(x) = 0, \quad y \in \partial\Omega,
\end{cases}
\end{equation}
where $0<q<1$ and $\Omega \subseteq \mathbb{R}^n$ is an arbitrary domain (possibly unbounded)
which possesses a positive Green function.

We observe that the existence and uniqueness of {\it bounded} solutions to \eqref{eq:laplacian} on
$\Omega = \mathbb{R}^n$ were established by Brezis and Kamin in \cite{BK}.

In \cite{CV1} (see also \cite{SV1}), it was proved that the condition
\begin{equation}\label{cond:dsigma-wolff}
{\bf{W}}_{1,p}\sigma \in L^{\frac{(\gamma+q)(p-1)}{p-1-q}}(\mathbb{R}^n, d\sigma)
\end{equation}
with $\gamma =1$, is necessary and sufficient for the existence of a positive {\it finite energy} solution $u \in \dot{W}_{0}^{1,p}(\mathbb{R}^n)$ to \eqref{eq:p-laplacian}. Moreover, such a solution
is unique in $\dot{W}_{0}^{1,p}(\mathbb{R}^n)$.

Here, for $1 \leq p<\infty$ and a nonempty open set $\Omega \subseteq \mathbb{R}^n$,
by $\dot{W}^{1,p}_{0}(\Omega)$ we denote the homogeneous Sobolev (or Dirichlet)
space defined \cite{HKM} as the closure of $C_{0}^{\infty}(\Omega)$ with respect to the (semi) norm
\[
\| u \|_{\dot{W}^{1,p}_{0}(\Omega)} =  \| \nabla u \|_{L^{p}(\Omega)}.
\]

In this  study, we will show that for a given $0<r<\infty$,
condition \eqref{cond:dsigma-wolff}
with a suitable value of $\gamma=\gamma(r)>0$ yields the existence of a positive
$p$-superharmonic solution $u \in L^{r}(\mathbb{R}^n)$ to \eqref{eq:p-laplacian}.

Our main results are new even in the case $p=2$, or when $\sigma$ is a nonnegative locally integrable function on $\mathbb{R}^n$.

\begin{Thm}\label{thm:main-p-laplacian}
Let $1<p<n$, $0<q<p-1$, $\frac{n(p-1)}{n-p}<r<\infty$, and $\sigma \in \mathcal{M}^{+}(\mathbb{R}^n)$ with $\sigma \not\equiv 0$. If \eqref{cond:dsigma-wolff} holds with
$\gamma = \frac{r(n-p)-n(p-1)}{n}$, then there exists a positive $p$-superharmonic solution $u \in L^r(\mathbb{R}^{n})$ to \eqref{eq:p-laplacian}.

If $n \leq p < \infty$, or $1<p<n$
and $0<r\le \frac{n(p-1)}{n-p}$, then
there is only a trivial supersolution to \eqref{eq:p-laplacian}.
\end{Thm}

A necessary  (but generally not sufficient) condition
for the existence of a nontrivial solution $u \in L^{r}(\mathbb{R}^n)$ to \eqref{eq:p-laplacian},
where $1<p<n$ and $\frac{n(p-1)}{n-p}<r<\infty$,  follows from \eqref{cond-two-sided} (see also Theorem \ref{thm:lower-est-p} below),
\begin{equation}\label{cond:dx-wolff-1}
{\bf{W}}_{1,p} \sigma \in L^{\frac{r(p-1)}{p-1-q}}(\mathbb{R}^n).
\end{equation}
Here we use $L^s(\mathbb{R}^n)$  with respect to Lebesgue measure, whereas in condition \eqref{cond:dsigma-wolff}, $L^s (\mathbb{R}^n, d \sigma)$ is with respect to the measure
$\sigma$.

A simple sufficient condition for \eqref{cond:dsigma-wolff} with $\gamma =
\frac{r(n-p)-n(p-1)}{n}$
is given by
\begin{equation}\label{cond:sigma-suff}
\sigma \in L^{s_1}(\mathbb{R}^n),   \quad s_{1}=\tfrac{nr}{n(p-1-q)+ p r},
\end{equation}
where  $1<p<n$ and $\frac{n(p-1)}{n-p}<r<\infty$ (see Proposition \ref{prop:wolff} with $\alpha=1$).
Thus, in light of Theorem \ref{thm:main-p-laplacian}, we have the following result.

\begin{Cor}\label{cor:main1}
Under the assumptions of Theorem \ref{thm:main-p-laplacian}, if
 \eqref{cond:sigma-suff} holds, then there exists a positive $p$-superharmonic solution
$u \in L^r(\mathbb{R}^n)$ to \eqref{eq:p-laplacian}, provided $1<p<n$ and $\frac{n(p-1)}{n-p}<r<\infty$.
\end{Cor}

\begin{Rem}\label{rem-main1}
 The sufficient conditions of Theorem \ref{thm:main-p-laplacian} and Corollary
\ref{cor:main1}, along with the necessary condition \eqref{cond:dx-wolff-1}, hold for the $\mathcal{A}$-Laplacian in place of $\Delta_{p}$, under the standard structural  assumptions on $\mathcal{A}$ (see \cite{CV2, HKM, MZ}).
\end{Rem}

The main new feature in our approach is the use of
two weight  inequalities involving Wolff potentials (see Lemma
\ref{lem1-wolff}  and Lemma \ref{lemma-link-p} below).

Our methods in the case $p=2$ are applicable to analogous existence results for the fractional Laplace problem
\eqref{eq:frac-laplacian}.

\begin{Thm}\label{thm:main2}
Let $0<q<1$, $0<\alpha<\frac{n}{2}$ and $\sigma \in \mathcal{M}^{+}(\mathbb{R}^n)$ with
$\sigma \not\equiv 0$. Suppose that $\frac{n}{n-2\alpha}<r<\infty$ and
\begin{equation}\label{cond:dsigma-riesz}
{\bf{I}}_{2\alpha}\sigma \in L^{\frac{\gamma+q}{1-q}}(\mathbb{R}^n, d\sigma),
\end{equation}
where $\gamma = \frac{r(n-2\alpha)-n}{n}$. Then there exists a positive  solution $u \in L^r(\mathbb{R}^{n})$ to \eqref{eq:frac-laplacian}.

If $0<r\le \frac{n}{n-2\alpha}$, then
there is only a trivial supersolution to \eqref{eq:frac-laplacian}.
\end{Thm}

Here, a solution $u\ge 0$ to problem \eqref{eq:frac-laplacian} is understood in the sense that
$u \in L^{q}_{loc}(\mathbb{R}^n, \sigma)$ satisfies the corresponding
integral equation
\begin{equation}\label{eq:int2}
u = {\bf{I}}_{2\alpha}(u^{q} d\sigma) \quad \text{in} \;\; \mathbb{R}^{n}.
\end{equation}

A necessary condition for the existence of such a solution $u \in L^{r}(\mathbb{R}^n)$,  to \eqref{eq:frac-laplacian}, which follows from Theorem \ref{lower_ptwise_est} below with $p=2$, is given by
\begin{equation}\label{cond:dx-riesz-1}
{\bf{I}}_{2\alpha} \sigma \in L^{\frac{r}{1-q}}(\mathbb{R}^n).
\end{equation}

As above, a simple sufficient condition for \eqref{cond:dsigma-riesz} with $\gamma = \frac{r(n-2\alpha)-n}{n}$
is provided by
\begin{equation}\label{cond:sigma-suff2}
\sigma \in L^{s_2}(\mathbb{R}^n), \quad  s_{2}=\tfrac{nr}{n(1-q)+2\alpha  r},
\end{equation}
where $\frac{n}{n-2\alpha}<r<\infty$.

\begin{Cor}\label{cor:main2}
Under the assumptions of Theorem \ref{thm:main2},
if \eqref{cond:sigma-suff2} is valid, then there exists a positive solution
$u \in L^r(\mathbb{R}^n)$ to \eqref{eq:frac-laplacian}.
\end{Cor}

We observe that  condition \eqref{cond:dsigma-riesz}
with $\gamma=1$ is necessary and sufficient for the existence of a positive finite energy solution
$u \in \dot{H}^{\alpha}(\mathbb{R}^n)$, $0<\alpha <\frac{n}{2}$, to \eqref{eq:frac-laplacian};  uniqueness of such a solution was proved in the case $0<\alpha\leq1$, see \cite{SV1}.

Here, the homogeneous Sobolev space $\dot{H}^{\alpha}(\mathbb{R}^n)$ ($0<\alpha<\frac{n}{2}$) can be defined by means of Riesz potentials,
\[
\dot{H}^{\alpha}(\mathbb{R}^n) =  \big\lbrace u\colon \,  u = {\bf{I}}_{\alpha}f,\; f \in L^{2}(\mathbb{R}^n) \big\rbrace,
\]
equipped with  norm  $\|u \|_{\dot{H}^{\alpha}(\mathbb{R}^n)}= \| f \|_{L^{2}(\mathbb{R}^n)}$.

Finally, we consider the sublinear elliptic problem \eqref{eq:laplacian} on an arbitrary domain
$\Omega \subseteq \mathbb{R}^n$ (possibly unbounded) with positive Green function
$G(x,y)$ on $\Omega \times \Omega$.
The corresponding Green potential of a measure $\sigma \in \mathcal{M}^{+}(\Omega)$ is defined by
\[
{\bf{G}}\sigma(x) = \int_{\Omega} G(x,y) \, d\sigma(y), \quad x \in \Omega.
\]
Our main results for this problem read as follows.

\begin{Thm}\label{thm:main3}
Let $0<q<1$, $\sigma \in \mathcal{M}^{+}(\Omega)$ with
$\sigma \not\equiv 0$, and let $G$ be the positive Green function
associated with $-\Delta$ in $\Omega \subseteq \mathbb{R}^n$, $n \ge 3$.
Suppose $\frac{n}{n-2}<r<\infty$ and
\begin{equation}\label{cond:dsigma-green}
{\bf{G}}\sigma \in L^{\frac{\gamma+q}{1-q}}(\Omega, d\sigma),
\end{equation}
where $\gamma = \frac{r(n-2)-n}{n}$. Then there exists a positive superharmonic solution
$u \in L^r(\Omega)$ to \eqref{eq:laplacian}.
\end{Thm}

Invoking Theorem \ref{Thm:lowerbound} below, in the special case where $G$ is a positive Green function associated with $-\Delta$, yields a necessary condition for the existence of such a solution,
\begin{equation}\label{cond:dx-green}
{\bf{G}} \sigma \in L^{\frac{r}{1-q}}(\Omega).
\end{equation}

A sufficient condition for \eqref{cond:dsigma-green} with $\gamma = \frac{r(n-2)-n}{n}$  is given by
\begin{equation}\label{cond:sigma-suff3}
\sigma \in L^{s_3}(\Omega), \quad  s_{3}=\tfrac{n r}{n(1-q)+2r},
\end{equation}
where $\frac{n}{n-2}<r<\infty$
(see Proposition \ref{prop3}, when $G$ is a positive Green function associated with
$-\Delta$).
Thus, the following corollary is easily deduced from Theorem \ref{thm:main3}.

\begin{Cor}\label{cor:main3}
Under the assumptions of Theorem \ref{thm:main3},
if condition \eqref{cond:sigma-suff3} holds, then there exists a positive superharmonic solution
$u \in L^r(\Omega)$ to \eqref{eq:laplacian}.
\end{Cor}

Corollary \ref{cor:main3}, when $\Omega \subseteq \mathbb{R}^n$ is a  bounded domain, is due to Boccardo and Orsina \cite{BO}, with a different proof.

It was recently shown by the authors \cite{SV2} that, for an arbitrary $0<\gamma<\infty$, condition
\eqref{cond:dsigma-green} is necessary and sufficient for the existence a positive
superharmonic solution $u$ to \eqref{eq:laplacian} with {\it finite generalized energy}:
\[
\mathcal{E}_{\gamma}[u] = \int_{\Omega} |\nabla u|^{2} u^{\gamma - 1}\,dx < +\infty.
\]
Theorem \ref{thm:main3}  can be deduced from this result, but we
give an independent proof below.

\begin{Rem}\label{rem3}
The above results hold, with the same proofs,  for  linear
uniformly elliptic operators with bounded measurable coefficients
$L=-\text{div}(\mathcal{A} \nabla \cdot)$ in place of $-\Delta$, on an
arbitrary domain $\Omega \subseteq \mathbb{R}^n$ with positive Green function, as well as
more general (nonlocal) integral equations $u= \mathbf{G}(u^q d\sigma)$
with positive kernel $G$.

The restrictions on $G$ imposed below are the {\it weak maximum principle} (WMP), and {\it quasi-symmetry}
(see Sec. \ref{sect:kernel}). We also sometimes add  the condition that $G$ is bounded above by the Riesz kernel
$I_{2\alpha}(x-y)= |x-y|^{2\alpha-n}$
of order $2\alpha$ for some $0<\alpha<\frac{n}{2}$,
\begin{equation}\label{cond:Riesz0}
G(x, y) \leq c I_{2\alpha}(x-y), \quad  \forall x, y \in \Omega,
\end{equation}
where $c$ is a positive constant independent of $x$ and $y$.
\end{Rem}

This paper is organized as follows.
In Sec. 2, we recall the necessary background and state preliminary results
concerning nonlinear potentials and related weighted norm inequalities.
In Sec. 3, we simultaneously prove the existence results for positive solutions
$u \in L^{r}(\mathbb{R}^n)$
to both problems \eqref{eq:p-laplacian} and \eqref{eq:frac-laplacian}.
Our main results for problem \eqref{eq:laplacian} will be established in Sec. 4.

\section{Preliminaries}
\subsection{Wolff potentials}\label{sect:wolff}
Recall that for $1<p<\infty$ and $0<\alpha<\frac{n}{p}$,
the (homogeneous) potential ${\bf{W}}_{\alpha, p}\sigma$ of a measure $\sigma \in \mathcal{M}^{+}(\mathbb{R}^n)$  is defined by \cite{HM, HW}
\[
{\bf{W}}_{\alpha, p}\sigma(x) := \int_{0}^{\infty} \left[ \frac{\sigma \left( B(x, r) \right)}{r^{n-\alpha p}}\right]^{\frac{1}{p-1}}\; \frac{dr}{r}, \quad x \in \mathbb{R}^n,
\]
where $B(x,r):= \{ y \in \mathbb{R}^n : |x-y|< r \}$ is a ball centered at $x \in \mathbb{R}^n$ of radius
$r >0$.

It is well-known that
\begin{equation}\label{inq:HM}
{\bf{W}}_{\alpha, p}\sigma(x) \leq c \mathbf{V}_{\alpha, p}\sigma(x), \quad x \in \mathbb{R}^n,
\end{equation}
where $\mathbf{V}_{\alpha, p}\sigma := {\bf{I}}_{\alpha}[({\bf{I}}_{\alpha}\sigma)^{\frac{1}{p-1}}dx]$ is
the Havin-Maz'ya potential, and $c=c(\alpha, p, n) $ is a positive constant.

Havin-Maz'ya potentials satisfy the weak maximum (or boundedness)
principle (see \cite[Theorem 2.6.3]{AH}). A similar weak maximum principle holds
for Wolff potentials: if $\sigma \in \mathcal{M}^{+}(\mathbb{R}^n)$, then
\begin{equation}\label{wolff-max}
{\bf{W}}_{\alpha, p}\sigma(x) \leq c  \, \sup \left \{{\bf{W}}_{\alpha, p}\sigma(y): \, y \in \textrm{supp}(\sigma)\right\}, \quad \forall x \in \mathbb{R}^n,
\end{equation}
where  $c=c(\alpha, p, n)$ is a positive constant.

The following pointwise iterated inequality for Wolff potentials will be used below.

\begin{Thm}[\cite{CV1}]\label{Thm:iterated-wolff}
Let $1<p< \infty$, $0<\alpha< \frac{n}{p}$, and $\sigma \in \mathcal{M}^{+}(\mathbb{R}^{n})$
with $\sigma \not\equiv 0$. Then the following estimates hold:
\begin{itemize}
	\item[(i)]If $t\geq1$, then
	\begin{equation}\label{iterated1-wolff}
({\bf{W}}_{\alpha, p}\sigma)^{t}(x) \leq c \, {\bf{W}}_{\alpha, p} \left( ({\bf{W}}_{\alpha, p}\sigma)^{(t-1)(p-1)} d\sigma \right)(x), \quad  \forall x \in \mathbb{R}^{n}.
	\end{equation}
	\item[(ii)] If $0<t\leq1$, then
	\begin{equation}\label{iterated2-wolff}
({\bf{W}}_{\alpha, p}\sigma)^{t}(x) \geq c \, {\bf{W}}_{\alpha, p} \left( ({\bf{W}}_{\alpha, p}\sigma)^{(t-1)(p-1)} d\sigma \right)(x), \quad  \forall x \in \mathbb{R}^{n}.
	\end{equation}
\end{itemize}
Here, $c= c(t, p, n, \alpha) > 0$.
\end{Thm}

The following weak continuity result will be used to prove the existence of
$p$-superharmonic solutions to quasilinear equations. We denote by
$\omega[u]$ the $p$-measure associated with a $p$-superharmonic function $u$ (see \cite{HKM}).

\begin{Thm} [\cite{TW}] \label{weak_cont_p-Laplacian}
Suppose $\lbrace u_{j} \rbrace_{j=1}^{\infty}$ is a sequence of nonnegative $p$-superharmonic functions in $\Omega$ such that $u_j \rightarrow u$ a.e. as $j \rightarrow \infty$, where $u$ is a $p$-superharmonic function in $\Omega$. Then $\omega[u_{j}]$ converges weakly to $\omega[u]$, that is,
\[
\lim_{j \rightarrow \infty} \int_{\Omega} \varphi \; d\omega[u_j]
= \int_{\Omega} \varphi \; d\omega[u]
\]
for all $\varphi \in C_{0}^{\infty}(\Omega)$.
\end{Thm}

We shall use the following lower bounds for supersolutions to \eqref{eq:p-laplacian}.

\begin{Thm}[\cite{CV1}] \label{thm:lower-est-p}
Let $1<p<n$, $0<q<p-1$ and $\sigma \in \mathcal{M}^{+}(\mathbb{R}^n)$.
Suppose $u$ is a nontrivial supersolution to \eqref{eq:p-laplacian}. Then $u$ satisfies the inequality
\[
u \geq c \left( {\bf{W}}_{1, p}\sigma \right)^{\frac{p-1}{p-1-q}},
\]
where $c = c(n, p, q) >0$. If $p \geq n$ there is only a trivial supersolution $u=0$ on $\mathbb{R}^n$.
\end{Thm}

\begin{Thm} [\cite{CV2}] \label{lower_ptwise_est}
Let $1<p<n$, $0<q<p-1$, $0<\alpha<\frac{n}{p}$ and $\sigma \in \mathcal{M}^{+}(\mathbb{R}^n)$. Suppose $u \in L^{q}_{loc}(\mathbb{R}^n, d\sigma)$, $u\ge 0$,  is a nontrivial
function satisfying
\[
u \geq  {\bf{W}}_{\alpha, p} (u^{q}d\sigma) \quad \text{in} \;\; \mathbb{R}^n.
\]
Then $u$ satisfies the inequality
\begin{equation}\label{inq:lowerbound}
u \geq c \left( {\bf{W}}_{\alpha, p}\sigma \right)^{\frac{p-1}{p-1-q}} \quad \text{in} \;\; \mathbb{R}^n,
\end{equation}
where $c = c(\alpha, n, p, q)$ is a positive constant.
\end{Thm}

The following important result is concerned with pointwise estimate of nonnegative $p$-superharmonic functions in terms of Wolff's potential.

\begin{Thm} [\cite{KM}]\label{pointwise_est_p-superharmonic}
Let $1<p<n$ and $\sigma \in \mathcal{M}^{+}(\mathbb{R}^n)$ Suppose
$u$ is a $p$-superharmonic function in $\mathbb{R}^n$
satisfying
\[
\begin{cases}
-\Delta_{p} u = \sigma \quad \text{in} \;\; \mathbb{R}^n, \\
\liminf\limits_{\vert x \vert \rightarrow \infty} u(x) = 0.
\end{cases}
\]
Then
\[
K^{-1}{\bf{W}}_{1,p} \sigma \leq u \leq K {\bf{W}}_{1,p}\sigma,
\]
where $K=K(n,p) \geq 1$.
\end{Thm}
\subsection{Kernels and potentials}\label{sect:kernel}
Let $G: \Omega \times \Omega \rightarrow (0, \infty]$ be a positive lower semicontinuous kernel. The potential of a measure $\sigma \in \mathcal{M}^{+}(\Omega)$ is defined by
\[
{\bf{G}} \sigma (x) := \int_{\Omega} G(x,y) \;d\sigma (y), \quad x \in \Omega.
\]

Note that ${\bf{G}}\sigma(x)$ is lower semi-continuous on $\Omega \times \mathcal{M}^{+}(\Omega)$ if and only if $G(x,y)$ is lower semi-continuous on $\Omega \times \Omega$ (see \cite{Br}).

A positive kernel $G$ on $\Omega \times \Omega$ is said to satisfy the weak maximum principle (WMP) with constant $h \geq 1$ if for any $\sigma \in \mathcal{M}^{+}(\Omega)$,
\begin{equation}\label{WMP}
\sup \{ {\bf{G}}\sigma(x) : x \in \text{supp}(\sigma) \} \leq 1
\Longrightarrow
\sup \{ {\bf{G}}\sigma(x) : x \in \Omega \} \leq h.
\end{equation}
Here we use the notation supp$(\sigma)$ for the support of $\sigma$.

When $h=1$ in \eqref{WMP}, the kernel $G$ is said to satisfy the strong maximum principle.  It holds for Green functions associated with the classical Laplacian $-\Delta$, or more
generally the linear uniformly elliptic operator in divergence form $L$, as well as
the fractional Laplacian $(-\Delta)^{\alpha}$ in the case $0<\alpha\leq 1$, in every domain
$\Omega \subset \mathbb{R}^n$ which possesses a Green function.

The WMP holds for Riesz kernels on $\mathbb{R}^n$ associated with $(-\Delta)^{\alpha}$ in the full range
$0< \alpha < \frac{n}{2}$, and more generally for all radially nonincreasing kernels on $\mathbb{R}^n$,
see \cite{AH}.

We say that a positive kernel $G$ on $\Omega \times \Omega$ is quasi-symmetric if there exists
a constant $a\geq 1$ such that
\begin{equation}\label{quasi-sym}
a^{-1}G(y,x) \leq G(x,y) \leq a G(y,x), \quad x,y \in \Omega.
\end{equation}
When $a=1$ in \eqref{quasi-sym}, the kernel $G$ is said to be symmetric.
There are many kernels associated with elliptic operators that are quasi-symmetric and
satisfy the WMP, see \cite{An}.

We begin with the following pointwise estimates (see \cite[Theorem 1.3]{GV}) for supersolutions
to sublinear integral equations such that
\begin{equation}\label{eq:intG}
u (x) \ge  {\bf{G}}(u^{q} d\sigma) (x) \quad \forall x \in \Omega.
\end{equation}
\begin{Thm}[\cite{GV}]\label{Thm:lowerbound}
Let $0<q<1$, $\sigma \in \mathcal{M}^{+}(\Omega)$, and let
$G$ be a positive lower semicontinuous kernel on
$\Omega \times \Omega$, which satisfies the WMP with constant $h \geq 1$.
If $u \in L^{q}_{loc}(\Omega, d\sigma)$ is a positive supersolution satisfying \eqref{eq:intG}, then
\begin{equation}\label{lowerbound}
u(x)  \geq (1-q)^{\frac{1}{1-q}} h^{-\frac{q}{1-q}} [ {\bf{G}} \sigma (x) ]^{\frac{1}{1-q}}, \quad  \forall x \in \Omega.
\end{equation}
\end{Thm}

The following iterated inequalities are used in our argument, see \cite[Lemma 2.5]{GV}.

\begin{Thm}[\cite{GV}]\label{Thm:iterated}
Let $\sigma \in \mathcal{M}^{+}(\Omega)$ with $\sigma \not\equiv 0$, and let $G$ be a positive lower semicontinuous kernel on
$\Omega \times \Omega$, which satisfies the WMP with constant $h\geq1$.
Then the following estimates hold:
\begin{itemize}
	\item[(i)]If $t\geq1$, then
	\begin{equation}\label{iterated1}
({\bf{G}}\sigma)^{t}(x) \leq t h^{t-1} \, {\bf{G}} \left( ({\bf{G}}\sigma)^{t-1} d\sigma \right)(x), \quad  \forall x \in \Omega.
	\end{equation}
	\item[(ii)] If $0<t\leq1$, then
	\begin{equation}\label{iterated2}
({\bf{G}}\sigma)^{t}(x) \geq t h^{t-1} \, {\bf{G}} \left( ({\bf{G}}\sigma)^{t-1} d\sigma \right)(x), \quad  \forall x \in \Omega.
	\end{equation}
\end{itemize}
\end{Thm}

The next result   characterizes explicitly
weighted norm inequalities of the $(s ,r)$-type  in the case $0< r < s$ and $1< s < \infty$,
 for
operators ${\bf{G}}$ (see \cite[Theorem 1.1]{V}):
\begin{equation}\label{weighted_norm_ineq}
\big\| {\bf{G}}(f d\sigma) \big\|_{L^{r}(\Omega,\,d\sigma)} \leq c \| f \|_{L^{s}(\Omega,\,d\sigma)},
\quad \forall  f \in L^{s}(\Omega, d\sigma),
\end{equation}
where $c$ is a positive constant independent of $f$, for an arbitrary measure $\sigma \in \mathcal{M}^{+}(\Omega)$, under certain assumptions on $G$.

\begin{Thm}[\cite{V}]\label{thm_Ver17}
Let $\sigma \in \mathcal{M}^{+}(\Omega)$ with $\sigma \not\equiv 0$, and let $G$ be a positive,
quasi-symmetric, lower semicontinuous kernel on $\Omega \times \Omega$, which satisfies the WMP.
\begin{itemize}
\item[(i)] If $1 < s < \infty$ and $0 < r < s$, then the weighted norm inequality
\eqref{weighted_norm_ineq} holds if and only if
\begin{equation}\label{energy1}
{\bf{G}}\sigma \in L^{\frac{sr}{s-r}}(\Omega, d\sigma).
\end{equation}

\item[(ii)] If $0 < q < 1$ and $0<\gamma<\infty$, then there exists a positive
(super)solution
$u \in L^{\gamma + q}(\Omega, d\sigma)$ to equation \eqref{eq:intG}
if and only if the weighted norm inequality \eqref{weighted_norm_ineq} holds with
$r=\gamma +q$ and $s= \frac{\gamma +q}{q}$, i.e.,
\begin{equation}\label{weighted_norm_ineq_special}
\big\| {\bf{G}}(f d\sigma) \big\|_{L^{\gamma + q}(\Omega,\;d\sigma)}
\leq c \| f \|_{L^{\frac{\gamma + q}{q}}(\Omega, d\sigma)},
\quad \forall f \in L^{\frac{\gamma + q}{q}}(\Omega,  d\sigma),
\end{equation}
or equivalently,
\begin{equation}\label{energy2}
{\bf{G}}\sigma \in L^{\frac{\gamma + q}{1-q}}(\Omega, d\sigma).
\end{equation}
\end{itemize}
\end{Thm}

The existence of solutions $u\in L^{q}(\Omega, d\sigma)$, and more generally
$u\in L^{q}_{{\rm loc}}(\Omega, d\sigma)$,
 is characterized in \cite{QV} via $(1,q)$-weighted norm inequalities, which corresponds to $\gamma=0$ in \eqref{weighted_norm_ineq_special}.

\section{Solutions  to problems \eqref{eq:p-laplacian} and \eqref{eq:frac-laplacian}}
In this section, we prove the existence results stated in the introduction for positive solutions  $u\in L^{r}(\mathbb{R}^{n})$, $0<r< \infty$, to both problems \eqref{eq:p-laplacian} and \eqref{eq:frac-laplacian}.

When $p \geq n$, it readily follows from Theorem \ref{thm:lower-est-p}
that there is only a trivial supersolution to \eqref{eq:p-laplacian}. Henceforth, we assume $1<p<n$.

We begin with an investigation of solvability of the corresponding
nonlinear integral equations involving Wolff potentials,
\begin{equation}\label{eq:int-wolff}
u = {{\bf{W}}}_{\alpha,p} (u^{q}d\sigma) \quad \text{in} \;\; \mathbb{R}^n
\end{equation}
where $1<p<n$, $0<q<p-1$, $0<\alpha<\frac{n}{p}$ and
$\sigma \in \mathcal{M}^{+}(\mathbb{R}^n)$.

The following theorem will be used in the construction of such solutions
to \eqref{eq:p-laplacian} and \eqref{eq:frac-laplacian}
in the cases $\alpha = 1$ and $p=2$, respectively.

\begin{Thm}\label{thm:wolff-existence}
Let $0<\gamma<\infty$, $1<p<n$, $0<q<p-1$, $0<\alpha<\frac{n}{p}$ and $\sigma \in \mathcal{M}^{+}(\mathbb{R}^n)$ such that $\sigma \not\equiv 0$. Then there exists a positive solution  $u \in L^{\gamma+q}(\mathbb{R}^n, d\sigma)$ to
\eqref{eq:int-wolff} if and only if
 \begin{equation}\label{cond:sigma-gamma}
 {{\bf{W}}}_{\alpha, p}\sigma \in L^{\frac{(\gamma+q)(p-1)}{p-1-q}}(\mathbb{R}^n, d\sigma).
 \end{equation}
\end{Thm}

The proof of this theorem is based on the following lemma,
which  is an extension of \cite[Lemma 3.3]{CV1} in the case $\gamma =1$.

\begin{Lem}\label{lem:wolff}
Let $0<\gamma<\infty$, $1<p<\infty$, $0<\alpha<\frac{n}{p}$, $0<q<p-1$, and $\sigma \in \mathcal{M}^{+}(\mathbb{R}^n)$. If \eqref{cond:sigma-gamma} holds then
\begin{equation}\label{winq:wolff-sigma-sigma}
\big\| {\bf{W}}_{\alpha, p}(f d\sigma) \big\|_{L^{\gamma + q}(\mathbb{R}^{n}, d\sigma)}
\leq c \| f \|_{L^{\frac{\gamma+q}{q}}(\mathbb{R}^{n}, d\sigma)}^{\frac{1}{p-1}},
\, \,  \forall f \in L^{\frac{\gamma+q}{q}}(\mathbb{R}^n, d\sigma),
\end{equation}
where $c$ is a positive constant independent of $f$.
\end{Lem}

\begin{proof}
Without loss of generality, we may assume $f \in L^{\frac{\gamma+q}{q}}(\mathbb{R}^n, d\sigma)$
with $f \geq 0$. Observe that
\[
{\bf{W}}_{\alpha, p}(f d\sigma) \leq \left( M_{\sigma}f \right)^{\frac{1}{p-1}} {\bf{W}}_{\alpha, p}\sigma
\]
where the centered maximal operator $M_{\sigma}$ is defined by
\[
M_{\sigma}f (x) := \sup_{r>0} \frac{1}{\sigma(B(x,r))} \int_{B(x,r)} |f| \, d\sigma, \quad x \in \mathbb{R}^{n}.
\]
By H\"older's inequality with exponents  $\frac{p-1}{q}$ and $\frac{p-1}{p-1-q}$, and the
boundedness of the maximal operator $M_{\sigma}: L^{s}(\mathbb{R}^{n}, d\sigma) \rightarrow L^{s}(\mathbb{R}^{n}, d\sigma) $ with $s=\frac{\gamma +q }{q}$, we obtain
\[
\begin{split}
& \int_{\mathbb{R}^n} \left({\bf{W}}_{\alpha, p}(f d\sigma) \right)^{\gamma+q}\,d\sigma
 \leq \int_{\mathbb{R}^n} \left(M_{\sigma}f \right)^{\frac{\gamma+q}{p-1}}
\left({\bf{W}}_{\alpha, p}\sigma \right)^{\gamma+q}\,d\sigma \\
& \leq \left(  \int_{\mathbb{R}^n} \left(M_{\sigma}f \right)^{\frac{\gamma+q}{q}}\,d\sigma \right)^{\frac{q}{p-1}}
\left( \int_{\mathbb{R}^n} \left( {\bf{W}}_{\alpha, p}\sigma \right)^{\frac{(\gamma+q) (p-1)}{p-1-q}}\,d\sigma \right)^{\frac{p-1-q}{p-1}} \\
& \leq c \left( \int_{\mathbb{R}^n} f^{\frac{\gamma+q}{q}}\,d\sigma \right)^{\frac{q}{p-1}}
 \left( \int_{\mathbb{R}^n} \left({\bf{W}}_{\alpha, p}\sigma \right)^{\frac{(\gamma+q)(p-1)}{p-1-q}}\,d\sigma\right)^{\frac{p-1-q}{p-1}}.
\end{split}
\]
Thus, by assumption \eqref{cond:sigma-gamma}, this yields \eqref{winq:wolff-sigma-sigma}.
\end{proof}

We now prove Theorem \ref{thm:wolff-existence}.

\begin{proof}[Proof of Theorem \ref{thm:wolff-existence}] To prove the sufficiency part, we construct a sequence of functions $\lbrace  u_j \rbrace_{j=0}^{\infty}$ by setting
\[
u_0 := c_{0}\left({{\bf{W}}}_{\alpha, p}\sigma\right)^{\frac{p-1}{p-1-q}}
\quad \text{and} \quad
u_{j+1} := {{\bf{W}}}_{\alpha, p} (u^{q}_{j} d\sigma), \quad \;\; j=0,1,2,\dots
\]
where $c_{0}$ is a small positive constant to be determined later.
Observe that $u_{0}>0$ since $\sigma \not\equiv 0$. Moreover,
\[
u_{1} = {{\bf{W}}}_{\alpha, p} (u^{q}_{0} d\sigma)
= c_{0}^{\frac{q}{p-1}} {{\bf{W}}}_{\alpha, p} ( ({{\bf{W}}}_{\alpha, p}\sigma)^{\frac{q(p-1)}{p-1-q}} d\sigma)
\geq c_{0}^{\frac{q}{p-1}}C ({{\bf{W}}}_{\alpha, p}\sigma)^{\frac{p-1}{p-1-q}}
\]
where $C$ is the positive constant in inequality \eqref{iterated1-wolff}
with $t=\frac{p-1}{p-1-q}$.
Choosing $c_{0}$ small enough so that $c_{0}^{\frac{q}{p-1}}C \geq c_{0}$,
we obtain  $u_{0} \leq u_{1}$.
By induction, we see that $\lbrace  u_j\rbrace_{j=0}^{\infty}$ is a nondecreasing sequence of positive functions.

Next, we show that each $u_j \in  L^{\gamma+q}(\mathbb{R}^n, d\sigma)$. Note that assumption \eqref{cond:sigma-gamma} yields
\[
u_0 = c_{0}\left({{\bf{W}}}_{\alpha, p}\sigma\right)^{\frac{p-1}{p-1-q}} \in L^{\gamma+q}(\mathbb{R}^n, d\sigma).
\]

By Lemma \ref{lem:wolff},  weighted norm inequality \eqref{winq:wolff-sigma-sigma} holds.
Suppose that $u_0, \dots, u_j \in L^{\gamma +q}(\mathbb{R}^n, d\sigma)$ for some $j \in \mathbb{N}$.
Applying \eqref{winq:wolff-sigma-sigma} with $f:=u^{q}_{j} \in
L^{\frac{\gamma+q}{q}}(\mathbb{R}^n, d\sigma)$, we deduce
\begin{equation}\label{est2_thm_existence}
\begin{split}
 \Vert u_{j+1} \Vert_{L^{\gamma+q}(\mathbb{R}^n,\,d\sigma)}&  = \big\| {\bf{W}}_{\alpha, p} (u^{q}_{j} d\sigma) \big\|_{L^{\gamma+q}(\mathbb{R}^n,\,d\sigma)} \\
&\leq c \Vert u_{j}\Vert^{\frac{q}{p-1}}_{L^{\gamma+q}(\mathbb{R}^n,\,d\sigma)} \\
&\leq c \Vert u_{j+1}\Vert^{\frac{q}{p-1}}_{L^{\gamma+q}(\mathbb{R}^n,\,d\sigma)}.
\end{split}
\end{equation}
Hence, by induction,
\[
\Vert u_{j+1} \Vert_{L^{\gamma+q}(\mathbb{R}^n,\,d\sigma)}
\leq c^{\frac{p-1}{p-1-q}} < +\infty, \quad j=0, 1, \ldots .
\]
Finally, applying the Monotone Convergence Theorem to $\lbrace u_{j} \rbrace_{j=0}^{\infty}$, we see that the pointwise limit $u=\lim_{j \rightarrow \infty} u_j$ exists so that $u>0$,
$u \in L^{\gamma+q}(\mathbb{R}^n, d\sigma)$, and $u$ satisfies equation \eqref{eq:int-wolff}.

The necessity part follows immediately from the global pointwise lower bound \eqref{inq:lowerbound} of supersolutions to \eqref{eq:int-wolff}.
\end{proof}

Let us consider the following weighted norm inequality with Lebesgue measure on the left-hand side:
\begin{equation}\label{winq:wolff-dx-sigma}
\Vert {\bf{W}}_{\alpha, p}(f d\sigma)\Vert_{L^r(\mathbb{R}^{n})} \leq C \,  \Vert f \Vert^{\frac{1}{p-1}}_{L^{s} (\mathbb{R}^{n},\, d\sigma)}, \quad \forall f \in L^{s} (\mathbb{R}^{n}, d \sigma),
\end{equation}
where $r>p-1$, $s>1$, and $0<\alpha  p< n$.

It is well-known (see \cite{HJ}, \cite{JPW}) that
\begin{equation}\label{mw}
\Vert {\bf{W}}_{\alpha, p}(f d\sigma)\Vert_{L^r(\mathbb{R}^{n})}  \approx  \Vert {\bf{I}}_{\alpha p}(f d\sigma)\Vert^{\frac{1}{p-1}}_{L^{\frac{r}{p-1}}(\mathbb{R}^{n})},
\end{equation}
with the constants of equivalence that depend on $\alpha, p, q, r, n$,
but not on $f$ and $\sigma$.

Hence, for $r>p-1$ and $s>1$,
 \eqref{winq:wolff-dx-sigma} is   equivalent by duality to
\begin{equation}\label{winq:wolff-dx-sigma-dual}
\Vert  {\bf{I}}_{\alpha p} (g dx)\Vert_{L^{s'}(\mathbb{R}^{n},\,d\sigma)}
\leq c \, C^{p-1}\, \Vert g \Vert_{L^{\frac{r}{r-p+1}} (\mathbb{R}^{n})}, \quad \forall g \in L^{\frac{r}{r-p+1}} (\mathbb{R}^{n}),
\end{equation}
where $(\frac{r}{p-1})'=\frac{r}{r-p+1}$, $s'=\frac{s}{s-1}$ are the dual exponents, and $c=c(\alpha, p, r, s, n)$
is a positive constant.

The next lemma gives a sufficient condition for the validity of \eqref{winq:wolff-dx-sigma}.

\begin{Lem}\label{lem1-wolff}
Let $1< p< \infty$, $0<\alpha< \frac{n}{p}$, $0<q<p-1$, and let  $\sigma \in \mathcal{M}^{+}(\mathbb{R}^{n})$. Suppose that $\frac{n(p-1)}{n-\alpha p}<r<\infty$
and
\begin{equation}\label{cond:dx-wolff}
{\bf{W}}_{\alpha, p} \sigma \in L^{\frac{r(p-1)}{p-1-q}}(\mathbb{R}^n).
\end{equation}
Then \eqref{winq:wolff-dx-sigma} is valid with $s= \frac{r(n-\alpha p)- (p-1) n}{nq}+1$.
Moreover, there exists a positive constant $c=c(\alpha, p, q, r, n)$ such that
\begin{equation}\label{winq:dx-wolff-2}
\Vert {\bf{W}}_{\alpha, p}(f d\sigma)\Vert_{L^r(\mathbb{R}^{n})}
\leq c  \Vert {\bf{W}}_{\alpha, p}\sigma \Vert^{\frac{1}{s'}}_{L^{\frac{r(p-1)}{p-1-q}} (\mathbb{R}^{n})}
\Vert f \Vert^{\frac{1}{p-1}}_{L^{s} (\mathbb{R}^{n}, d \sigma)},
\end{equation}
for all $ f \in L^{s} (\mathbb{R}^{n}, d \sigma)$.
\end{Lem}

\begin{proof}
As was shown above,  \eqref{winq:dx-wolff-2} is equivalent to
\begin{equation}\label{dual-2}
\Vert {\bf{I}}_{\alpha p}(g\,dx)\Vert_{L^{s'}(\mathbb{R}^{n},\,d\sigma)}
\leq c  \Vert {\bf{W}}_{\alpha, p}\sigma \Vert^{\frac{p-1}{s'}}_{L^{\frac{r(p-1)}{p-1-q}} (\mathbb{R}^{n})}
\Vert g\Vert_{L^{\frac{r}{r-p+1}} (\mathbb{R}^{n})},
\end{equation}
for all $g \in L^{\frac{r}{r-p+1}} (\mathbb{R}^n)$, where
without loss of generality we may assume that $g \ge 0$.

To prove the preceding inequality, we
apply the iterated inequality \eqref{iterated1-wolff} in the special case $p=2$, to
the Riesz potential  ${\bf{I}}_{\alpha p} (g \, d x)$,
and $t=s'$, so that
\[
\Big( {\bf{I}}_{\alpha p} (g \, d x) \Big)^{s'} \le c \, {\bf{I}}_{\alpha p} \left [ \Big( {\bf{I}}_{\alpha p} (g \, d x) \Big)^{s'-1} \, g\, dx \right].
\]
Hence, by the preceding inequality and Fubini's theorem,
\[
\begin{split}
\int_{\mathbb{R}^{n}} \Big( {\bf{I}}_{\alpha p} (g \, d x) \Big)^{s'}\,d\sigma
& \leq c \int_{\mathbb{R}^{n}}  {\bf{I}}_{\alpha p} \Big[ \Big( {\bf{I}}_{\alpha p} (g \, d x) \Big)^{s'-1} \, g\, dx
  \Big] d \sigma
  \\
&= c \int_{\mathbb{R}^{n}} \Big({\bf{I}}_{\alpha p} (g \, d x) \Big)^{s'-1} ({\bf{I}}_{\alpha p} \sigma) \, g\,dx.
\end{split}
\]
By H\"older's inequality with the three exponents
$\frac{r}{q}$,
$\frac{r}{p-1-q}>1$ and
$\frac{r}{r-p+1}$, we estimate
\[
\begin{split}
 \int_{\mathbb{R}^{n}} \Big({\bf{I}}_{\alpha p} (g \, d x) \Big)^{s'-1} ({\bf{I}}_{\alpha p} \sigma) \, g\,dx
& \le \big\| {\bf{I}}_{\alpha p}(g dx) \big\|^{s'-1}_{L^{\frac{r(s'-1)}{q}}(\mathbb{R}^n)}
\\
& \times\Vert {\bf{I}}_{\alpha p} \sigma\Vert_{L^\frac{r}{p-1-q}(\mathbb{R}^n)}  \Vert g \Vert_{L^\frac{r}{r-p+1}(\mathbb{R}^n)}.
\end{split}
\]
We observe that
$\frac{r(s'-1)}{q} = \frac{rn}{r(n-\alpha p)- (p-1) n}>1$ and
\[
\tfrac{q}{r(s'-1)} = \tfrac{r(n-\alpha p) -(p-1)n}{rn} = \tfrac{r-p+1}{r} - \tfrac{\alpha p }{n}.
\]
Hence, by
 Sobolev's inequality for the Riesz potential  ${\bf{I}}_{\alpha p}$,
\[
\big\| {\bf{I}}_{\alpha p}(g dx) \big\|_{L^{\frac{r(s'-1)}{q}}(\mathbb{R}^n)}
\leq c \big\| g \big\|_{L^{\frac{r}{r-p+1}}(\mathbb{R}^n)}.
\]
Combining the above estimates proves
\begin{equation*}\label{inq:weight-wolff-dual-2}
\Vert {\bf{I}}_{\alpha p}(g\,dx)\Vert^{{s'}}_{L^{s'}(\mathbb{R}^{n},\,d\sigma)}
\leq c  \Vert {\bf{I}}_{\alpha p}\sigma \Vert_{L^{\frac{r}{p-1-q}} (\mathbb{R}^{n})}
\Vert g\Vert^{s'}_{L^{\frac{r}{r-p+1}} (\mathbb{R}^{n})}.
\end{equation*}
Notice that  by \eqref{mw} with $f\equiv 1$ and $\frac{r(p-1)}{p-1-q}$ in place of $r$, we have
\begin{equation*}
 \Vert {\bf{I}}_{\alpha p}\sigma \Vert_{L^{\frac{r}{p-1-q}} (\mathbb{R}^{n})}
 \approx \Vert {\bf{W}}_{\alpha, p}\sigma \Vert^{p-1}_{L^{\frac{r(p-1)}{p-1-q}} (\mathbb{R}^{n})}.
 \end{equation*}
This completes the proof of  \eqref{dual-2}, and hence
\eqref{winq:dx-wolff-2}.
\end{proof}

The following lemma provides a crucial link between conditions
\eqref{cond:sigma-gamma} and \eqref{cond:dx-wolff}.
In particular, it shows that \eqref{cond:sigma-gamma} yields both weighted norm inequalities \eqref{winq:wolff-sigma-sigma} and \eqref{winq:wolff-dx-sigma}
for suitable $r$ and $s$, in light of Lemma \ref{lem:wolff} and Lemma \ref{lem1-wolff}.

\begin{Lem}\label{lemma-link-p}
Let $1< p< \infty$, $0<\alpha< \frac{n}{p}$, $0<q<p-1$, and let  $\sigma \in \mathcal{M}^{+}(\mathbb{R}^{n})$.
If $\frac{n(p-1)}{n-\alpha p}<r<\infty$, then
\eqref{cond:sigma-gamma} with $\gamma = \frac{r(n-\alpha p)-(p-1)n}{n}$
implies \eqref{cond:dx-wolff}.
\end{Lem}

\begin{proof}
For  $k \in \mathbb{N}$, let $\sigma_{k}:= \chi_{\Omega_{k}}\sigma$, where $\chi_{\Omega_{k}}$ is the characteristic
function of the set
\[
\Omega_{k} =\{ x \in \mathbb{R}^{n}: \,\, {\bf{W}}_{\alpha, p}\sigma(x) \leq k
\}\cap B(0, k).
\]
Then ${\bf{W}}_{\alpha, p}\sigma_k (x)\le k$ on
the support of $\sigma_k$, and  by the weak maximum principle
\eqref{wolff-max}, we have
${\bf{W}}_{\alpha, p}\sigma_k (x)\le  c(\alpha, p, n) \, k$ for all $x \in \mathbb{R}^n$.

Moreover, if $|x|\ge 2 k$, then $\sigma_k(B(x, \rho))=0$ for $0<\rho<\frac{|x|}{2}$, and
hence, for a positive constant $c=c(\alpha, p, n)$,
\[
{\bf{W}}_{\alpha, p}\sigma_k (x)\le \int_{\frac{|x|}{2}}^\infty \left[\frac{\sigma_k(B(x, \rho))}
{\rho^{n-\alpha p}}
\right]^{\frac{1}{p-1}} \frac{d \rho}{\rho}\le  \, c \,
\sigma(B(0, k))^{\frac{1}{p-1}}
|x|^{-\frac{n-\alpha p}{p-1}}.
\]
Thus, $
{\bf{W}}_{\alpha, p}\sigma_k (x) \leq C (|x|+1)^{-\frac{n-\alpha p}{p-1}}$, for all $x \in \mathbb{R}^n$.
Since $r>\frac{n(p-1)}{n-\alpha p}$ and $0<q<p-1$, it follows that
\[
\int_{\mathbb{R}^n} ({\bf{W}}_{\alpha, p}\sigma_{k})^{\frac{r(p-1)}{p-1-q}} dx < + \infty.
\]

Applying the iterated inequality \eqref{iterated1-wolff} with $t=\frac{p-1}{p-1-q}>1$, along with the
weighted norm inequality \eqref{winq:dx-wolff-2} with
$f=({\bf{W}}_{\alpha, p}\sigma_{k})^{\frac{q(p-1)}{p-1-q}}$, we deduce
\begin{align*}
& \int_{\mathbb{R}^n} ({\bf{W}}_{\alpha, p}\sigma_{k}) ^{\frac{r(p-1)}{p-1-q}}\,dx
\leq c \int_{\mathbb{R}^n} \left[ {\bf{W}}_{\alpha, p} (f d \sigma_{k}) \right]^r dx\\
&\leq  c \big\| {\bf{W}}_{\alpha, p}\sigma_{k} \big\|^{\frac{r}{s'}}_{L^{\frac{r(p-1)}{p-1-q}} (\mathbb{R}^n)}
             \big\| f \big\|^{\frac{r}{p-1}}_{L^{s} (\mathbb{R}^n,\,d\sigma_{k})} \\
&\leq c \Big[ \int_{\mathbb{R}^n} ({\bf{W}}_{\alpha, p}\sigma_{k})^{\frac{r(p-1)}{p-1-q}}\,dx \Big]^{\frac{p-1-q}{(p-1)s'}}
\Big[ \int_{\mathbb{R}^n} (\mathbf{W}_{\alpha, p} \sigma)^{\frac{qs(p-1)}{p-1-q}} d \sigma\Big]^{\frac{r}{s(p-1)}},
\end{align*}
where $s= \frac{r(n-\alpha p)- (p-1) n}{nq}+1$. In other words,
\[
\Big[ \int_{\mathbb{R}^n} ({\bf{W}}_{\alpha, p}\sigma_{k})^{\frac{r(p-1)}{p-1-q}}\,dx \Big]^{1-\frac{p-1-q}{s'(p-1)}}
\leq c \Big[ \int_{\mathbb{R}^{n}} ({\bf{W}}_{\alpha, p} \sigma)^{\frac{qs(p-1)}{p-1-q}} d \sigma\Big]^{\frac{r}{s(p-1)}},
\]
where the right-hand side integral is finite by assumption \eqref{cond:sigma-gamma},  since
$\frac{qs(p-1)}{p-1-q}=\frac{(\gamma+q)(p-1)}{p-1-q}$. Notice that $1-\frac{p-1-q}{(p-1)s'}>0$.
Passing to the limit as $k \rightarrow \infty$ in the preceding estimate yields \eqref{cond:dx-wolff}.
\end{proof}

We now prove Theorem \ref{thm:main-p-laplacian} and Theorem \ref{thm:main2}.

\begin{proof}[Proof of Theorem \ref{thm:main-p-laplacian}] It is well known
\cite{HKM} that,
for $p \ge n$,  there is no positive $p$-superharmonic functions on $\mathbb{R}^n$. Hence,
\eqref{eq:p-laplacian}  has no nontrivial solutions.

In the case $1<p<n$,  suppose \eqref{eq:p-laplacian}  has a nontrivial solution
$u \in L^r(\mathbb{R}^n)$, $r>0$. Let $d \omega =u^q d\sigma$,  and
pick $R>0$ so that $\omega(B(0, R))>0$. Then by Theorem \ref{pointwise_est_p-superharmonic}, we have, for $|x|> 2R$,
\[
u(x) \ge C \,  {\bf{W}}_{1, p}\omega (x)\ge C \, \omega(B(0, R))^{\frac{1}{p-1}} \, |x|^{-\frac{n-p}{p-1}}.
\]
Hence,  $u \in L^r(\mathbb{R}^n)$ yields
$r>\frac{n(p-1)}{n-p}$.

Suppose now that $1<p<n$, $r>\frac{n(p-1)}{n-p}$, and
\eqref{cond:dsigma-wolff} is valid. In view of Theorem \ref{thm:wolff-existence} with $\alpha=1$, there exists a positive solution
$v \in L^{\gamma +q}(\Omega, d\sigma)$, where $\gamma = \frac{r(n-p)-n(p-1)}{n}>0$, to the integral equation
\begin{equation}\label{eq:int-wolff-1}
v = {\bf{W}}_{1,p}(v^{q}d\sigma) \quad \text{in} \;\; \mathbb{R}^n.
\end{equation}
Moreover, by Lemma \ref{lem1-wolff} and Lemma \ref{lemma-link-p}, the weighted
norm inequality \eqref{winq:wolff-dx-sigma} holds with $s=\frac{\gamma+q}{q}$.
Letting $f=v^{q} \in L^{\frac{\gamma + q}{q}}(\mathbb{R}^n)$ in \eqref{winq:wolff-dx-sigma} yields
\[
\begin{split}
\| v \|_{L^{r}(\mathbb{R}^n)} &= \big\| {\bf{W}}_{1, p}(v^{q}d\sigma) \big\|_{L^{r}(\mathbb{R}^n)}
\leq c \| v^q \|^{\frac{1}{p-1}}_{L^{\frac{\gamma + q}{q}}(\mathbb{R}^n, d \sigma)} \\
&= c \| v \|_{L^{\gamma +q}(\mathbb{R}^n, d \sigma)}^{\frac{q}{p-1}}
< +\infty.
\end{split}
\]
This shows that the solution $v \in L^{r}(\mathbb{R}^n)$ as well.

Using an iterative procedure  as in the proof of \cite[Theorem 1.1]{CV2},
based on  Theorem \ref{weak_cont_p-Laplacian} and Theorem \ref{pointwise_est_p-superharmonic},
we construct a (minimal) positive $p$-superharmonic solution $u$
to \eqref{eq:p-laplacian} so that $u(x) \le C \, v(x)$. Consequently,
$u \in L^{r}(\mathbb{R}^n)$.
\end{proof}

\begin{proof}[Proof of Theorem \ref{thm:main2}]
Since a superharmonic solution $u$ to \eqref{eq:frac-laplacian} is understood in the sense that
\begin{equation}\label{eq:int-riesz}
u = {\bf{I}}_{2\alpha}(u^{q}d\sigma) \quad \text{in} \;\; \mathbb{R}^n,
\end{equation}
then by using an argument similar to that  in the proof of Theorem \ref{thm:main-p-laplacian} with
${\bf{I}}_{2\alpha} = {\bf{W}}_{\alpha, 2}$ in place of ${\bf{W}}_{1,p}$, we see that condition \eqref{cond:dsigma-riesz} yields the existence of a (minimal) positive solution $u \in L^{r}(\mathbb{R}^n)$
to \eqref{eq:frac-laplacian} for $r>\frac{n}{n- 2 \alpha}$. If $0<r \le \frac{n}{n- 2 \alpha}$, then there is no nontrivial solution $u \in L^r(\mathbb{R}^n)$, since
$u(x)={\bf{I}}_{2\alpha}(u^{q}d\sigma)(x)\ge C \, |x|^{2 \alpha-n}$ for $|x|$ large, as in   the
$p$-Laplace case considered above.
\end{proof}

The next proposition is used to deduce Corollaries \ref{cor:main1} and \ref{cor:main3}.

\begin{Prop}\label{prop:wolff}
Let $1<p<n$, $0<\alpha<\frac{n}{p}$ and $\beta>0$.
If $\sigma \in L^{s}(\mathbb{R}^n)$ is a nonnegative function, where $s=\frac{n(\beta+p-1)}{n(p-1)+p\alpha\beta}$, then
\begin{equation}\label{cond:sigma-gen}
{\bf{W}}_{\alpha, p}\sigma \in L^{\beta}(\mathbb{R}^{n}, d\sigma).
\end{equation}
\end{Prop}

\begin{proof}
Since $s>1$, by H\"{o}lder's inequality, we have
\[
\int_{\mathbb{R}^n} ({\bf{W}}_{\alpha, p}\sigma)^{\beta}\,d\sigma
\leq \big\| {\bf{W}}_{\alpha, p}\sigma \big\|_{L^{\beta s'}(\mathbb{R}^n)}^{\beta}
\| \sigma \|_{L^{s}(\mathbb{R}^n)}.
\]
Notice that
$
\frac{p-1}{\beta s'} = \frac{1}{s} - \frac{\alpha p}{n}>0.
$
Estimating again   the norm of ${\bf{W}}_{\alpha, p} \sigma$ in terms of the norm of
the Riesz potential ${\bf{I}}_{\alpha p} \sigma$ by \eqref{mw}, along with Sobolev's inequality
for  ${\bf{I}}_{\alpha p} \sigma$, we  obtain
\[
\big\| {\bf{W}}_{\alpha, p}\sigma \big\|_{L^{\beta s'}(\mathbb{R}^n)}
\leq c \big\| {\bf{I}}_{\alpha p}\sigma \big\|_{L^{\frac{\beta s'}{p-1}}(\mathbb{R}^n)}^{\frac{1}{p-1}}
\leq c \| \sigma \|_{L^{s}(\mathbb{R}^n)}^{\frac{1}{p-1}}.
\]
Combining the preceding estimates yields \eqref{cond:sigma-gen}.
\end{proof}

Letting  $\beta=\frac{(\gamma+q)(p-1)}{p-1-q}$, where  $0<q<p-1$, and
$\gamma=\frac{r(n-\alpha p)-n(p-1)}{n}$, where $r>\frac{n(p-1)}{n-\alpha p}$,
we obtain $\beta>0$ and $s=\frac{nr}{n(p-1-q)+ \alpha p r}>1$ in Proposition \ref{prop:wolff}. Notice that condition \eqref{cond:sigma-gen} coincides with \eqref{cond:sigma-gamma}.
Consequently, $\sigma \in L^{s}(\mathbb{R}^n)$ is sufficient for  \eqref{cond:sigma-gamma}
by Proposition \ref{prop:wolff}.

In the case $\alpha=1$, we obtain that
\eqref{cond:sigma-suff} implies \eqref{cond:dsigma-wolff}.
Thus,  Corollary \ref{cor:main1} follows from Theorem \ref{thm:main-p-laplacian}.

Similarly, applying Proposition \ref{prop:wolff} with $p=2$, we see that
\eqref{cond:sigma-suff2} implies \eqref{cond:dsigma-riesz}.
Hence, Corollary \ref{cor:main3} holds in view of Theorem \ref{thm:main2}.

\section{Solutions to problem \eqref{eq:laplacian}}
As mentioned above, our method of proof of Theorem \ref{thm:main3}
is based on both the
weighted norm inequality \eqref{weighted_norm_ineq}, and a similar inequality
with Lebesgue measure on the left-hand side:
\begin{equation}\label{inq:weighted}
\Vert \mathbf{G} (f d\sigma)\Vert_{L^r(\Omega)} \leq c \Vert f \Vert_{L^{s} (\Omega,\, d\sigma)}, \quad \forall f \in L^{s} (\Omega, d \sigma),
\end{equation}
for $1<r,s<\infty$. Notice that by duality, \eqref{inq:weighted} is equivalent to
\begin{equation}\label{inq:weight-dual}
\Vert \mathbf{G} (g dx)\Vert_{L^{s'}(\Omega,\,d\sigma)}
\leq c \Vert g \Vert_{L^{r'} (\Omega)}, \quad \forall g \in L^{r'} (\Omega),
\end{equation}
where $r'=\frac{r}{r-1}$ and $s'=\frac{s}{s-1}$.

We observe that such inequalities were characterized when ${\bf{G}}$ is
the Riesz potential on $\mathbb{R}^n$, in terms of Riesz capacities, by Maz'ya and Netrusov,
and in terms of nonlinear  potentials,  by Cascante, Ortega, and Verbitsky (see \cite{COV}, \cite{Maz}).
Our first lemma gives sufficient conditions for the validity of
 \eqref{inq:weighted}.

\begin{Lem}\label{lemma1}
Let $0<q<1$ and $\sigma \in \mathcal{M}^{+}(\Omega)$. Let $G$ be a positive lower semicontinuous kernel on
$\Omega \times \Omega$, which satisfies the WMP and \eqref{cond:Riesz0} for some
$0<\alpha<\frac{n}{2}$. Suppose that $\frac{n}{n-2\alpha}<r<\infty$
and that condition \eqref{cond:dx-green} holds.
Then \eqref{inq:weighted} is valid with $s= \frac{r(n-2\alpha) - n(1-q)}{nq}$.
Moreover,
\begin{equation}\label{inq:weighted-2}
\Vert \mathbf{G}(f d\sigma)\Vert_{L^r(\Omega)}
\leq c  \Vert \mathbf{G} \sigma \Vert^{\frac{1}{s'}}_{L^{\frac{r}{1-q}} (\Omega)}
\Vert f \Vert_{L^{s} (\Omega, d \sigma)}, \quad \forall f \in L^{s} (\Omega, d \sigma),
\end{equation}
where $c$ is a positive constant independent of $f$ and $\sigma$.
\end{Lem}

\begin{proof}
Observe that the hypotheses ensure
\[
1<s<\infty, \quad \tfrac{(s'-1)r}{q}>1 \quad  \text{and} \quad
\tfrac{1}{r'} - \tfrac{q}{(s'-1)r} = \tfrac{2\alpha}{n}.
\]
By duality, \eqref{inq:weighted-2}  is equivalent to
\begin{equation}\label{inq:weight-dual-2}
\Vert \mathbf{G} (g\,dx)\Vert_{L^{s'}(\Omega,\,d\sigma)}
\leq c  \Vert \mathbf{G} \sigma \Vert^{\frac{1}{s'}}_{L^{\frac{r}{1-q}} (\Omega)}
\Vert g \Vert_{L^{r'} (\Omega)}, \quad \forall g \in L^{r'} (\Omega),
\end{equation}
where $c$ is a positive constant independent of $g$  and $\sigma$.
Without loss of generality, suppose $g \in L^{r'} (\Omega)$ with $g \ge 0$.
Applying iterated inequality \eqref{iterated1}
with $t=s'$, followed by Fubini's theorem and H\"older's inequality with exponents $\frac{r}{1-q}$ and
$\frac{r}{r-1+q}$, we deduce
\[
\begin{split}
 \int_\Omega \Big(\mathbf{G} (g \, d x) \Big)^{s'}\,d\sigma
&  \leq c \int_\Omega \mathbf{G} \Big[ \Big(\mathbf{G} (g \, d x) \Big)^{s'-1} \, g dx
  \Big] d \sigma \\
&= c \int_\Omega \Big(\mathbf{G} (g \, d x) \Big)^{s'-1} (\mathbf{G} \sigma) \, g dx.
\end{split}
\]
Using  H\"older's inequality with the exponents $\frac{r}{q}$,
$\frac{r}{1-q}$ and $r'$, we estimate
\[
\begin{split}
\int_\Omega \Big(\mathbf{G} (g \, d x) \Big)^{s'-1} (\mathbf{G} \sigma) \, g dx
& \leq \Vert \mathbf{G} (g \, d x) \Vert^{s'-1}_{L^{\frac{(s'-1)r}{q}} (\Omega)} \\ &
 \times \Vert \mathbf{G} \sigma \Vert_{L^\frac{r}{1-q}(\Omega)}
\Vert g \Vert_{L^{r'}(\Omega)} .
\end{split}
\]
Denote by $\tilde{g}$ the zero extension of $g$ to $\mathbb{R}^n$. By assumption
\eqref{cond:Riesz0} and Sobolev's inequality for Riesz potentials, we have
\[
\| {\bf{G}}(g\,dx) \|_{L^{\frac{(s'-1)r}{q}}(\Omega)}
\leq c \| {\bf{I}}_{2\alpha}\tilde{g}\|_{L^{\frac{(s'-1)r}{q}}(\mathbb{R}^n)}
\leq  c \| \tilde{g}\|_{L^{r'}(\mathbb{R}^n)}
= c \| g \|_{L^{r'}(\Omega)}.
\]
Combining the preceding inequalities yields \eqref{inq:weight-dual-2} as desired.
\end{proof}

The following lemma provides an important link between conditions \eqref{cond:dx-green} and \eqref{cond:dsigma-green}.
In particular, it shows that \eqref{cond:dsigma-green} yields both weighted norm inequalities \eqref{weighted_norm_ineq} and \eqref{inq:weighted} for suitable $r$ and $s$, in light of Theorem \ref{thm_Ver17} and Lemma \ref{lemma1}.

\begin{Lem}\label{lemma2}  Let $0<q<1$ and $\sigma \in \mathcal{M}^{+}(\Omega)$ with $\sigma \not\equiv 0$. Suppose $G$ is a positive lower semicontinuous kernel on
$\Omega \times \Omega$, which satisfies the WMP and \eqref{cond:Riesz0} for some
$0<\alpha<\frac{n}{2}$. If $\frac{n}{n-2\alpha}<r<\infty$, then
\eqref{cond:dsigma-green} with $\gamma=\frac{r(n-2 \alpha)-n}{n}$
implies \eqref{cond:dx-green}.
\end{Lem}

\begin{proof}
Let $\lbrace \Omega_{k}\rbrace_{k \in \mathbb{N}}$ be an  increasing exhaustive sequence of relatively compact open sets in $\Omega $ such that ${\bf{G}}\sigma (x)\leq k$ for $x\in \Omega_k$. Set $\sigma_{k}:= \chi_{\Omega_{k}}\sigma$. By the WMP, ${\bf{G}}\sigma_k \leq h k$ for all
$x \in \Omega$. Moreover, if $\Omega_{k}\subset B(0, R)$, and $|x|>2R$, then
${\bf{G}}\sigma_k (x)\leq C \, |x|^{2 \alpha-n}$ by \eqref{cond:Riesz0}. Consequently, ${\bf{G}}\sigma_k (x)
\leq C \, (1+|x|)^{2 \alpha-n}$ for all $x\in\Omega$. Since $r>\frac{n}{n-2 \alpha}$ and
$0<q<1$, we have
\[
\int_{\Omega} ({\bf{G}}\sigma_{k})^{\frac{r}{1-q}} dx < + \infty.
\]
Applying the iterated inequality \eqref{iterated1} with $t=\frac{1}{1-q}$, followed by
the weighted norm inequality \eqref{inq:weighted-2} with $f=({\bf{G}}\sigma_{k})^{\frac{q}{1-q}}$, we estimate
\begin{align*}
& \int_\Omega ({\bf{G}}\sigma_{k})^{\frac{r}{1-q}}\,dx
\leq c \int_\Omega \left[ \mathbf{G} (f d \sigma_{k}) \right]^r dx\\
&\leq  c \big\| {\bf{G}}\sigma_{k} \big\|^{\frac{r}{s'}}_{L^{\frac{r}{1-q}} (\Omega)}
             \big\| f \big\|^r_{L^{s} (\Omega,\,d\sigma_{k})}\\
&\leq c \Big[ \int_{\Omega} ({\bf{G}}\sigma_{k})^{\frac{r}{1-q}} dx \Big]^{\frac{1-q}{s'}}
\Big[ \int_\Omega (\mathbf{G} \sigma)^{\frac{qs}{1-q}} d \sigma\Big]^{\frac{r}{s}},
\end{align*}
where $s= \frac{r(n-2\alpha) - n(1-q)}{nq}$. In other words,
\[
\Big[ \int_{\Omega} ({\bf{G}}\sigma_{k})^{\frac{r}{1-q}} dx\Big]^{1-\frac{1-q}{s'}}
\leq c \Big[ \int_\Omega (\mathbf{G} \sigma)^{\frac{qs}{1-q}} d \sigma\Big]^{\frac{r}{s}}
< +\infty,
\]
where the right-hand side integral is finite by assumption \eqref{cond:dsigma-green} with $\gamma=\frac{r(n-2 \alpha)-n}{n}$,  since
$qs=\gamma+q$. Here obviously $1-\frac{1-q}{s'}>0$.
Passing to the limit as  $k \rightarrow \infty$
in the preceding estimate yields \eqref{cond:dx-green}.
\end{proof}
We are now ready to prove Theorem \ref{thm:main3}.
\begin{proof}[Proof of Theorem \ref{thm:main3}]

Suppose \eqref{cond:dsigma-green} holds.
In view of Theorem \ref{thm_Ver17}, there exists a positive solution
$u \in L^{\gamma +q}(\Omega, d\sigma)$ to the integral  equation
\begin{equation}\label{eq:int-green}
u = {\bf{G}}(u^{q}d\sigma) \quad \text{in} \;\; \mathbb{R}^n.
\end{equation}
Further, in light of Lemma \ref{lemma1} and Lemma \ref{lemma2},
the weighted norm inequality \eqref{inq:weighted} holds with $s=\frac{\gamma + q}{q}$.
Hence, applying \eqref{inq:weighted} with $f =u^{q} \in L^{\frac{\gamma+q}{q}}(\Omega, d\sigma)$ yields
\[
\| u \|_{L^{r}(\Omega)}
\leq \| {\bf{G}}(u^{q}d\sigma) \|_{L^{r}(\Omega)}
\leq c \| u^{q} \|_{L^{\frac{\gamma + q}{q}}(\Omega)}
= c \| u \|_{L^{\gamma +q}(\Omega)}^{q}
< + \infty.
\]
This shows that  $u\in L^{r}(\Omega)$ is a positive superharmonic solution to \eqref{eq:laplacian}.\end{proof}

To deduce Corollary \ref{cor:main3}, we need the following proposition
(\cite[Proposition 4.8]{SV2}).

\begin{Prop}\label{prop3} Let $\Omega \subseteq \mathbb{R}^n$.
Let $G$ be a positive kernel on $\Omega\times \Omega$ which satisfies \eqref{cond:Riesz0} with $0<2 \alpha<n$.
 If $\sigma \in L^{s}(\Omega)$ is a positive function, where
$
s=\frac{n(\beta+1)}{n+2\alpha\beta},
$ and $\beta>0$,
then
\begin{equation}\label{cond_general}
{\bf{G}}\sigma \in L^{\beta}(\Omega, d\sigma).
\end{equation}
\end{Prop}

Proposition \ref{prop3} with $\alpha=1$ and
$\beta=\frac{r(n-2)-n(1-q)}{n(1-q)}=\frac{\gamma+q}{1-q}$
shows that  \eqref{cond:sigma-suff3} implies \eqref{cond:dsigma-green}, and
consequently  Theorem \ref{thm:main3} yields Corollary \ref{cor:main3}.


\section{Introduction}
\label{intrd}

Glottal closures are the instants of significant excitation to the vocal tract system occurs during closure of vocal folds for every cycle of vocal fold vibration during phonation. For healthy vocal folds, the glottal closures can be approximated with a sequence of strong impulses~\cite{rabiner1978digital} due to abrupt closure of vocal folds, where as a disorder vocal fold produces a weak or smeared excitation due to incomplete closure of vocal folds~\cite{yamauchi2016quantification}. Hence, the methods developed to detect GCI from speech signal recorded from healthy vocal folds cannot be directly applied on the pathological speech. Therefore, we can find significantly less attempts to directly extract the GCI from pathological speech. An alternative to the acoustic speech signal is the laryngeal (Electroglottograph or EGG) signal, which is a low frequency signal recorded by measuring the impedance across the larynx by passing a weak electrical current with the help of two electrodes~\cite{childers1985critical}. The EGG signal thus captures the vocal folds activity free from vocal tract resonances which is prominent in raw speech. The negative peaks in the derivative of EGG co-insides with the GCI hence, we can find significantly many works based on EGG for detecting GCI from healthy and disordered vocal folds~\cite{deshpande2018effective, awan2013effect, thomas2009sigma, rothenberg1988monitoring, mittal2014study, mecke2012comparing, chen2013speech}. It is not always possible to record the EGG from the larynx, requires a dedicated device and requires at least basic skills to record EGG from the device. Where as acoustic speech can be easily recorded even with a microphone in a hand held device such as a smart phone, which can be processed in device to give the preliminarily assessment before consulting a doctor or the processed data can be sent to a pathologist for further analysis. Hence, there is a great need for GCI detection directly from speech signal.

\textbf{Studies based on Modal and Pathological Speech.} We can find many algorithms for GCI detection from modal or normal speech. Most of the available GCI detection methods relies on a representative signal from speech which emphasizes the locations of glottal closure instants. Further, the GCI are detected from the peaks of representative signal either by eliminating spurious instants or by picking genuine GCI using hand crafted heuristics. Most of the methods prefer linear prediction residual (LPR) as a representative signal~\cite{naylor2007estimation, thomas2012estimation, prathosh2016cumulative, koutrouvelis2016fast}, which is a correlate of glottal source signal based on source filter model of speech production~\cite{fant2012acoustic}. Other methods exploit the properties of impulsive nature of excitation~\cite{murty2008epoch, drugman2009glottal, d2011glottal, khanagha2014detection} to derive the representative signal. Recently, \textit{Rachel et.al}~\cite{rachel2017estimation} and \textit{Kebin Wu et. al}~\cite{wu2017gmat} showed their GCI detection method tuned for modal speech as an application to pathological speech on a limited data, which consists of one type of disorder: dysphonia. It should be noted that afore mentioned methods are developed for modal/normal speech and heavily depends on the choice of representation signal, model assumptions to extract the GCI. Also, they depend heavily on signal processing pipe lines requires manual tunning of parameters and hand crafted heuristics specific to dataset to reliably detect the GCI from speech. Recently, classification based data driven model is proposed for GCI detection from modal speech~\cite{matouvsek2017classification, matouvsek2018glottal}. It should be noted that the model is trained with hand crafted features extracted from speech signal and voting classifiers to detect the GCI. In this paper, we propose a deep CNN model, which requires no manual parameter tuning, trained with both raw acoustic and LPR to predict the GCI from pathological speech.

\section{Vocal Disorder Speech and EGG Dataset}
\label{dataset}

For carrying out this study, we collected simultaneous recordings of both EGG and speech signals from the patients who had pathological disorder in vocal folds from B. C. Roy Technology Hospital, Indian Institute of Technology Kharagpur, India. We have collected data from 78 patients registered for diagnosis and treatment from the same hospital. There were 45 male and 33 female patients with the age group of 22-72 years and 12-63 years respectively. The patients suffering from vocal fold disorders are categorized into six types namely vocal nodule (\textbf{N}), vocal polyp (\textbf{P}), laryngitis (\textbf{L}), thickened vocal folds (\textbf{T}), cancer (\textbf{C}) and paralysis (\textbf{PV}) vocal folds. The speech samples are captured using a high quality microphone. The EGG signals are recorded using clinical grade Electroglottograph (procured from TechCadenza, India) device. The acoustic speech and EGG of the sustained vowels /a/, /e/, /o/ are recorded each of them for three times from the participating patients. Sample speech and EGG signals for cancer (\textbf{C}), nodule (\textbf{N}), laryngitis (\textbf{L}) and  thickened vocal folds (\textbf{T}) are shown in Fig~\ref{fig:speec_egg} (due to page constraints, signals are shown for four disorders).

\begin{figure}[h]
        \centering
		\resizebox{14cm}{5cm}{\includegraphics{speech_egg-eps-converted-to.pdf}}
		\vspace{-0.5cm}
        \caption{Pathological speech and the corresponding EGG for various vocal disorders.}
        \vspace{-0.4cm}
        \label{fig:speec_egg}
\end{figure}

\section{Proposed GCI detection Method}
\label{prop_method}

We propose a classification based GCI detection for pathological speech data by training multi-column (parallel) deep CNN models. The proposed CNN based GCI detection architecture is shown in Fig.~\ref{fig:cnn_model}.

\begin{figure}[h]
        \centering
		\resizebox{14cm}{5.5cm}{\includegraphics{cnn_architecture.pdf}}
		\vspace{-0.5cm}
        \caption{Multi-column CNN based GCI detection (FS = Filter size, FC = fully connected).}
        \vspace{-0.4cm}
        \label{fig:cnn_model}
\end{figure}

\subsection{Feature Representation, Training and Testing Dataset}
\label{feat_rep}
To leverage the advantage of both raw speech and Linear prediction residual (LPR), we trained the deep models on the following input representations. 1) The low-pass filtered pathological speech signal (\textbf{LPF\_S}). Speech is low-passed since high frequency components do not contribute to the GCI. 2) The low-pass filtered LPR (\textbf{LPF\_LPR}) with LP order 12. The LPR is a noisy residual signal contains unwanted high frequency content hence low-pass filtered. 3) The positive clipped low-pass filtered speech (\textbf{PC\_LPF\_S}) and LPR (\textbf{PC\_LPF\_LPR}). Positive clipped since most of the GCI information is present in the negative portion of the signal~\cite{rengaswamy2016robust}. Low-pass filtering is performed with zero phase, order six low-pass Butterworth filter with cut-of-frequency of 1000Hz. The contemporaneous EGG signal recorded along with the speech is used as reference to manually mark the GCI locations on the speech signal. Speech signals are downsampled to 16 KHz and switched to negative polarity before labeling. The negative peaks in the derivative of the EGG signal is taken as reference to place the GCI markers after compensating the delay between the EGG and corresponding Speech signal. The input representation signals are chunked into frames of 16 samples and assigned a label for presence or absence of GCI. The frame of 16 samples around the GCI (captures GCI slope, shape and amplitude) shown in Fig.~\ref{fig:feature_rep} as dotted box plot is labeled as GCI frame and rest as non-GCI frames for training the model. The training and testing data consist of 62 and 18 pathological speech samples respectively. It should be noted that each sample in the dataset is from unique patient and more than one patient can have same type of vocal disorder.

\begin{figure}[h]
        \centering
		\resizebox{14cm}{5.5cm}{\includegraphics{featRep_gciMarkers_features-eps-converted-to.pdf}}
		\vspace{-0.5cm}
        \caption{Input representation, ground truth GCI markers and GCI frame (samples around the GCI).}
        \vspace{-0.4cm}
        \label{fig:feature_rep}
\end{figure}

\subsection{Multi-Column CNN Model}

The proposed multicolumn deep CNN model for GCI detection is shown in Fig.~\ref{fig:cnn_model}. The final model consists of three CNN networks. Each network is trained with a input representation discussed in~\ref{feat_rep}. Further, each CNN network consists of five convolution layers, each convolution layer is followed by a batch normalization. In our model, we have skipped max pooling layers because we want the model to capture the variations of GCI regions due to stochastic nature of input signal. The resulting feature vector from the CNN layer is connected densely to a sigmoid activation function to predict the posterior classification probability for each frame. The binary cross entropy loss function is optimized by ADAM optimizer with learning rate of 0.0001. The model is trained for 30 epochs until there is no change in validation loss. The weights are initialized from He normal distribution~\cite{kubitschek1962normal}.

\subsection{Experiments}

Initially, we trained single column CNN model for each input feature representation discussed in~\ref{feat_rep} separately to evaluate its significance for pathological GCI detection. We denote the models trained with \textbf{LPF\_S}, \textbf{LPF\_LPR}, \textbf{PC\_LPF\_S} and \textbf{PC\_LPF\_LPR} as \textbf{Model 1}, \textbf{Model 2}, \textbf{Model 3} and \textbf{Model 4} respectively. The F1-scores of \textbf{Model 1}, \textbf{Model 2}, \textbf{Model 3} and \textbf{Model 4} are 86.82, 82.94, 86.02, 85.52 respectively. We can observe that models trained with \textbf{LPF\_S}, \textbf{PC\_LPF\_S}, \textbf{PC\_LPF\_LPR} have better F1-score than \textbf{Model 2} trained with \textbf{LPF\_LPR}. A further investigation into the predicted class probabilities from \textbf{Model 1}, \textbf{Model 3} and \textbf{Model 4} revealed that \textbf{Model 3} assigned a little high probability to the non-GCI frames where secondary excitations are prominent, results in false alarms. \textbf{Model 4} assigned very low probability at the low voiced and transition region frames results in high missed rates. Also, \textbf{Model 4} assigns very low probability for the frames with dominant secondary excitation shown in Fig.~\ref{fig:post_prob}, which is good for reducing false alarms. Hence, we trained a joint acoustic-residual model to reap the benefit of both models shown in Fig.~\ref{fig:cnn_model}. In this joint model, we extract the features trained from \textbf{Model 3} and \textbf{Model 4} and append them to a densely connected sigmoid activation function to predict the class probability. Since the model trained with features from low pass filtered speech signal \textbf{LPF\_S} also gave good results, we combine the posterior probabilities of joint acoustic-residual model and model trained with \textbf{LPF\_S} i.e, \textbf{Model 1} in a maximum likelihood sense to predict the final class probability for classification. The frames which achieves class probability greater than or equal to 0.1 are classified as GCI frames. The location of maximum negative peak in the classified frame is considered as glottal closure instant.

\begin{figure}[h]
        \centering
        \label{tab:results}
		\resizebox{12cm}{3.5cm}{\includegraphics{predClassProb_LPFSpeech-eps-converted-to.pdf}}
		\vspace{-0.2cm}
        \caption{Posterior class probabilities predicted from \textbf{Model 3} and \textbf{Model 4}.}
        \vspace{-0.4cm}
        \label{fig:post_prob}
\end{figure}

\section{Evaluation, Results and Summary}
The proposed method is assessed with the reliability and accuracy measures given in~\cite{thomas2009sigma}. Identification Rate (\textbf{IDR}): measures the percentage of GCI for which exactly one GCI is detected. Miss Rate (\textbf{MR}): the percentage of GCI for which no GCI is detected. False Alarm (\textbf{FAR}): the percentage of GCI for which more than one GCI is detected. Identification accuracy (\textbf{IDA}): the standard deviation of the timing error between the detected and the reference GCI location (lower IDA is the better). We compared our method with the popular state-of-the-art GCI detection methods: SEDREAMS~\cite{drugman2009glottal}, DYPSA~\cite{naylor2007estimation} and ZFF~\cite{murty2008epoch}. The evaluation results of the proposed method compared with the other state-of-the-art methods shown in Table~\ref{sample-table}. From Table~\ref{sample-table}, we can observe that the proposed method is significantly better than other methods in detecting GCI with high \textbf{IDR}, low miss rate and false alarms compared to other methods. In summary, we proposed a combined joint acoustic-residual classification based GCI detection method for pathological speech. The evaluation results showed that the proposed method performs significantly better than other methods, which gave the hope to future research on acoustic speech for vocal disorder speech GCI detection and classification.

\begin{table}
  \caption{Comparison of proposed GCI (final model) detection method with other methods.}
  \label{sample-table}
  \centering
  \begin{tabular}{lllll}
    \toprule
    Method     & \textbf{IDR}     & \textbf{MR}     & \textbf{FAR}     & \textbf{IDA}(msec)\\
    \midrule
    Proposed   &89.30         &5.97        &4.70         &0.54     \\
    SEDREAMS   &84.98         &9.06        &6.01         &0.42     \\
    DYPSA      &79.33         &13.42        &7.25         &0.69     \\
    ZFF        &80.88         &10.64        &8.48         &0.56     \\

    \bottomrule
  \end{tabular}
\end{table}

\small

\bibliographystyle{unsrt}


\section{Introduction}
The interaction between light and matter is a central and fundamental problem in quantum optics. The simplest light--matter interaction model is the Jaynes--Cummings (JC) model, proposed by Jaynes and Cummings within the rotation-wave approximation (RWA). The JC model describes the interaction between a two-level atom and a single-mode field~\cite{JC1962}, which is an exactly
solvable model with applications in many fields~\cite{GB1991, RE2011, HW2012, Mb2015, FS2018}. The JC model is experimentally realized
by the squeezed radiation field in a $^{85}$Rb-atom micromaser ~\cite{Rempe1987} and by the collapse and revival of atomic inversion in a $^{138}$Ba-atom microlaser ~\cite{An1994}.

The JC model has been generalized in different ways. For example, the Dicke model (also known as the Tavis--Cummings model) couples the single mode with multiple two-level atoms. The coupling of multi-modes and multi-atoms can be constructed by Yang--Baxter algebra and is integrable through the Bethe ansatz~\cite{Kundu2003}. An extra counter-rotating term in the JC model gives the famous Rabi model, which has been only analytically solved by Braak~\cite{Braak2011} till date. The JC model can also be generalized by including nonlinear couplings between a single atom and the radiation field~\cite{MH2014, FS2014, LL2018, A2012}. An example is the intensity-dependent JC
model, which exhibits significantly different dynamics in the absence and presence of RWA, owing to the dramatically enhanced field-squeezing effect in the latter case~\cite{Ng2000}.
The Buck--Sukumar (BS) model~\cite{BS1981} is a more general JC model considering the interaction of a two-level system with a field under intensity-dependent
coupling, which is still analytically solvable.
In contrast, the BS model with the counter-rotating term (sometimes called the nonlinear Rabi model) is non-analytic. Although the nonlinear term in the intensity-dependent JC or the BS model involves multiphoton interactions that are not physically realizable in current cavity- or circuit-QED setups, they are still of interest to the quantum optics and cold-atom fields. For instance, the cold atoms in optical lattices may provide a means of realizing the nonlinear light--atom coupling appearing in BS-type Hamiltonians, via the engineering of adequate optical lattices~\cite{Celi2017}.

The dynamical nature of the quantum Rabi-like model has been recently discussed. Through time-dependent correlation functions, the dynamical correlation functions of the quantum Rabi model involve the JC scheme in the ultrastrong-coupling and deep strong-coupling regimes~\cite{Wolf2013,Xie2017}. The nonlinear dynamics of trapped-ion models have been proposed for blocking the propagation of quantum information along the Hilbert space of the JC and quantum Rabi models~\cite{Cheng2017}.
Exploiting the parity symmetry of the Rabi model, Hu et al. ~\cite{Hu2017} recursively evolved the
corresponding positive and negative parities using the expansion coefficients in the dynamical equations.

In this paper, quantum spin-boson systems expressed by time-independent Hamiltonians $\hat{H}$ are evolved by a non-perturbation numerical approach. When the stationary eigenvalues of the Hamiltonians are known, the non-perturbation numerical approach recovers the existing results~\cite{HW2012, LW2009} obtained by evolving the known eigenvalues and eigenvectors of $\hat{H}$~\cite{ZX2014, Braak2011} from the initial state. However, when the eigenvalues are unknown, our method can evolve the properties from different initial states. Herein, we apply this method to the BS model and discuss the evolution of the mean photon number, the atomic inversion, and the ground-state properties in the absence and presence of RWA. Ng {\it et al.} proved that the BS Hamiltonian is non-real positive-definite when the strengths of the rotating- and counter-rotating-wave terms are equal~\cite{Ng2000}. Such a system has no physical meaning and no definite time evolution. However, Rodr\'{i}guez-Lara {\it et al.}~\cite{Rodrguez2013} temporally evolved the equivalent nonlinear JC photonic lattice using quantum optics-based methods via the analogy between the transport of single-photon states and propagation of a classical field. We resolve this obvious contradiction through numerical calculations.

The remainder of this paper is organized as follows. Section ~\ref{sec:model} introduces the BS model, and Sec.~\ref{sec:method} presents our non-perturbation numerical method. Sections~\ref{sec:dis1} and ~\ref{sec:dis2} present the BS results with and without RWA, respectively. Conclusions are presented in Sec.~\ref{sec:con}.

\section{The model}~\label{sec:model}
The BS model generalizes the JC model to allow nonlinear couplings~\cite{Ng2000}. The Hamiltonian $(\hbar=1)$ of the BS model is given by
\begin{equation}
\label{H}
\begin{split}
\hat{H} &=\omega_{f}\hat{a}^{\dagger}\hat{a}
+\frac{\omega_{0}}{2}\hat{\sigma}_{z}
+g_{-}(\hat{a}\sqrt{\hat{a}^{\dagger}\hat{a}}\hat{\sigma}_{+}
+\sqrt{\hat{a}^{\dagger}\hat{a}}\hat{a}^{\dagger}\hat{\sigma}_{-})\\
&+g_{+}(\hat{a}\sqrt{\hat{a}^{\dagger}\hat{a}}\hat{\sigma}_{-}
+\sqrt{\hat{a}^{\dagger}\hat{a}}\hat{a}^{\dagger}\hat{\sigma}_{+}),
\end{split}
\end{equation}
where $ \omega_{f}$ and $ \omega_{0}$ are the frequencies of the field and the two-level transition, respectively; $\hat{a}^\dagger (\hat{a})$ is the creation (annihilation) operator of the field; and $\hat{\sigma}_{i} $ ($ i= z,+,-$) are Pauli matrices. The nonlinear coupling is split into a rotating-wave term
described by the coupling strength $ g_{-}$ and a counter-rotating term described by the coupling strength $ g_{+}$.

When $g_+=0$, we restore the original BS model without the counter-rotating term.
In this case, the energy is lower-bound
only if $g_- \le 1/(2S)$, where $S$ denotes the eigenvalues of the constant of motion of the Casimir operator $\hat{S}^2$. In terms of the Pauli operators, $\hat{S}^2 = \hat{S}_{\mp}\hat{S}_{\pm} + \hat{S}^2_z \pm \hat{S}_z$.
When $g_- > 1/(2S)$, the eigenenergies can be arbitrarily low. As there is no stable ground state, the Hamiltonian is unbounded and no ground state exists. In this sense, the model appears to be incomplete for $g_- > 1/(2S)$. The ground-state energy decreases indefinitely when the eigenvalues of $\hat{C}=\hat{a}^\dagger\hat{a}+2\hat{S}_z$, denoted as $c$, increase above a critical value. This effect illustrates a dramatic phase transition
between the normal phase and the super-radiant phase~\cite{Cordeiro2008}.

When the counter-rotation term is included ($g_+\ne 0$), the Hamiltonian cannot be analytically diagonalized, but the eigenstates and eigenenergies can be computed by a numerical approach.
After expanding the initial state in terms of the numerically determined eigenstates and eigenenergies, one can temporally evolve the physical quantities.
The following section proposes another formal and systematic numerical method, which does not require precalculation of the stationary eigenstates and eigenenergies.

\section{Non-Perturbative Time Evolution}
\label{sec:method}

The dynamics of a Hamiltonian system are given by the time-dependent Schr\"{o}dinger equation,
\begin{equation}
i\hbar \frac{\partial}{\partial t} \vert  t\rangle=\hat{H}\vert  t\rangle,
\end{equation}
which is formally solved as
\begin{equation}\label{TaylorE}
\vert  t\rangle=e^{-i\hat{H}t}\vert 0\rangle
=\sum_{n=0}^{\infty} \frac{1}{n!}(-i\hat{H}t)^n \vert 0\rangle.
\end{equation}
Here, $\vert 0\rangle$ denotes the initial ground state $\vert t= 0\rangle$.
Defining
 \begin{equation}
 \vert B _{n}\rangle=\frac{1}{n!}(-i\hat{H}t)^n \vert 0\rangle,
 \end{equation}
we notice that $\vert B _{n}\rangle$ can be concisely iterated as
 \begin{equation}
 \label{iter-1}
 \vert B _{n+1}\rangle=-\frac{it}{n+1}\hat{H}\vert B _{n}\rangle.
 \end{equation}
$\hat{H}$ can be easily operated on the initial state $\vert 0\rangle$ when $\vert 0\rangle$ is expanded in the Fock state space $\{\vert p\rangle\}$, $p=0, 1, 2,...$ or the coherent space $\{\vert \alpha \rangle\}$.

In the following analysis, we illustrate the proposed numerical recipe on a BS model. The time-dependent wave function of the two-level system is expanded as an excited state $\vert e\rangle$ and the ground state $\vert g\rangle$, namely, as $\vert  t\rangle =\vert t\rangle_{1}\vert e\rangle+\vert t\rangle_{2}\vert g\rangle $, where $\vert t \rangle_{1}$ and $\vert t \rangle_{2}$ are the time-dependent wavefunctions in the $\vert e\rangle$ and $\vert g\rangle$ states, respectively. Furthermore, $\vert t \rangle$ is expanded as the following Taylor series:
\begin{equation}
\label{tstate}
\begin{split}
\vert t\rangle &=e^{-i\hat{H}t}(\vert t\rangle_{1}\vert e\rangle+\vert t\rangle_{2}\vert g\rangle)\\
&=\sum_{n=0}^\infty
\frac{1}{n!}(-i\hat{H}t)^n(\vert 0\rangle_{1}\vert e\rangle+\vert 0\rangle_{2}\vert g\rangle).
\end{split}
\end{equation}
Defining
 \begin{equation}
 \label{bn1}
 \vert B _{n}\rangle=\frac{1}{n!}(-i\hat{H}t)^n(\vert 0\rangle_{1}\vert e\rangle+\vert 0\rangle_{2}\vert g\rangle),
 \end{equation}
we recover the iteration relation Eq.~(\ref{iter-1}).
Expanding $ \vert B _{0}\rangle$ as $\vert 0\rangle_{1}\vert e\rangle+\vert 0\rangle_{2}\vert g\rangle $, Eq.~(\ref{bn1}) can be rewritten as the following operator iteration relation:
\begin{equation}
\label{bn3}
\vert B_{n}\rangle= \frac{1}{n!}(-i\hat{H}t)^n\vert B _{0}\rangle.
\end{equation}

Once $\vert B_{n}\rangle$ is known, the time evolution of the state is clearly given by
\begin{equation}
\label{bn3-1}
\vert t\rangle=\sum_n \vert B_n \rangle.
\end{equation}
The initial states can then be expanded in the orthogonal complete Fock basis $\{ \vert p\rangle\}$ as
\begin{equation}
\label{bn2}
\vert B _{n}\rangle= \sum_{p=0}^\infty (f_{p}^{1,n}(t)\vert e\rangle+f_{p}^{2,n}(t)\vert g\rangle)\otimes \vert p\rangle.
\end{equation}
Substituting $\vert B _{0}\rangle= \sum_{p=0}^\infty (f_{p}^{1,0}\vert e\rangle+f_{p}^{2,0}\vert g\rangle)\otimes \vert p\rangle,$ into Eq.~(\ref{bn3}),
applying $\hat{a}\vert p\rangle=\sqrt{p-1}\vert p\rangle$ and $\hat{a}^\dagger \vert p\rangle=\sqrt{p}\vert p\rangle$ (which appear in $\hat{H}^n$ over the Fock basis $\vert p\rangle$), and comparing the result with Eq.~(\ref{bn2}), we obtain the following iterative relationship between the corresponding expansion coefficients $ f_{p}^{i,n} (i=1,2; p=0,1,2,...) $,
\begin{equation}
\label{ite-10}
\begin{bmatrix} f_{0}^{1,n}(t) \\ f_{1}^{1,n}(t)\\.\\.\\.\\ f_{0}^{2,n}(t) \\ f_{1}^{2,n}(t)\\.\\.\\. \end{bmatrix} = \frac{(-it)^n}{n!} Q^n \begin{bmatrix} f_{0}^{1,0}(0) \\ f_{1}^{1,0}(0)\\.\\.\\. \\f_{0}^{2,0}(0)\\f_{1}^{2,0}(0)\\.\\.\\. \end{bmatrix},
\end{equation}
where $Q$ is the transfer matrix, given by
\begin{equation}
Q=\begin{bmatrix} Q_{11} & Q_{12} \\ Q_{21} & Q_{22} \end{bmatrix}, \quad
\end{equation}
and
\begin{equation}\label{Qg}
\begin{split}
 &Q_{11}=(\omega_{0}/2+\omega_{f}p)\delta_{p,p^{'}},\\
 &Q_{22}=(-\omega_{0}/2+\omega_{f}p)\delta_{p,p^{'}},\\
 &Q_{12}=g_{+}(p+1)\delta_{p+1,p^{'}}+g_{-}p\delta_{p-1,p^{'}},\\ &Q_{21}=g_{-}(p+1)\delta_{p+1,p^{'}}+g_{+}p\delta_{p-1,p^{'}}.
\end{split}
\end{equation}
This formulation transforms the operator iteration relations between the $\vert B _{n}\rangle $ into matrix relations between the $ f_{p}^{i,n} (i=1,2)$. As the Hamiltonian multiplication $\hat{H}^n$ is absorbed in the matrix multiplication $Q$ as $Q^n$, the number of numerical calculations is greatly reduced, gaining an additional advantage. In this formulation, the matrix multiplication $Q^n$ is uniquely determined by the Hamiltonian and is independent of the $ f_{p}^{i,0} (i=1,2; p=0,1,2,...) $ determined by the initial states. The time-dependent properties of the system under different initial states depend on the same $Q^n$. Therefore, the matrix multiplication $Q^n$ can be calculated and stored prior to defining the arbitrary initial states.

To avoid error accumulation during the time evolution, we discrete the time $t$ into $ K $ steps, where each unit period $ \Delta t $ satisfies $ t=K\Delta t $. The term
\begin{equation}
\label{mt}
M(t)=\left[\frac{(-i \Delta t)^n}{n!}Q^n\right]^{K}, K=0,1,2,...
\end{equation}
is saved after each calculation. When $ M(t) $ acts on $ f_{p}^{i,0} (i=1,2; p=0,1,2,...)$ for different $ K $, $ f_{p}^{i,n}(t) (i=1,2; p=0,1,2,...) $ at each time is obtained as
\begin{equation}
\label{Mt1}
\begin{bmatrix} f_{0}^{1,n}(t) \\ f_{1}^{1,n}(t)\\.\\.\\.\\ f_{0}^{2,n}(t) \\ f_{1}^{2,n}(t)\\.\\.\\. \end{bmatrix} = M(t) \begin{bmatrix} f_{0}^{1,0}(0) \\ f_{1}^{1,0}(0)\\.\\.\\. \\f_{0}^{2,0}(0)\\f_{1}^{2,0}(0)\\.\\.\\. \end{bmatrix}.
\end{equation}

In summary, once the initial state is given, $ f_{p}^{i,0}$, $i=1,2$; $p=0,1,2,...$ can be determined, and the subsequent $f_{p}^{i,n}(t), i=1,2; p=0,1,2,...$ can be calculated by Eqs.~(\ref{ite-10})--(\ref{Mt1}). Numerically, we must truncate $n$ and $p$. Different combinations of the truncated forms, denoted as $N$ and $P$, respectively, were tested in the present study. The results with truncation errors below $10^{-9}$ are presented below.

Now, assume a coherent initial state of the system:
\begin{equation}
\vert 0\rangle=e^{-\alpha^2/2}\sum_{p} \frac{\alpha^p}{\sqrt{p!}}\vert p\rangle (\sin\theta \vert e\rangle + \cos\theta \vert g\rangle),
\end{equation}
where the amplitude $\alpha $ of the coherent state and the mean photon number of the initial state are related: $\alpha^2=\langle p \rangle$. Given $ \vert B _{0}\rangle=\vert 0\rangle_{1}\vert e\rangle+\vert 0\rangle_{2}\vert g\rangle $ and Eq.~(\ref{bn2}), $ f_{p}^{i,0}(0)$ for $i=1, 2$ are obtained as
\begin{equation}\label{initialf}
\begin{split}
&f_{p}^{1,0}(0)=e^{-\alpha^2/2}\frac{\alpha^p}{\sqrt{p!}}\sin\theta,\\
&f_{p}^{2,0}(0)=e^{-\alpha^2/2}\frac{\alpha^p}{\sqrt{p!}}\cos\theta.
\end{split}
\end{equation}

The resulting wavefunction $ \vert t\rangle $, given by
\begin{equation}
\label{wf}
\vert t\rangle=\sum_{n}\sum_{p} (f_{p}^{1,n}(t) \vert p\rangle \vert e\rangle +f_{p}^{2,n}(t) \vert p\rangle \vert g\rangle),
\end{equation}
can be numerically calculated by the above procedure.
In the following analysis, we investigate two typical measurements of interest in the quantum optics community. The first is the mean photon number $ \langle \hat{n} \rangle $, expressed as
\begin{eqnarray}
\label{photonum}
\langle \hat{n} \rangle &=& \langle t \vert \hat{a}^{\dagger}\hat{a} \vert t\rangle \nonumber\\
&=& \sum_{m}\sum_{q}\left[f_{q}^{1,m*}(t) \langle q \vert \langle e \vert +f_{q}^{2,m*}(t) \langle q \vert \langle g \vert \right]\hat{a}^{\dagger}\hat{a} \nonumber\\
&&\times\sum_{n}\sum_{p} [f_{p}^{1,n}(t) \vert p\rangle \vert e\rangle +f_{p}^{2,n}(t) \vert p\rangle \vert g\rangle]\nonumber\\ &=&\sum_{m,n}\sum_{p}p[f_{p}^{1,m*}(t)f_{p}^{1,n}(t)
+f_{p}^{2,m*}(t)f_{p}^{2,n}(t)].
\end{eqnarray}

The second is the atomic inversion $ \langle \hat{\sigma}_{z} \rangle $:
\begin{eqnarray}
\label{inversion}
\langle \hat{\sigma}_{z} \rangle &= &\langle t \vert \hat{\sigma}_{z} \vert t\rangle \nonumber\\
&=& \sum_{m}\sum_{q}[f_{q}^{1,m*}(t) \langle q \vert \langle e \vert +f_{q}^{2,m*}(t) \langle q \vert \langle g \vert ]\nonumber\\
&&\times(\vert e\rangle\langle e \vert - \vert g\rangle\langle g \vert )\nonumber\\
&&\times\sum_{n}\sum_{p} [f_{p}^{1,n}(t) \vert p\rangle \vert e\rangle +f_{p}^{2,n}(t) \vert p\rangle \vert g\rangle]\nonumber\\ &=&\sum_{m,n}\sum_{p}[f_{p}^{1,m*}(t)f_{p}^{1,n}(t)-f_{p}^{2,m*}(t)
f_{p}^{2,n}(t)].
\end{eqnarray}
Physically, the mean photon number and the atomic inversion may periodically collapse and revive, as observed in the JC model.

For comparison, we briefly describe the time evolution from the known eigenvalues and eigenvectors of $\hat{H}$ (hereafter, this method is abbreviated as the TEEE method). The stationary state of Eq.~(\ref{H}) is assumed as $ \vert \psi\rangle=\sum_{p=0}^{\infty}[a_{n}\vert p \rangle\vert e \rangle+b_{n}\vert p \rangle\vert g \rangle $, where $ \vert p \rangle$ is the Fock state in the Fock representation. Solving $ \hat{H}\vert \psi\rangle=E\vert \psi \rangle $, we get
\begin{equation}
\label{traditional}
\begin{split}
&(\omega_{f}p+\frac{\omega_{0}}{2})a_{p}+g_{-}(p+1)b_{p+1}+g_{+}pb_{p-1}=Ea_{p},\\
&(\omega_{f}p-\frac{\omega_{0}}{2})b_{p}+g_{+}(p+1)a_{p+1}+g_{-}pa_{p-1}=Eb_{p}.
\end{split}
\end{equation}
After numerically solving the above equations, we obtain the eigenenergy spectrum $ \{E_{j}\} $ and the corresponding stationary-state set $ \{a_{p}^{j},b_{p}^{j}\} $, where $j=0,1,2,...,2P+2$ with $P$ being the truncation of the integer number $p$.
Then, expanding the given initial state $ \vert 0\rangle $ in terms of the stationary-state set $ \vert 0\rangle=\sum_{j=0}^{2P+2} F_{j}\vert  \psi_j\rangle $, where $ F_{j} $ is the expansion coefficient, the time evolution of the wavefunction is obtained as
\begin{eqnarray}\label{Eq:TEEE}
\vert t\rangle &=&e^{-i\hat{H}t}\vert  0\rangle\nonumber\\
&=&\sum_{j=0}^{\infty}F_{j}e^{-iE_{j}t}\sum_{p=0}^{\infty}\left[a_{n}\vert p \rangle\vert e \rangle+b_{n}\vert p \rangle\vert g \rangle\right]\nonumber\\
&\approx&\sum_{j=0}^{2P+2}F_{j}e^{-iE_{j}t}\sum_{p=0}^{P}\left[a_{n}\vert p \rangle\vert e \rangle+b_{n}\vert p \rangle\vert g \rangle\right].
\end{eqnarray}
As in Eqs.~(\ref{wf})--(\ref{photonum}), we can calculate the time evolutions of the mean photon number $ \langle \hat{n} \rangle=\langle t \vert \hat{a}^{\dagger}\hat{a} \vert t\rangle $ and the atomic inversion $ \langle \hat{\sigma}_{z} \rangle=\langle t \vert \hat{\sigma}_{z} \vert t\rangle $.

Here, we emphasize that both the Taylor expansion in Eq.~(\ref{TaylorE}) and the basis expansion in Eq.~(\ref{bn2}) include an infinite sum of terms. Therefore, the accuracy of the time evolution can be systematically refined by including more Taylor series and Fock states, respectively. The expansion converges quickly due to the factorial of $n$ in the denominator. This approach, which we call the non-perturbation time-evolution method, offers several advantages over expanding the eigensolutions (i.e., the TEEE method). First, this method is directly applicable when the eigenenergy and eigenstate are unknown. Second, it systematically obtains the long-time evolution. Finally, the expansion coefficients $M(t)$ in Eq.~(\ref{Mt1}) need not be recalculated in the time evolutions of different initial states.

\section{BS model with RWA $ (g_{+}=0) $}
\label{sec:dis1}
We now apply the above numerical method to the BS model and compare the results with those of the TEEE method and Rodr\'{i}guez-Lara {\it et al.}~\cite{Rodrguez2013}. We first study the BS model with RWA ({\it i.e.}, $g_{+}=0$) at resonance, setting $ \omega_{f}=\omega_{0}=1 $. In this case, the BS model is exactly solvable. Figure~\ref{Fig1} shows the time evolutions of the mean photon number $\langle \hat{n} \rangle $ and the atomic inversion $ \langle \hat{\sigma}_{z} \rangle $, starting from a coherent initial state with $ \alpha=5 $.
Physically, the atomic inversion $ \langle \hat{\sigma}_{z} \rangle $ is equivalent to the atomic intensity difference between the two levels $\langle e \rangle $ and $\langle g \rangle $.
The atomic inversion collapses and revives with the variable photon number. Under the truncations $ N=20 $ and $ P=50 $, our results accorded with those of Rodr\'{i}guez-Lara {\it et al.}~\cite{Rodrguez2013} and the TEEE method (truncated with $P=60$). Similar agreement was observed in other parameter spaces (data not shown). Note that our method can be extended to very long-time evolution. Our tests confirmed highly accurate time-evolution behavior up to $t=600$ (in units of $\omega_f$), and the evolution time is extendible if necessary.
\begin{figure}[htbp]
\begin{center}
 \centering
\includegraphics[width=0.5\textwidth]{Fig1a.eps}
\includegraphics[width=0.5\textwidth]{Fig1b.eps}
\end{center}
 \caption{Evolutions of (a) mean photon number $\langle \hat{n} \rangle $ and (b) atomic inversion $ \langle \hat{\sigma}_{z} \rangle $ in the BS model with RWA. The initial state is the coherent state $ \vert t=0\rangle = \vert \alpha_{+},g\rangle + \vert \alpha_{-},e\rangle $ with $ \alpha=5 $ in the BS model when $ g_{+}=0 $ at resonance $ \omega_{f}=\omega_{0}=1$.  The coupling strength is $ g_{-}=0.1\omega_{f} $. For comparison, the results of Eqs.~(\ref{ite-10})--(\ref{Mt1}) (black line) and the TEEE method (blue points) are also plotted.
} \label{Fig1}
\end{figure}

\section{BS model without RWA $(g_{+}\ne 0) $}
\label{sec:dis2}
In this section, we confirm that our method is also applicable to RWA-void BS model, {\it i.e.}, when $g_{+}\ne 0$. In this case, the BS model cannot be solved analytically. As before, the non-perturbation method is implemented by Eqs.~(\ref{ite-10})--(\ref{Mt1}). While our method must only recalculate $M(t)$ for the different $Q$s in Eq.~(\ref{Qg}), the TEEE method must numerically recalculate the eigenvalues and eigenvectors of the Hamiltonian.
Figure~\ref{Fig2} shows the evolutions of the atomic inversion $ \langle \hat{\sigma}_{z} \rangle $ and mean photon number $ \langle \hat{n} \rangle $ in RWA-void BS model, starting from the vacuum initial state $ \vert t=0\rangle = \vert 0,e\rangle $ with parameters $ \omega_{f}=1 $, $ \omega_{0}=3/4\omega_{f} $ and $ g_{-}=g_{+}=0.4\omega_{f} $ (corresponding to $\alpha=0$ and $\theta=\pi/2$ in Eq.~(\ref{initialf})).

Our results (truncated by $ N=30 $ and $ P=50 $) agree well with those of Rodr\'{i}guez-Lara {\it et al.}~\cite{Rodrguez2014} and the TEEE method (truncated by $P=60$),
demonstrating the applicability of our method to the RWA-void BS model. Note that the collapse and revival phenomena are absent in Fig.~\ref{Fig2}, as they were removed by the counter-rotating-wave term.
\begin{figure}[htbp]
\begin{center}
 \centering
\includegraphics[width=0.5\textwidth]{Fig2a.eps}
\includegraphics[width=0.5\textwidth]{Fig2b.eps}
\end{center}
 \caption{Evolution of (a) mean photon number $\langle \hat{n} \rangle $ and (b) atomic inversion $ \langle \hat{\sigma}_{z} \rangle $ in the RWA-void BS model. The initial state is a vacuum $ \vert t=0\rangle =\vert 0,e\rangle $ corresponding to $\alpha=0$ and $\theta=\pi/2$ in Eq.~(\ref{initialf}) (meaning that no photons are excited at $t=0$). Other parameters are $ \omega_{f}=1 $, $ \omega_{0}=3/4\omega_{f} $ and $ g_{-}=g_{+}=0.4\omega_{f} $. The results of Eqs.~(\ref{ite-10})--(\ref{Mt1}) (black line) and the TEEE method (blue points) are plotted for comparison.
 } \label{Fig2}
\end{figure}

We subsequently consider the RWA-void BS model under the resonance condition $ \omega_{f}=\omega_{0}=1 $ with parameter values $ g_{-}=g_{+}=2\omega_{f} $. In this case, the Hamiltonian has no ground state and therefore has no physical meaning~\cite{Ng2000}. The aim is to confirm the behaviors of the mean atomic inversion $ \langle \hat{\sigma}_{z} \rangle $ and the mean photon number $\langle \hat{n} \rangle $, and the applicability of the numerical method. Rodr\'{i}guez-Lara {\it et al.}~\cite{Rodrguez2013} calculated the physics of this system via quantum optics methods. They assumed an analogy between the transport of single-photon states and the propagation of a classical field. The numerical propagation was stable on a 2000 $\times$ 2000 photonic lattice. However, the ground-state energy of the Hamiltonian is unbounded in this parameter space. In the na\"{\i}ve sense, this result implies an unstable mean atomic inversion and mean photon number. In the following, we first reproduce Rodr\'{i}guez-Lara {\it et al.'s} results by the TEEE method, then reveal the contradiction between the two arguments.

In practice, Eq.~(\ref{Eq:TEEE}) is truncated by the limit of $P$ instead of being implemented in infinite dimensional Fock space. The numerical procedures of the TEEE method should therefore converge as the truncation number increases. In this case, the final results can be reliably obtained while satisfying the given precision requirements. As shown in panels (a) and (b) of Fig.~\ref{Fig3}, the mean photon number and mean atomic inversion results that were truncated by $P=2000$ matched those obtained by the TEEE method of Rodr\'{i}guez-Lara {\it et al.}~\cite{Rodrguez2013} with the same truncation factor.
However, when the truncation limit was increased to 3000 (see Fig.~\ref{Fig3}(c) and (d)), the results of our method and TEEE were very different (In order to see clearly the difference, we showed only the results for $t\in [0,40]$ comparing to Fig.~\ref{Fig3}(a) and (b)). Hence, whether the truncation requires to be increased or whether the discrepancy naturally arises from the non-physicality of the unbounded ground-state energy due to the incompleteness of the model must be elucidated in further analysis.
\begin{figure}[htbp]
\begin{center}
 \centering
\includegraphics[width=0.5\textwidth]{Fig3a.eps}
\includegraphics[width=0.5\textwidth]{Fig3b.eps}
\includegraphics[width=0.5\textwidth]{Fig3c.eps}
\includegraphics[width=0.5\textwidth]{Fig3d.eps}
\end{center}
 \caption{Evolutions of (a) mean photon number $\langle \hat{n} \rangle$ and (b) atomic inversion $ \langle \hat{\sigma}_{z} \rangle $ for $P=2000$. (c) and (d) represent the same data for (a) and (b), respectively, but for $P=2000$ and $3000$ simultaneously. The BS model begins from a vacuum state $ \vert t=0\rangle =\vert 0,e\rangle $ with $ \omega_{f}=\omega_{0}=1 $ and $ g_{-}=g_{+}=2\omega_{f} $.
 } \label{Fig3}
\end{figure}
\begin{figure}[htbp]
\begin{center}
 \centering
\includegraphics[width=0.5\textwidth]{Fig4a.eps}
\includegraphics[width=0.5\textwidth]{Fig4b.eps}
\end{center}
 \caption{Top panel: Difference between the j-th ($ j=1,2,3,... $) excited-state energy and the ground-state energy ($\Delta E=E_{j}-E_{0} $) in the RWA-void BS model with parameters $ \omega_{f}=\omega_{0}=1 $ and $ g_{-}=g_{+}=2\omega_{f} $, calculated by the TEEE method truncated by $P= 2000$ (red line) and $P=3000$ (blue line). Bottom panel: Energy-level differences for $P=2000$ (red line in bottom left) and $P=3000$ (blue line in bottom right panel), respectively.
 } \label{Fig4}
\end{figure}

To resolve the above doubt, Fig.~\ref{Fig4}(a) shows the energy-level differences ($ \Delta E=E_{j}-E_{0} $) between the $j$-th ( $j=1,2,3,...$) excited state and the ground state in the RWA-void BS model. Figure~\ref{Fig4}(b) compares the ($ \Delta E=E_{j}-E_{0}$ in the models truncated by $P=2000$ and $P=3000$. Comparison of the two results clearly shows the invalidity of simply increasing the truncation limit, indicating that the non-convergence of the physical quantities originates from the incomplete Hamiltonian and the unbounded ground-state energy.

To further confirm this judgment, we studied the ground-state energy as a function of truncation. Figure~\ref{Fig5}(a) shows the results for $ g_{-}=g_{+}=0.4\omega_{f}$. After an initial rapid decrease, the ground-state energy stabilized at $P>6$.
Moreover, as shown in Fig.~\ref{Fig5}(b), the ground-state energy decreased linearly with increasing truncation limit for $ g_{-}=g{+}=2\omega_{f}$. In other words, the investigated system lacked a ground state and was physically infeasible.

\begin{figure}[htbp]
\begin{center}
 \centering
\includegraphics[width=0.5\textwidth]{Fig5a.eps}
\includegraphics[width=0.5\textwidth]{Fig5b.eps}
\end{center}
 \caption{Ground-state energy versus truncation factor in the BS model with parameters (a) $ \omega_{f}=1 $ and $ \omega_{0}=3/4\omega_{f} $, $ g_{-}=g_{+}=0.4\omega_{f} $, and (b) $ \omega_{f}=\omega_{0}=1 $ and $ g_{-}=g{+}=2\omega_{f} $.
 } \label{Fig5}
\end{figure}

To incorporate the environmental influences in real-system dynamics, we must consider the irreversible processes. For this purpose, we included a purely imaginary coupling term comprising two phenomenological dissipative terms~\cite{non-Hermitian1} in the Hamiltonian~(\ref{H}) of the BS model as follows:
\begin{equation}
\hat{H}_{\rm dissipative}=\hat{H}-i\beta\hat{a^{+}}\hat{a}-i\gamma \vert e \rangle \langle e \vert,
\end{equation}
where $ \beta $ and $ \gamma $ are the photon leakage rate from the cavity and the spontaneous emission rate, respectively. The last (non-Hermitian) term renders the Hamiltonian as non-Hermitian. The reliability of the non-Hermitian Hamiltonian approach has been already tested in the dynamics of the pumped-dissipative JC model~\cite{non-Hermitian2}.
As our numerical method expands the initial state in Fock space, it is directly applicable to systems of dissipations. This property is absent in the TEEE method. Figure~\ref{Fig6} shows the evolution of the atomic inversion $ \langle \hat{\sigma}_{z} \rangle $ and the mean photon number $ \langle \hat{n} \rangle $ in this system, starting from a vacuum state $ \vert t=0\rangle = \vert 0,e\rangle $ (corresponding to $\alpha=0$ and $\theta=\pi/2$ in Eq.~(\ref{initialf})). Here, we set $ \beta=\gamma=0.01\omega_{f} $, and the other parameters as $ \omega_{f}=1 $, $ \omega_{0}=3/4\omega_{f} $, and $ g_{-}=g_{+}=0.4\omega_{f} $. Owing to finite damping of the photons and atoms, the oscillating photon number and atomic inversion in the dissipation case tended to steady values (see Fig.~\ref{Fig6}(a) and (b), respectively). The damping of the oscillations is attributable to interactions with the environment.

\begin{figure}[htbp]
\begin{center}
 \centering
\includegraphics[width=0.5\textwidth]{Fig6a.eps}
\includegraphics[width=0.5\textwidth]{Fig6b.eps}
\end{center}
 \caption{Evolutions of (a) mean photon number $\langle \hat{n} \rangle $ and (b) atomic inversion $ \langle \hat{\sigma}_{z} \rangle $ of the system for $ \beta=\gamma=0.01\omega_{f}$, starting from a vacuum initial state $ \vert t=0\rangle =\vert 0,e\rangle $. The other parameters are $ \omega_{f}=1 $, $ \omega_{0}=3/4\omega_{f} $ and $ g_{-}=g_{+}=0.4\omega_{f} $. For comparison, the results in the no-dissipation case $\beta=\gamma=0$ are plotted in black).
 } \label{Fig6}
\end{figure}

\section{\label{sec:level5}Conclusion}
\label{sec:con}
We proposed a systematic numerical method that calculates the time evolution of the physical quantities in the BS model of the interactions between a two-level atomic system and a field. Applying this method, we successfully transformed the Hamiltonian multiplication of $\hat{H}^n$ into a matrix multiplication of $Q^n$. This transformation greatly reduced the numerical calculation involved. Moreover, the matrix multiplication $Q^n$ is independent of the initial states, meaning that the time-dependent properties of the system can be computed for arbitrary initial states without recalculating the matrix multiplication $Q^n$.

The performance of the method was tested in numerical comparisons with the TEEE method. The method performed well in the whole coupling regime, and captured the collapses and revivals of the atomic inversion along with the photon-number variations. In the parameter space of a non-real positive-definite Hamiltonian, we also resolved the uncertainty regarding whether the evolution was stabilized by increasing the truncation limit over the dimension of the Fock space or whether it was a natural consequence of the unbounded ground-state energy (a nonphysical situation caused by the incompleteness of the Hamiltonian). We found that when the ground-state energy of the Hamiltonian is unbounded in a certain parameter space, the time evolution of the physical quantities is naturally unstable. Furthermore, the performance of the numerical method was tested in a dissipative BS model, which represents open quantum systems.

\begin{acknowledgments}
This work was supported by the NSF of China under Grant Nos. 11774316 and 11835011.
\end{acknowledgments}

\newpage
\small