\section{Introduction and the main results}\label{sec1}
In this paper we deal with  the nonexistence of nonnegative solutions of the following parabolic  problems
\begin{equation}\label{main}
\begin{cases}
&u_t-L_\mathcal A u\geq  (K\ast u^p)u^q \quad\mbox{ in } \mathbb R^N \times \mathbb (0,\infty),\, N\geq 1\\
&u(x,0) = u_0(x)\ge0 \,\, \text{ in } \mathbb R^N,
\end{cases}
\end{equation}
and
\begin{equation}\label{main2}
\begin{cases}
&u_t+L_\mathcal A u\geq  (K\ast u^p)u^q \quad\mbox{ in } \mathbb R^N \times \mathbb (0,\infty),\\
&u(x,0) = u_0(x)\ge0 \,\, \text{ in } \mathbb R^N, N\geq 1,\end{cases}
\end{equation}
where $u_0\in L^1_{loc}({\mathbb R}^N)$, $u_0\geq 0$  and the differential operator $L_{\mathcal{A}}u= \text{div} \mathcal{A}(x,u,\nabla u)$ is such that the mapping  $\mathcal{A}: \mathbb R^N \times \mathbb{R}\times \mathbb{R}^N \to \mathbb{R}^N$ is a Carath\'{e}odory function with
\begin{equation}\label{A_cal}\mathcal{A}(x,z,0) = 0,\qquad  \mathcal{A}(x,z,\xi)\cdot \xi\geq 0\end{equation}
for every $(x,z,\xi) \in \mathbb{R}^N\times [0, \infty)\times\mathbb{R}^N$. In addition, $L_{\mathcal{A}}$  is assumed to be  a weakly $m$-coercive operator, that is,  there exist  a  constant $c_0> 0$ and an exponent $m>1$
such that the inequality \begin{equation} \label{weak_m}
 \mathcal A(x,z,\xi)\cdot \xi\ge c_0 |\mathcal A(x,z,\xi)|^{m'}
\end{equation}
holds for all $(x,z,\xi) \in \mathbb{R}^N\times [0, \infty)\times\mathbb{R}^N$ with $m'= m /(m-1)$.

The main prototypes for $\mathcal{A}$ are:
\begin{itemize}
\item the $m$-Laplace operator given by $\mathcal{A}(x,z,\xi)=|\xi|^{m-2}\xi$, $m>1$;
\item the $m$-mean curvature operator, or the generalized mean curvature operator, given by
$$
\mathcal{A}(x, z, \xi)= \frac{|\xi|^{m-2}}{\sqrt{1+|\xi|^m}}\xi.
$$
\end{itemize}
Furthermore, on the right-hand side  of \eqref{main} and \eqref{main2} we assume the exponents $p,q>0$. The function  $K\in C(\mathbb{R}^N\setminus\{0\})$, $K>0$ satisfies $\liminf_{x\to 0}K(x)>0$  and there exists $\rho>0$ and $0<\beta<m/2$ such that
\begin{equation}\label{la}
K(x)\geq c|x|^{-\beta} \quad \mbox{ for all }\; x\in \mathbb{R}^N, \, |x|>\rho,
\end{equation}
where $c>0$ is a positive constant. Also, by $K\ast u^p$ we denote the standard convolution operator defined by
$$
(K\ast u^p)(x,t)=\int_{\mathbb R^N}K(x-y)u^p(y,t) dy\quad\mbox{ for all }(x,t)\in\mathbb R^N\times (0, \infty).
$$
We are interested in  {\it nonnegative weak solutions} of \eqref{main},  that is nonnegative functions $u(x,t)$, belonging to the class $\mathcal S$ given by those  $u\in  {W}^{1,m}_{\rm loc}(\mathbb R^N\times (0, \infty))$ which fulfil the two conditions below
\begin{enumerate}
\item[(i)] ${\mathcal{A}(x,u,\nabla u)  \in  [{L}^{m'}_{\rm loc}(\mathbb R^N\times (0, \infty))]^N}$,
\item[(ii)] $(K\ast u^p)u^q  \in  {L}^1_{\rm loc}(\mathbb R^N\times (0, \infty))$,
\end{enumerate}
and such that  for any nonnegative test function  $\varphi \in C^{1}_{c}(\mathbb R^N\times \mathbb{R})$, we have
\begin{equation} \label{weak_form2}
\int_0^\infty\!\int_{\mathbb R^N} (K\ast u^p)u^q  \varphi \, dx \, dt \leq  \int_0^\infty\!\int_{\mathbb R^N} u_t \; \varphi \, dx \, dt
 + \int_0^\infty\!\int_{\mathbb R^N} \mathcal{A}(x,u,\nabla u) \cdot \nabla \varphi \, dx \, dt,
\end{equation}
or equivalently
\begin{equation} \label{weak_form1}
\begin{aligned}
\int_0^\infty\!\!\int_{\mathbb R^N} (K\ast u^p)u^q \varphi \, dx \, dt \leq & - \int_{{\mathbb{R}^N}} \, u_0(x) \varphi(x,0) \, dx - \int_0^\infty\!\!\int_{\mathbb R^N} u \, \varphi_{t} \, dx \, dt \\
& + \int_0^\infty\!\!\int_{\mathbb R^N} \mathcal{A}(x,u,\nabla u) \nabla\varphi \, dx \, dt. \end{aligned}
\end{equation}

For {\it nonnegative weak solutions} of \eqref{main2} the functions $u$ in class $\mathcal{S}$ satisfy (i)-(ii) above and we replace \eqref{weak_form2} by
\begin{equation} \label{weak_form3}
\int_0^\infty\!\int_{\mathbb R^N} (K\ast u^p)u^q  \varphi \, dx \, dt \leq  \int_0^\infty\!\int_{\mathbb R^N} u_t \; \varphi \, dx \, dt
 - \int_0^\infty\!\int_{\mathbb R^N} \mathcal{A}(x,u,\nabla u) \cdot \nabla \varphi \, dx \, dt.
\end{equation}

The pioneering paper in this subject is  due to Fujita \cite{Fuj} in $1966$, where he investigated the fol\-lo\-wing Cauchy problem for the semilinear heat equation
\begin{equation}\label{fujita66}
\begin{cases}
&u_t - \Delta u = u^{q} \,\,\, \text{in }\mathbb{R}^{N} \times (0, \infty),\\
&u(x,0) = u_0(x) \,\, \text{in } \mathbb{R}^{N}.\\
\end{cases}
\end{equation}
Fujita obtained the critical exponent $ q_{F} = 1 + 2/N$ on the existence versus nonexistence of nonnegative nontrivial solutions. Precisely,  nonexistence of solutions, that is blow up, holds when $1<q< 1 + 2/N$ and $u_0$ is bounded and nonnegative, while blow up can occur when
 $q> 1 + 2/N$ depending on the size of $u_0$. Since then, there have been a number of extensions of Fujita results in many directions.
In particular, the result obtained by Fujita was completed, relatively to the critical case, in \cite{Hay} for $N=1,2$ and in \cite{KST} for $N\geq 3$.

When one moves from the equality to the inequality case, in contrast to elliptic problems, the critical exponents for parabolic equations and inequalities coincide,  at least in the semilinear case.  This observation  is due to Mitidieri and Pokhozev in \cite[Part II]{MPbook} where a powerful tool on a priori estimates is devised. Specifically,
the technique of the nonlinear capacity developed by the authors in \cite{MPbook} is based on the derivation of suitable a priori bounds on the solution of the problem under consideration by carefully choosing special test functions and scaling argument (other types of a priori estimates for quasilinear parabolic problems may be found in \cite{LS2009, LS2013, LSS2011}). This fact enables one to analyze by the same method nonexistence of solutions, not necessarily nonnegative, of essentially different problems, namely noncoercive parabolic problems
$$u_t-\Delta u\ge |u|^q \quad\text{in }\mathbb{R}^{N} \times (0, \infty),$$
and coercive parabolic problems respectively given by
\begin{equation}\label{coer_mod}u_t+\Delta u\ge |u|^q \quad\text{in }\mathbb{R}^{N} \times (0, \infty).
\end{equation}
As noted in \cite[Remark 3.8]{MP_JEv}, it is possible to show that in this case the problem \eqref{coer_mod} has no non negative weak solutions for every $q > 1$.

We also refer the reader to \cite[Section 18]{QS} where the authors propose a modification of the above arguments, based on rescaling of test functions.

Later on, several extensions to quasilinear problems were developed in literature: see e.g., the book by
Samarskii, Galaktionov, Kurdyumov and Mikhailov \cite{SGKM} and the work by Levine, Lieberman and Meier \cite{LLM} on parabolic problems involving the mean curvature operator.

For Fujita type results for degenerate parabolic inequalities and for systems of quasilinear parabolic inequalities  we refer  the reader to \cite{MP_JEv} and to \cite{MP}. In particular, the Fujita exponent for the $m$-Laplacian parabolic case
\begin{equation}\label{fujita661}
\begin{cases}
&u_t - \Delta _mu \ge u^{q} \,\,\, \text{in }\mathbb{R}^{N} \times (0, \infty),\\
&u(x,0) = u_0(x)\ge0 \,\, \text{in } \mathbb{R}^{N},\\
\end{cases}
\end{equation}
is obtained in  \cite{MP}  where it is proved that \eqref{fujita661} has non nonnegative nontrivial solutions for
$$\max\{1,m-1\}<q\le m-1+\frac mN.$$ This last condition forces $m>2N/(N+1)$.
For further generalizations we refer the reader to
\cite{zhou} for inhomogeneous pseudoparabolic equations, to \cite{KK} for  parabolic problems in halfspace, or to \cite{FL} for parabolic inequalities with weights and  nonlinearites depending on the gradient.

Nonlocal models describe many natural phenomena, such as the non-Newton flux in the mechanics of fluid, population of biological species and filtration.
Concerning a nonlocal sources problem, Galaktionov and Levine \cite{GL} investigated positive solutions of a Cauchy problem for the following semilinear parabolic equation with weighted nonlocal sources
\begin{equation}\label{Gal_Lev}
\begin{cases}
&u_t=\Delta u^m +\biggl( \int_{\mathbb R^N}K(x)u^p(x,t)dx\biggr)^{(r-1)/p}u^q\,\,\, \text{in }\mathbb{R}^{N} \times (0, \infty),\\
&u(x,0) = u_0(x)\ge0 \,\, \text{in } \mathbb{R}^{N},\\
\end{cases}
\end{equation}
where $p, q, r\ge1$, $m>1$ and $K$ is a  function not necessarily in $L^1(\mathbb R^N)$. Other semilinear parabolic problems with nonlocal terms are discussed in \cite[Part V]{QS}.

This paper is a first attempt in solving quasilinear parabolic problems of type \eqref{main} and \eqref{main2} which feature  nonlocal terms defined by the convolution operation. Such convolution terms are used to model various quantities in gravity and quantum physics. For instance, the equation
\begin{equation}\label{hartree}
i\psi_t-\Delta \psi=(|x|^{\alpha-N}\ast \psi^2)\psi\quad\mbox{in }{\mathbb R}^N, \alpha\in (0, N), N\geq 1,
\end{equation}
was introduced in 1928 by Hartree \cite{H1,H2,H3}, shortly after the publication of the Schr\"odinger equation, in order to study the  non-relativistic atoms, using the concept of self-consistency. The stationary case of \eqref{hartree} for $N = 3$ and $\alpha = 2$ is known in the literature
as the {\it Choquard} (or {\it Choquard-Pekar}) equation and was introduced in \cite{P} as a model in quantum theory (see also \cite{MS} for a mathematical account on this equation). Stationary singular solutions of the Choquard inequality are discussed in \cite{GT} while quasilinear elliptic equations and inequalities with convolution terms are considered in \cite{FG, GKS2020, GKS2021, GMM2021}.

Our main result concerning \eqref{main} is the following.

\begin{theorem}\label{thmain}
Assume
\begin{equation}\label{cond_beta_m} 0<\beta<m/2\qquad\mbox{ and } \qquad m >\frac{2N+1}{N+1}\bigl(\in (1,2)\bigr).\end{equation}
If
\begin{equation}\label{cond1}
2\max\{1,m-1\}<p+q\le  m-1+\frac{N-\beta+m}{N+\beta} ,
\end{equation}
then problem \eqref{main} does not have nonnegative nontrivial solutions  $u\in \mathcal{S}$.
\end{theorem}

In particular, condition \eqref{cond1} yields the following upper bounds for $\beta$
\begin{equation}\label{cond2_beta}
\beta<
\begin{cases}
\displaystyle 1-\frac{m-2}{m}N &\quad\mbox{ if } 2<m<\frac{2N}{N-1},\\[0.1in]
\displaystyle \frac{m-(2-m)N}{4-m} &\quad\mbox{ if }\frac{2N+1}{N+1}<m<2.
\end{cases}
\end{equation}

The proof of Theorem \ref{thmain} is carried out through nonlinear capacity estimates specifically adapted to the nonlocal setting of our problem and the weak $m$-coercivity of the differential operator $L_{\mathcal{A}}$. Precisely, we derive integral estimates in time for the new quantity
$$
J(t)=\int_{\mathbb R^N}u^\ell (x,t)\psi^k(x, t)dx\, , \,\,\, t\geq 0,
$$
where $\ell=(p+q)/2$ and $\psi$ is a carefully selected test function  with compact support (see \eqref{testfunction}). As corollaries of our main results we obtain:

\begin{corollary} Assume
$$
0<\beta<1 \quad\text{and}\quad
2<p+q\le 2+\frac{2(1-\beta)}{N+\beta}.
$$
Then problem
$$\begin{cases}
&\displaystyle u_t-\mbox{\rm div} \biggl(\frac{\nabla u}{\sqrt{1+|\nabla u|^2}}\biggr)\geq  (K\ast u^p)u^q \quad\mbox{ in } \mathbb R^N \times \mathbb (0,\infty),\\[0.2in]
&u(x,0) = u_0(x)\ge0 \,\, \text{ in } \mathbb R^N,\end{cases}
$$
does not have nonnegative nontrivial solutions.
\end{corollary}

\begin{corollary}
Assume
$$0<\beta<m/2 \quad \mbox{ and }\quad m >\frac{2N+1}{N+1}\bigl(\in (1,2)\bigr).$$
If  \eqref{cond1} holds,
then the problem
$$\begin{cases}
&u_t-\Delta_m u\geq  (K\ast u^p)u^q \quad\mbox{ in } \mathbb R^N \times \mathbb (0,\infty),\\
&u(x,0) = u_0(x)\ge0 \,\, \text{ in } \mathbb R^N,\end{cases}
$$
does not have nonnegative nontrivial solutions.
\end{corollary}

The critical Fujita exponent also depends on the asymptotic behavior at infinity of the initial function $u_0(x)$. The result below is an extension of Lemma 26.2 in \cite{MP} to the setting of nonlocal problems.

\begin{corollary}\label{cornu} Assume that $u_0\in L^1_{loc}(\mathbb{R}^N)$ and for all $R>0$ large enough $u_0$ satisfies
\begin{equation}\label{u_0_nu}
\int_{B_R}u_0(x)dx\ge c R^\nu\,,
\end{equation}
for some positive constant $c$ and an exponent $0\le \nu< N+\beta$, $0<\beta<(m+\nu)/2$.
If
 \begin{equation}\label{critic_nu}
2\max\{1,m-1\}<p+q<m-1+\frac{N-\beta+m}{N+\beta-\nu},
\end{equation}
then problem  \eqref{main} does not have nonnegative nontrivial solutions $u\in \mathcal{S}$.
\end{corollary}
In this case, condition \eqref{critic_nu} yields the following upper bound for $\beta$
\begin{equation}\label{cond2_beta_nu}
\beta<
\begin{cases}
\displaystyle 1-\frac{(m-2)N-(m-1)\nu}{m} &\quad\mbox{ if } 2<m<\frac{2N}{N-1},\\[0.1in]
\displaystyle \frac{m-(2-m)N+(3-m)\nu}{4-m} &\quad\mbox{ if }\frac{2N+1}{N+1}<m<2.
\end{cases}
\end{equation}
In particular, when $m>2$, if we restrict the range of $\nu$ to the set
$$\frac{N(m-2)-m}{m-1}<\nu<N+\beta,$$
then the upper bound $2N/(N-1)$  on  $m$ in \eqref{cond2_beta_nu} can be removed.

Next, we turn to the study of \eqref{main2}. In this setting we derive our a priori estimates employing similar arguments to those we used in the proof of Theorem \ref{thmain}. Unlike the approach for \eqref{main} where the  a priori estimates are obtained by choosing two classes of test functions in \eqref{weak_form2}, for the counterpart problem \eqref{main2} we can only use a single class of test functions (this was already emphasised in \cite[Section 26]{MPbook} for local coercive parabolic problems). We shall overcome this fact by imposing a higher  locally integrability on the solution. The precise solution space will be defined in what follows.

Assume
\begin{equation}\label{new_beta_coerc}0<\beta<m/2\quad\text{and}\quad p+q>2\max\{1,m-1\}.
\end{equation}
Let $d>0$ be such that
\begin{equation}\label{eqd}
2<\Gamma(d):=\max\Big\{\frac{p+q+d}{d+1}\, , \, \frac{p+q+d}{d+m-1}  \Big\}<\frac{2N+m}{N+\beta}.
\end{equation}
Such a value $d>0$ always exists since $\Gamma$ is decreasing as a function of $d$ and
$$
\Gamma(0)>2 \quad\mbox{ and }\quad \lim_{d\to \infty}\Gamma(d)=1<\frac{2N+m}{N+\beta}.
$$

Our main result concerning the problem \eqref{main2} is the following.

\begin{theorem}\label{thmain2}
If \eqref{new_beta_coerc} and \eqref{eqd} hold,
then problem \eqref{main2} does not have nonnegative nontrivial solutions
 $u$ in the class
$$
\mathcal{S}\cap \Big\{(K\ast u^p)u^{q+d} \in  L^{1}_{loc}(\mathbb R^N\times [0, \infty)) \Big\}.
$$
\end{theorem}
In particular, if \eqref{new_beta_coerc} holds, then \eqref{main2} has no nonnegative nontrivial solutions $u\in \mathcal{S}\cap C(\mathbb R^N\times (0, \infty))$.

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

Suppose by contradiction that \eqref{main} admits a nonnegative nontrivial solution
$u\in \mathcal{S}$. We start with the following result which provides extra local integrability of $u$.

\begin{lemma}\label{lp0}
Let $u\in \mathcal{S}$ be a nonnegative solution of \eqref{main}. Then,
\begin{equation}\label{uLp+q_2}
u^{(p+q)/2}\in L^1_{loc}(\mathbb{R}^N\times [0,\infty)).
\end{equation}
\end{lemma}
\begin{proof} Let $R>\rho$ be large, where $\rho>0$ appears in \eqref{la}. For $x\in B_R(0)$, using \eqref{la} we have
\begin{equation}\label{KKK1}
\begin{aligned}
(K\ast u^p)(x,t)& \geq \int_{\mathbb B_{R}} K(x-y)u^p(y,t) dy\\
&=\int_{|x-y|\leq \rho} K(x-y)u^p(y,t))dy+\int_{\substack{|x-y|>\rho\\  |x|,\, |y|<R}} K(x-y)u^p(y,t)dy\\
&\geq \inf_{z\in B_1(0)} K(z) \int_{|x-y|\leq \rho} u^p(y,t)dy+c\int_{\substack{|x-y|>\rho \\  |x|,\, |y|<R}}  |x-y|^{-\beta} u^p(y,t) dy\\
&\geq \inf_{z\in B_1(0)} K(z)  \int_{|x-y|\leq \rho} u^p(y,t) dy+c(2R)^{-\beta}
\int_{|x-y|>\rho ,\, |y|<R}  u^p(y,t) dy\\
&\geq CR^{-\beta}\left\{ \int_{|x-y|\leq \rho\, , \,  |y|<R} u^p(y,t)dy+\int_{|x-y|>\rho ,\, |y|<R}  u^p(y,t) dy \right\} \\
&\geq CR^{-\beta}\int_{B_R(0) } u^p(y,t) dy,
\end{aligned}
\end{equation}
provided $R>\rho$ is large enough.

Next, using \eqref{KKK1} and H\"older's inequality we find
$$
\begin{aligned}
\infty&>\int_0^T\!\!\int_{B_{R}(0)} \big(K\ast u\big)(x,t)u^q(x,t) dxdt \\
&\geq CR^{-\beta}\int_0^T\Big(\int_{B_R(0)} u^p(x) dx\Big)   \Big(\int_{B_R(0)} u^q(x) dx\Big) dt \qquad \mbox{  (by \eqref{KKK1})}\\
&\geq CR^{-\beta} \int_0^T \Big(\int_{B_R(0)} u^{(p+q)/2}(x) dx\Big)^2 dt \qquad \mbox{  (by H\"older's inequality on $B_R(0)$)}\\
&\geq \frac{CR^{-\beta}}{T} \Big( \int_0^T \int_{B_R(0)} u^{(p+q)/2}(x) dx \, dt \Big)^2 \qquad \mbox{  (by H\"older's inequality on $[0,T]$)}\\
\end{aligned}
$$
which shows that $u^{(p+q)/2}\in L^1(B_R(0)\times [0,T])$ and concludes our proof.
\end{proof}

From \eqref{cond1} we have $(p+q)/2>\max\{1, m-1\}$. We may thus choose $d>0$ small enough such that $(p+q)/2>(m-1)(d+1)$. Using  H\"older's inequality we deduce
\begin{equation}\label{Summ_with_d}
u^{1-d}\, , \, u^{(m-1)(1-d)}\,,\, u^{(m-1)(d+1)}\in L^1_{loc}(\mathbb{R}^N\times [0,\infty)).
\end{equation}
This will ensure that all integrals in this section are finite.

The proof of Theorem \ref{thmain} will be achieved along a several lemmas. First we want to precise the choice of the test functions in \eqref{weak_form2}. Take a standard cut off function $\xi \in C^{1}[0,\infty)$ such that
\begin{itemize}
\item $\xi = 1$ in $(0,1)$, $\xi= 0$  in $(2,\infty)$;
\item  $0\le \xi\le 1$ and $|\xi'| \leq C$ in $[0, \infty)$, for some $C>0$.
\end{itemize}

Now take $R > 0$ and consider the functions
\begin{equation} \label{cutoff}
\chi_R (x) = \xi \biggl(\frac{|x|}{R}\biggr) ,\quad \eta_{R}(t ) = \xi\biggl(\frac{t}{R^{\gamma}} \biggr),
\end{equation}
with $\gamma\ge1$. Clearly $\chi_R (x)=1$ in $B_R(0)$, where $B_R(0)$ denotes  the open ball in $\mathbb{R}^N$, centered at the origin and having radius $R > 0$.

Finally define, for all $R > 0$, the nonnegative cut off function $\psi:\mathbb R^N\times [0, \infty)\to [0, \infty)$, given by
\begin{equation} \label{testfunction}
\psi(x,t) = \chi_R(x) \, \eta_{R}(t).
\end{equation}

Note that,  $\psi$ is the restriction to  $\mathbb R^N\times [0, \infty)$ of a
$C^1_c(\mathbb R^N\times \mathbb R)$ function, and,  as in \cite[Lemma 3.1]{FL}, by the shape of $\psi$ in \eqref {testfunction} the following inequalities hold for $\varsigma >1$, $\gamma\ge1$, $k$ and $R$ sufficiently large, say $k>\varsigma$,
\begin{equation}\label{PSI}
\int_0^\infty\!\int_{\mathbb R^N}{\psi^{k-\varsigma}}{|\nabla\psi|^{\varsigma}}dxdt\le c R^{-\varsigma +N +\gamma},\end{equation}
\begin{equation}\label{PSII} \int_0^\infty\!\int_{\mathbb R^N}{\psi^{k-\varsigma}}{|\psi_t|^{\varsigma}} dxdt \le c R^{-\gamma\varsigma +N +\gamma},
\end{equation}
where $c>0$ is a constant.
\begin{lemma}\label{lp1}
Let $u\in \mathcal{S}$ be a nonnegative solution of \eqref{main} and let $\psi$ be defined by \eqref{testfunction}. Then,
\begin{equation}\label{ineq2.0}
\begin{aligned}
\int_0^\infty\int_{\mathbb{R}^N} & (K\ast u^p)u^{q-d} \psi^{k} \, dx \, dt+
\int_0^\infty\!\!\int_{\mathbb R^N}  u^{- d - 1} \psi^{k} |\mathcal{A}(x,u,\nabla u)|^{m'}\, dx \,  dt  \\
&\leq   c_1\int_0^\infty\!\!\int_{\mathbb R^N}  u^{1 - d}\psi^{k - 1}|\psi_t| \, dx  dt +
c_2 \int_0^\infty\!\!\int_{\mathbb R^N}  u^{ m- d- 1}\psi^{k - m} |\nabla \psi|^m  dx \, dt,
\end{aligned}\end{equation}
for some constants $c_1,c_2>0$ and $d\in(0,1)$.
\end{lemma}

\begin{proof} Let $\epsilon > 0$ be sufficiently small and let $\{\xi_{\epsilon}\}_{\epsilon>0}$ be a standard family of mollifiers. For $\tau > 0$ we define
$$
\tilde{u}_\epsilon(x,t)  :=  \tau + \int_{\mathbb R^N} \, \xi_{\varepsilon}(x-y,t) u(y,t) \, dy \quad \text{and} \qquad u_{\tau}(x,t) := \tau + u(x,t)
$$
for $(x,t)\in {\mathbb R}^N\times (0, \infty)$. Clearly $\tilde{u}_{\epsilon}, u_{\tau} \geq \tau > 0$. In particular, since
$u \in L^{1}_{\rm loc} (\mathbb R^N\times (0, \infty))$, we have $\tilde{u}_{\epsilon} \in C^1(\mathbb R^N\times (0, \infty))$ so that the function
$
\varphi(x,t) = \tilde{u}_{\epsilon}^{-d} \; \psi^{k}(x,t)\geq 0,
$
with $d\in(0,1)$ sufficiently small so that \eqref{Summ_with_d} holds, $k$ positive to be chosen and  $\psi(x,t)$ defined in $\eqref{testfunction}$,
can be used as test function in the weak formulation of $\eqref{main}$, given by $\eqref{weak_form2}$, so that
$$\begin{aligned}
\int_0^\infty\int_{\mathbb{R}^N} (K\ast u^p)&u^q\tilde{u}_\epsilon^{-d} \psi^{k} \, dx \, dt +d \int_0^\infty\int_{\mathbb{R}^N} \psi^{k}\tilde{u}_{\epsilon}^{-d-1}\,\mathcal{A}(x,u,\nabla u) \cdot \nabla\tilde{u}_{\epsilon} \, dx \, dt \\& \le \int_0^\infty\int_{\mathbb{R}^N}  u_t \;
\tilde{u}_{\epsilon}^{-d} \; \psi^{k} \, dx \, dt
 + k  \int_0^\infty\int_{\mathbb{R}^N} \psi^{k-1} \tilde{u}_{\epsilon}^{-d} \mathcal{A}(x,u,\nabla u) \cdot  \nabla\psi \, dx \, dt. \\
\end{aligned}$$

Since $\tilde{u}_{\epsilon} \to u_{\tau}$ in $W^1_{\rm loc}(\mathbb{R}^N)$ as $\varepsilon \to 0$,
using Lebesgue dominated convergence theorem, and being $\nabla u_\tau=\nabla u$, we arrive at
\begin{equation}\begin{aligned}
\int_0^\infty\int_{\mathbb{R}^N} &(K\ast u^p)u^q u_\tau^{-d} \psi^{k} \, dx \, dt +d\int_0^\infty\int_{\mathbb{R}^N} \psi^{k} u_\tau^{-d-1}\,\mathcal{A}(x,u,\nabla u) \cdot \nabla u \, dx \, dt
\\& \le \int_0^\infty\int_{\mathbb{R}^N}  u_t \,
u_\tau^{-d} \psi^{k} \, dx \, dt
 + k  \int_0^\infty\int_{\mathbb{R}^N}  \psi^{k-1}u_\tau^{-d} |\mathcal{A}(x,u,\nabla u)| |\nabla\psi| \, dx \, dt.
\end{aligned}\end{equation}
Using the weak $m$-coerciveness of $\mathcal{A}$, from \eqref{weak_m} we deduce
\begin{equation}\label{ineq0}\begin{aligned}
\int_0^\infty\int_{\mathbb{R}^N} &(K\ast u^p)u^q u_\tau^{-d} \psi^{k} \, dx \, dt +
dc_0  \int_0^\infty\!\!\int_{\mathbb{R}^N} \psi^{k} u_\tau^{-d-1}\,|\mathcal{A}(x,u,\nabla u)|^{m'} \, dx \, dt
\\ & \le \int_0^\infty \!\!\int_{\mathbb{R}^N}  u_t \,
u_\tau^{-d} \psi^{k} \, dx \, dt
 + k  \int_0^\infty\!\!\int_{\mathbb{R}^N} \psi^{k-1}u_\tau^{-d} |\mathcal{A}(x,u,\nabla u)| |\nabla\psi| \, dx \, dt.
\end{aligned}\end{equation}
Now consider the first term on the right-hand side  of \eqref{ineq0} and we claim that
\begin{equation} \label{part1}
\int_0^\infty\!\!\int_{\mathbb{R}^N}  u_t \, u_{\tau}^{-d} \, \psi^k \, dx \, dt \leq  \frac{k}{1 - d} \int_0^\infty\!\!\int_{\mathbb{R}^N} \, u_{\tau}^{1 - d}
 \psi^{k - 1}  |\psi_t| \, dx \, dt.
\end{equation}
Indeed, since $(u_{\tau})_t = u_t$, by definition of $u_\tau$, we have
$$
\begin{aligned}
\int_0^\infty\int_{\mathbb{R}^N} u_t \, u_{\tau}^{-d} \psi^k \, dx \, dt &=\int_0^\infty\int_{\mathbb{R}^N} (u_{\tau})_t \; u_{\tau}^{-d} \; \psi^{k} \, dx \, dt = \int_{S}
\frac 1{1- d}\bigl( {u_{\tau}^{1 - d}} \bigr)_t \psi^{k} \, dx \, dt \\
&= \frac{1}{1 - d}  \int_0^\infty\int_{\mathbb{R}^N} \bigg[ (u_{\tau}^{1 - d} \; \psi^k)_t - u_{\tau}^{1 - d}
(\psi^k)_t \bigg] \, dx \, dt  \\
&= - \, \frac{1}{1 - d}\int_{\mathbb R^N} \, u_{\tau}^{1 - d}(x,0) \psi^k(x,0) \, dx -
\frac{k}{1 - d} \int_0^\infty\int_{\mathbb{R}^N} \, u_{\tau}^{1 - d} \; \psi^{k - 1}  \psi_t \, dx \, dt \\
\end{aligned}
$$
where the last equality is due to $\psi \in C^{1}_{c}(\mathbb R^N\times \mathbb [0,\infty))$, hence $ \lim_{t \to \infty} \psi(x,t) = 0.$ Consequently, $\eqref{part1}$ follows immediately from
$u_{\tau}(x,0) =  u_0(x) + \tau > 0$, $\psi\ge0$ and $0<d < 1$.

Further, by Young inequality (see also \cite[Lemma 4.1]{FL}) we have
$$
u_\tau^{-d} \psi^{k-1} |\mathcal A(x,u,\nabla u)|| \nabla \psi| \leq
\frac{dc_0}{2 k} u_\tau^{-d-1}\psi^{k}
|\mathcal A(x,u,\nabla u)|^{m'}  + Cu_\tau^{ m-d - 1}\psi^{k - m} |\nabla \psi|^m.
$$
Thanks to property (i) in the definition of $\mathcal S$ and the fact that $u_\tau\ge\tau>0$ and $0\le\psi\le1$, we have
$$ u_\tau^{-d-1}\ \psi^{k}|\mathcal{A}(x,u,\nabla u)|^{m'}
\in L^1_{loc}(\mathbb{R}^N\times [0,\infty)).$$
Thus, a combination of   \eqref{ineq0} and \eqref{part1} yields
$$
\begin{aligned}
\int_0^\infty\!\!\int_{\mathbb{R}^N} (K\ast u^p)u^q & u_\tau^{-d} \psi^{k} \, dx \, dt +
\int_0^\infty\!\!\int_{\mathbb R^N}  u_{\tau}^{- d - 1} \psi^{k} |\mathcal{A}(x,u,\nabla u)|^{m'}\, dx \,  dt  \\
&\leq   c_1\int_0^\infty\!\!\int_{\mathbb R^N}  u_{\tau}^{1 - d}\psi^{k - 1}|\psi_t| \, dx  dt + c_2 \int_0^\infty\!\!\int_{\mathbb R^N}  u_\tau^{ m- d- 1}\psi^{k - m} |\nabla \psi|^m  dx \, dt,
\end{aligned}$$
with $c_1,c_2>0$. Since all the exponents of $u$ on the right hand side are positive and \eqref{Summ_with_d}, we let  $\tau\to0$ and apply Fatou and Lebesgue theorems, to obtain \eqref{ineq2.0}.
\end{proof}

\begin{lemma}\label{lp2}

Let $\ell>1$ and $u\geq 0$ be a solution of \eqref{main} such that $u^\ell\in L^1_{loc}(\mathbb{R}^N\times [0, \infty)).$

Define
\begin{equation}\label{Jdef}
J(t)=\int_{\mathbb R^N}u^\ell (x,t)\psi^k(x, t)dx,
\end{equation}
where $\psi$ is given by \eqref{testfunction}. Then,
\begin{equation}\label{ineqE}\begin{aligned}
\int_0^\infty\!\!\int_{\mathbb R^N} (K\ast u^p)u^q \psi^k \, dx \, dt \le & \; c_1\biggl(\int_0^{2R^\gamma} J(t)\biggr)^{1/\ell} \cdot
R^{\frac{N+\gamma}{\ell'}-\gamma}\\
&\;+c_2\biggl(\int_0^{2R^\gamma} J(t)\biggr)^{\frac{2(m-1)}{m\ell}}\cdot
R^{(N+\gamma)\bigl(1-\frac{2}{m'\ell}\bigr) -1-\frac{\gamma}{m'}}\\
&\;+c_3 \biggl(\int_0^{2R^\gamma} J(t)\biggr)^{\frac{m-1}{\ell}}
\cdot R^{(N+\gamma)\bigl(1-\frac{m-1}{\ell}\bigr)-m},
\end{aligned}
\end{equation}
for some constants $c_1,c_2,c_3>0$.
\end{lemma}

\begin{proof} We first choose $\varphi=\psi^k$ in the weak formulation \eqref{weak_form1}, with  $\psi$ given by \eqref{testfunction}. Since $u_0, \varphi\ge0$,  we find
\begin{equation}\label{weak_k}
 \int_0^\infty\!\!\int_{\mathbb R^N} (K\ast u^p)u^q \psi^k \, dx \, dt \leq  - \!\int_0^\infty\!\!\int_{\mathbb R^N} u \, (\psi^k)_{t} \, dx \, dt
 +\!\!  \iint\limits_{{\rm supp}\,(\nabla \psi)}  |\mathcal{A}(x,u,\nabla u)||\nabla\psi^k| dx \, dt.
\end{equation}
From now on, the constant $c_1,\, c_2$ and $c_3$ will assume  different values.

By   H\"older inequality  we obtain
\begin{equation}\label{weak_kk}
\begin{aligned}
\int_0^\infty\!\!\int_{\mathbb R^N} (K\ast u^p)u^q \psi^k \, dx \, dt \leq &
\, C \biggl( \iint\limits_{{\rm supp}\,  (\psi_t)}  u^\ell\psi^kdx dt\biggr)^{1/\ell}\cdot \biggl(\iint\limits_{{\rm supp}\,  (\psi_t) }
\psi^{k-\ell'}|\psi_t|^{\ell'}\biggr)^{1/\ell'} \\
& +\!\! \iint\limits_{{\rm supp}(\nabla \psi)} |\mathcal{A}(x,u,\nabla u)||\nabla\psi^k| dx \, dt.
\end{aligned}
\end{equation}
We next estimate the last integral in \eqref{weak_kk}. First, by H\"older's inequality we have
$$
\begin{aligned}
\iint\limits_{{\rm supp}(\nabla \psi)} |\mathcal{A}(x,u,\nabla u)||\nabla\psi^k| dx \, dt\leq &
k \biggl(\int_0^\infty\!\!\int_{\mathbb R^N}u^{-d-1}\psi^k|\mathcal{A}(x,u,\nabla u)|^{m'} dx\, dt \biggr)^{1/m'}\\
& \cdot
\biggl(\int_0^\infty\!\!\int_{\mathbb R^N}u^{(d+1)(m-1)}\psi^{k-m}|\nabla \psi |^m dx\, dt\biggr)^{1/m}.
\end{aligned}
$$
Next, we use the estimate \eqref{ineq2.0} from Lemma \ref{lp1} and the standard inequality $(a+b)^r\le c(a^r+b^r)$ for $a,b,r>0$. We deduce

\begin{equation}\label{ineq0_Hold}
\begin{aligned}\qquad\quad\int_0^\infty\!\!\int_{\mathbb R^N} |\mathcal{A}(x,u,\nabla u)|&|\nabla\psi^k| dx \, dt\le
c_1 \biggl(\int_0^\infty\!\!\int_{\mathbb R^N}  u^{1 - d}\psi^{k - 1}|\psi_t| \, dx  dt \biggr)^{1/m'}\\
&\,\,\cdot
\biggl(\int_0^\infty\!\!\int_{\mathbb R^N}u^{(d+1)(m-1)}\psi^{k-m}|\nabla \psi |^m dx\, dt\biggr)^{1/m}
\\&+ c_2\biggl(\int_0^\infty\!\!\int_{\mathbb R^N}  u^{ m- d- 1}\psi^{k - m} |\nabla \psi|^m  dx \, dt\biggr)^{1/m'}\\
&\,\,\cdot
\biggl(\int_0^\infty\!\!\int_{\mathbb R^N}u^{(d+1)(m-1)}\psi^{k-m}|\nabla \psi |^m dx\, dt\biggr)^{1/m}.
\end{aligned}\end{equation}
We now use H\"older inequality in all the factors of the right hand side, so that
\begin{equation}\label{ineqA}
\begin{aligned}
\int_0^\infty\!\!\int_{\mathbb R^N}  u^{1 - d}\psi^{k - 1}|\psi_t| \, dx  dt
&\leq \Big(\iint\limits_{{\rm supp}\,  (\psi_t)}  u^\ell \psi^k\Big)^{1/\sigma}
\Big(\iint\limits_{{\rm supp}\,  (\psi_t)} \psi^{k-\sigma'}|\psi_t|^{\sigma'}\Big)^{1/\sigma'},\\
\end{aligned}
\end{equation}
where
$$\sigma=\frac{\ell}{1-d},\qquad \sigma'=\frac \ell{\ell+d-1};$$
\begin{equation}\label{ineqB}
\begin{aligned}
\int_0^\infty\!\!\int_{\mathbb R^N}u^{(d+1)(m-1)}&\psi^{k-m}|\nabla \psi |^m dx\, dt\\
&\leq \Big( \iint\limits_{{\rm supp}\,  (\nabla \psi)}    u^\ell \psi^k\Big)^{1/\eta}
\Big(\int_0^\infty\!\!\int_{\mathbb R^N}   \psi^{k-m\eta'}|\nabla \psi|^{m\eta'}\Big)^{1/\eta'},
\end{aligned}
\end{equation}
 where
\begin{equation}\label{eta}
\eta=\frac{\ell}{(d+1)(m-1)},\quad \eta'=\frac{\ell}{\ell-m+1-d(m-1)};
\end{equation}
\begin{equation}\label{ineqC}
\begin{aligned}
\int_0^\infty\!\!\int_{\mathbb R^N}  u^{m-1-d} \psi^{k-m}|\nabla \psi |^m
&\leq \Big(\iint\limits_{{\rm supp}\,  (\nabla \psi)}  u^\ell \psi^k\Big)^{1/\theta}
\Big(\int_0^\infty\!\!\int_{\mathbb R^N}   \psi^{k-m\theta'}|\nabla \psi|^{m\theta'}\Big)^{1/\theta'},\\
\end{aligned}
\end{equation}
where
\begin{equation}\label{theta}
\theta=\frac{\ell}{m-d-1},\quad \theta'=\frac{\ell}{\ell-m+d+1}.
\end{equation}
We next replace \eqref{ineqA}, \eqref{ineqB} and \eqref{ineqC} in \eqref{ineq0_Hold}. With $J(t)$ defined in \eqref{Jdef} we obtain

$$
\begin{aligned}&\int_0^\infty\!\!\int_{\mathbb R^N} |\mathcal{A}(x,u,\nabla u)||\nabla\psi^k| dx \, dt\le\\
&\;\; c_1 \biggl(\int_0^\infty J(t)\biggr)^{\frac1{m'\sigma}+\frac1{m\eta}}
\Big(\int_0^\infty\!\!\int_{\mathbb R^N}  \psi^{k-\sigma'}|\psi_t|^{\sigma'}\Big)^{1/m'\sigma'}
\Big(\int_{\mathbb R^N}  \psi^{k-m\eta'}|\nabla \psi|^{m\eta'}\Big)^{1/m\eta'}\\
&+c_2 \biggl(\int_0^\infty J(t)\biggr)^{\frac1{m'\theta}+\frac1{m\eta}}
\Big(\int_0^\infty\!\!\int_{\mathbb R^N}  \psi^{k-m\theta'}|\nabla \psi|^{m\theta'}\Big)^{1/m'\theta'}
\Big(\int_0^\infty\!\!\int_{\mathbb R^N}  \psi^{k-m\eta'}|\nabla \psi|^{m\eta'}\Big)^{1/m\eta'}.
\end{aligned}
$$
Inserting the above inequality in \eqref{weak_kk}  and using \eqref{PSI} and \eqref{PSII} together with the fact that $J(t)=0$ for $t\geq 2R^\gamma$ {from the definition of the test function $\psi$},   we find
\begin{equation}\label{ineqD}
\begin{aligned}
\int_0^\infty\!\!\int_{\mathbb R^N} (K\ast u^p)u^q \psi^k \, dx \, dt \le & \; c_1\biggl(\int_0^{2R^\gamma} J(t)\biggr)^{1/\ell} \cdot
R^{\frac{N+\gamma-\gamma\ell'}{\ell'}}\\
&\;+c_2\biggl(\int_0^{2R^\gamma} J(t)\biggr)^{\frac1{m'\sigma}+\frac1{m\eta}}\cdot
R^{\frac{N+\gamma-\gamma\sigma'}{m'\sigma'}+\frac{N+\gamma-m\eta'}{m\eta'}}\\
&\; +c_3 \biggl(\int_0^{2R^\gamma} J(t)\biggr)^{\frac1{m'\theta}+\frac1{m\eta}}
\cdot R^{\frac{N+\gamma-m\theta'}{m'\theta'}+\frac{N+\gamma-m\eta'}{m\eta'}}\\
= & \;c_1\biggl(\int_0^{2R^\gamma} J(t)\biggr)^{1/\ell} \cdot
R^{\frac{N+\gamma-\gamma\ell'}{\ell'}}\\
&\; +c_2\biggl(\int_0^{2R^\gamma} J(t)\biggr)^{\frac1{m'\sigma}+\frac1{m\eta}}\cdot
R^{(N+\gamma)\bigl (\frac 1{m'\sigma'}+\frac1{m\eta'}\bigr) -1-\frac{\gamma}{m'}}\\
&\;+c_3 \biggl(\int_0^{2R^\gamma} J(t)\biggr)^{\frac1{m'\theta}+\frac1{m\eta}}
\cdot R^{(N+\gamma)\bigl(\frac1{m'\theta'}+\frac1{m\eta'}\bigr)-m}.
\end{aligned}
\end{equation}
Consequently, being
$$\frac 1{m'\sigma }+\frac 1{m\eta }=\frac{2(m-1)}{m\ell},  \qquad \frac 1{m'\theta } +\frac 1{m\eta }=\frac{m-1}{\ell}, $$
$$\frac1{m'\sigma'}+\frac1{m \eta'}=\frac{m\ell-2(m-1)}{m\ell},\qquad  \frac1{m'\theta'}+\frac1{m \eta'}=\frac{\ell-m+1}{\ell}, $$
 using \eqref{ineqD} and the values of the parameters involved, we get the estimate \eqref{ineqE}.
\end{proof}

\begin{lemma}\label{lp3}
Let $u\geq 0$ be a solution of \eqref{main} and $J$ be given by \eqref{Jdef} with $\ell=(p+q)/2>1$. Then
\begin{equation}\label{ineq_final}
\int_0^{2R^\gamma} J(t)^2 dt \le c\bigl(R^{\alpha_1}+R^{\alpha_2}+R^{\alpha_3}\bigr),
\end{equation}
where $c>0$ is a constant and
$$
\begin{aligned}
\alpha_1& =\frac{(N+\beta)(p+q)-2N-\gamma}{p+q-1}\,,\\
\alpha_2&=\frac{(N+\beta-1)m(p+q)-4N(m-1)+\gamma[p+q-2(m-1)]}{m(p+q)-2(m-1)}\,,\\
\alpha_3& =\frac{(N+\beta-m)(p+q)-2N(m-1)+\gamma(p+q-m+1)}{p+q-m+1}.
\end{aligned}
$$
\end{lemma}

\begin{proof}
First note that by \eqref{uLp+q_2}, we can choose $\ell=(p+q)/2$ so that the requirement of Lemma \ref{lp2} is satisfied.
Now, observe that supp$\, (\psi^k)\subset B_{2R}(0)\times [0, 2R^\gamma)$.
If $(x,t)\in B_{2R}(0)\times [0, \infty)$, then, in the same way as we estimated \eqref{KKK1} we find
\begin{equation}\label{KKK2}
(K\ast u^p)(x,t)
\geq CR^{-\beta}\int_{B_{2R}(0) } u^p(y,t)dy
\geq CR^{-\beta}\int_{\mathbb R^N } u^p(y,t)\psi^k(y,t)dy,
\end{equation}
provided $R>\rho$ is large enough.

Furthermore, for $\ell=(p+q)/2>1$, by H\"older's inequality we have
$$
\begin{aligned}
\biggl(\iint_{\mathbb R^N\times \mathbb R^N}u^p(y,t)\psi^k(y,t) & u^{q}(x,t) \psi^{k}(x,t)  \, dx \, dy\biggr) ^2\\
=\; &\biggl(\iint_{\mathbb R^N\times \mathbb R^N}u^p(y,t)\psi^k(y,t) u^{q}(x,t) \psi^{k}(x,t) \, dx \, dy\biggr)\\
\; &\cdot \biggl(\iint_{\mathbb R^N\times \mathbb R^N}u^p(x,t)\psi^k(x,t) u^{q}(y,t) \psi^{k}(y,t) \, dx \, dy\biggr) \\
\geq \; &
\biggl(\iint_{\mathbb R^N\times \mathbb R^N}u^{\frac{p+q}{2}} (x,t)u^{\frac{p+q}{2}} (y,t)\psi^k(x,t)  \psi^{k}(y,t) \, dx \, dy\biggr)^2\\
=\; & \biggl(\int_{\mathbb R^N}u^{\ell} (x,t)\psi^k(x,t)  \, dx \biggr)^4=J(t)^4,
\end{aligned}
$$
where $J$ is given by \eqref{Jdef} with $\ell=(p+q)/2>1$. Hence,
$$
\iint_{\mathbb R^N\times \mathbb R^N}u^p(y,t)\psi^k(y,t)  u^{q}(x,t) \psi^{k}(x,t)  \, dx \, dy\geq J^2(t).
$$
Hence, using the above estimate and \eqref{KKK2} we deduce
\begin{equation}\label{JJ1}\begin{aligned}
& \int_{\mathbb R^N} (K\ast u^p)u^q \psi^k \, dx \, dt\\
&\qquad\quad\ge C R^{-\beta}\iint_{\mathbb R^N\times \mathbb R^N}u^p(y,t)\psi^k(y,t)  u^{q}(x,t) \psi^{k}(x,t)  \, dx \, dy\\
&\geq  C R^{-\beta} J(t)^2.
\end{aligned}\end{equation}

As observed before inequality \eqref{ineqD}, we have $J(t)=0$ if $t\geq 2R^{\gamma}$.
Thus, \eqref{JJ1} and inequality \eqref{ineqE} in Lemma \ref{lp2} yields
\begin{equation}\label{JJ2}\begin{aligned}
\int_0^{2R^\gamma} J(t)^2 dt &\le c_1\biggl(\int_0^{2R^\gamma} J(t) dt\biggr)^{1/\ell} \cdot
R^{\frac{N+\gamma}{\ell'}-\gamma+\beta}\\
&\quad+c_2\biggl(\int_0^{2R^\gamma} J(t) dt\biggr)^{\frac{2(m-1)}{m\ell}}\cdot
R^{(N+\gamma)\bigl(1-\frac{2}{m'\ell}\bigr) -1-\frac{\gamma}{m'}+\beta}\\
&\quad+c_3 \biggl(\int_0^{2R^\gamma} J(t) dt\biggr)^{\frac{m-1}{\ell}}
\cdot R^{(N+\gamma)\bigl(1-\frac{m-1}{\ell}\bigr)-m+\beta}\\
\end{aligned}
\end{equation}
In addition, by H\"older's inequality we find
\begin{equation}\label{Holder}
\int_0^{2R^\gamma} J(t) dt\leq \sqrt{2}R^{\gamma/2}\Big(\int_0^{2R^\gamma} J(t)^2 dt\Big)^{1/2},
\end{equation}
so that inequality \eqref{JJ2} implies
 \begin{equation}\label{JJ2_1}\begin{aligned}
\int_0^{2R^\gamma} J(t)^2 dt &\le c_1\biggl(\int_0^{2R^\gamma} J(t)^2 dt\biggr)^{1/2\ell} \cdot
R^{\frac{N+\gamma}{\ell'}-\gamma+\beta+\frac{\gamma}{2\ell}}\\
&\quad+c_2\biggl(\int_0^{2R^\gamma} J(t)^2 dt\biggr)^{\frac{m-1}{m\ell}}\cdot
R^{(N+\gamma)\bigl(1-\frac{2}{m'\ell}\bigr) -1-\frac{\gamma}{m'}+\beta+\frac{\gamma}{m'\ell}}\\
&\quad+c_3 \biggl(\int_0^{2R^\gamma} J(t)^2 dt\biggr)^{\frac{m-1}{2\ell}}
\cdot R^{(N+\gamma)\bigl(1-\frac{m-1}{\ell}\bigr)-m+\beta+\frac{\gamma(m-1)}{2\ell}}.
\end{aligned}
\end{equation}
Now, applying Young inequality we arrive to
\begin{equation}\label{JJ2_2}\begin{aligned}
\int_0^{2R^\gamma} J(t)^2 dt &\le \frac 16\int_0^{2R^\gamma} J(t)^2 dt +c_4
R^{\bigl[\frac{N+\gamma}{\ell'}+\beta-\frac\gamma{(2\ell)'}\bigr](2\ell)'}\\
&\quad+\frac 16\int_0^{2R^\gamma} J(t)^2 dt +c_5
R^{\bigl[(N+\gamma)\bigl(1-\frac{2}{m'\ell}\bigr) -1-\frac{\gamma}{m'\ell'}+\beta\bigr]\bigl(\frac{m\ell}{m-1}\bigr)'}\\
&\quad+\frac 16 \int_0^{2R^\gamma} J(t)^2 dt
+c_6 R^{\bigl[(N+\gamma)\bigl(1-\frac{m-1}{\ell}\bigr)-m+\beta+\frac{\gamma(m-1)}{2\ell}\bigr]\bigl(\frac{2\ell}{m-1}\bigr)'},
\end{aligned}
\end{equation}
which yields \eqref{ineq_final},   being
$$\ell=\frac{p+q}2,\quad \biggl(\frac{m\ell}{m-1}\biggr)'=\frac{m(p+q)}{m(p+q)-2(m-1)}, \quad
\biggl(\frac{2\ell}{m-1}\biggr)'=\frac{p+q}{p+q-m+1}.
$$
\end{proof}

{\it Proof of Theorem \ref{thmain}.} The proof  will be carried out by taking a specific value of
$\gamma\geq 1$ in the definition of the cut off functions in \eqref{cutoff}. To this aim,
we write $\alpha_1$, $\alpha_2$ and $\alpha_3$ so that
$$\alpha_1=\frac{\Upsilon}{p+q-1},\quad \alpha_2=\frac{2(m-1)}{m(p+q)-2(m-1)}\Upsilon, \quad \alpha_3=\frac{m-1}{p+q-m+1}\Upsilon$$
for a certain real constant $\Upsilon$ to be determined.
Since
$$
\begin{aligned}
\alpha_2&=\frac{2(m-1)}{m(p+q)-2(m-1)}\biggl[\frac{N+\beta-1}{2(m-1)}m(p+q)-2N+\gamma\,\frac{p+q-2(m-1)}{2(m-1)}\biggr]\,,\\
\alpha_3 & =\frac{m-1}{p+q-m+1}\biggl[\frac{N+\beta-m}{m-1}(p+q)-2N+\gamma\,\frac{p+q-m+1}{m-1}\biggr]\,,
\end{aligned}
$$
we obtain that
\begin{equation}\label{gamma_iupsilon}\gamma=(m-2)(N+\beta)+m, \quad \Upsilon=(N+\beta)(p+q)-[Nm+\beta(m-2)+m].
\end{equation}
Consequently, for
$$p+q<m-1+\frac{N-\beta+m}{N+\beta}
$$
we obtain $\Upsilon<0$, that is, $\alpha_1,\,\alpha_2, \, \alpha_3<0$. By letting $R\to\infty$ in  \eqref{ineq_final}, it follows that
\begin{equation}\label{esst}
\int_0^{ \infty} J(t)^2 dt=0.
\end{equation}
In turn, using the definition of $J$ in \eqref{Jdef} and the fact that $\psi\equiv1$ in $B_R\times[0, R^\gamma)$,
this easily yields $u\equiv 0$, a contradiction.

If
\begin{equation}\label{critical}
p+q=m-1+\frac{N-\beta+m}{N+\beta},
\end{equation}
then inequality \eqref{ineq_final} for $R\to\infty$ implies $$J\in L^2(0, \infty).$$

Note that
$$
{\rm supp}\; (\nabla \psi)=\bigl(B_{2R}\setminus B_R\bigr)\times (0, 2R^\gamma)\quad \mbox{ and }\quad
{\rm supp}\, (\psi_t)=B_{2R}\times (R^\gamma,2R^\gamma).
$$
Define next
$$
\Theta_R=\max\left\{ \int_{R^\gamma}^{2R^\gamma}\biggl(\int_{B_{2R}}u^\ell\psi^kdx\biggr)^2dt\,,
\int_0^{2R^\gamma} \biggl(\int_{B_{2R}\setminus B_R} u^\ell\psi^kdx\biggr)^2dt
\right\}.
$$
Since $J\in L^2(0, \infty)$, it follows that $\Theta_R\to 0$ as $R\to \infty$.
We now retake the estimate \eqref{weak_kk} and all the calculations that follow up to
\eqref{JJ2_1} to derive
$$
\int_0^{2R^\gamma} J(t)^2 dt \le C\biggl \{\Theta_R^{\frac{1}{2\ell}} +\Theta_R^{\frac{m-1}{m\ell}} + \Theta_R^{\frac{m-1}{\ell}} \biggr\},
$$
being zero all the exponents of $R$ on the right hand side of \eqref{JJ2_1}, by \eqref{critical}
and \eqref{gamma_iupsilon} where $\Upsilon=0$.
Now, letting $R\to \infty$ in the above estimate it follows that $J\equiv 0$ and then $u\equiv 0$ which again contradicts our assumption.
\qed

\noindent{\bf Proof of Corollary \ref{cornu}.} In Lemma  \ref{lp3} we choose the value of $\gamma$ as follows
$$\gamma=\frac{(N-\beta)(m-2)+m(p+q-1)}{p+q-m+1}.$$
Then
$$\alpha_1=\alpha_2=\alpha_3:=\alpha=N+\beta-\frac{N-\beta+m}{p+q-m+1}.$$
Starting with \eqref{weak_form1} and going through the same estimates as above without removing the term $\int_{\mathbb R^N}u_0(x)\varphi(x,0)dx$, the inequality
\eqref{ineq_final} changes to
$$\int_0^{2R^\gamma} J(t)^2 dt \le c_1R^{\alpha}-c_2\int_{B_R}u_0(x)dx,$$
where we have used $\psi\equiv1$ in $B_R$.
Thanks to \eqref{u_0_nu} we arrive to
$$\int_0^{2R^\gamma} J(t)^2 dt \le c\bigl(R^{\alpha}-R^\nu\bigr).$$
Since condition \eqref{critic_nu} forces $\alpha<\nu$, the above estimate yields a contradiction as $R\to\infty$ in both cases $\nu>0$ and $\nu=0$.\qed

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

Suppose by contradiction that \eqref{main2} admits a nonnegative nontrivial solution $u\in \mathcal{S}$ such that
\begin{equation}\label{Kextra}
(K\ast u^p)u^{q+d} \in  L^{1}_{loc}(\mathbb R^N\times [0, \infty)).
\end{equation}
In the same way as in Lemma \ref{lp0} (we only have to replace the exponent $q$ by $q+d$) we deduce
$$
u^{(p+q+d)/2}\in L^1_{loc}(\mathbb{R}^N\times [0,\infty)).
$$
From \eqref{eqd} we have
$$
\frac{p+q+d}{2}>\max\{d+1\, , \, d+m-1\}.
$$
Thus, by H\"older's inequality we deduce
$$
u^{d+m-1}\,,\, u^{d+1}\in L^1_{loc}(\mathbb{R}^N\times [0,\infty)).
$$
This will ensure that all integrals in this section are finite.

We start with the following result which is a counterpart of Lemma \ref{lp1}.

\begin{lemma}\label{lpp1}
 Let $u\in \mathcal{S}$ be a nonnegative solution of \eqref{main2} satisfying \eqref{Kextra} and let $\psi$ be defined by \eqref{testfunction}. Then,
\begin{equation}\label{ineqq2.0}
\begin{aligned}
\int_0^\infty\int_{\mathbb{R}^N} & (K\ast u^p)u^{q+d} \psi^{k} \, dx \, dt+
\int_0^\infty\!\!\int_{\mathbb R^N}  u^{d - 1} \psi^{k} |\mathcal{A}(x,u,\nabla u)|^{m'}\, dx \,  dt  \\
&\leq   c_1\int_0^\infty\!\!\int_{\mathbb R^N}  u^{d+1}\psi^{k - 1}|\psi_t| \, dx  dt +
c_2 \int_0^\infty\!\!\int_{\mathbb R^N}  u^{ d+m- 1}\psi^{k - m} |\nabla \psi|^m  dx \, dt,
\end{aligned}\end{equation}
for some constants $c_1,c_2>0$ and $d>0$ satisfies \eqref{eqd}.
\end{lemma}
The proof follows line by line  that of Lemma \ref{lp1} in which we replace $d$ by $-d$.

Similar to the estimate \eqref{JJ1} we find
\begin{equation}\label{L1}\begin{aligned}
\int_{\mathbb R^N} (K\ast u^p)u^{q+d} \psi^k \, dx \, dt
\ge R^{-\beta} L(t)^2,
\end{aligned}\end{equation}
where
\begin{equation}\label{Ldef}
L(t)=\int_{\mathbb R^N}u^\tau (x,t)\psi^k(x, t)dx,
\end{equation}
where $\tau=(p+q+d)/2>1$ and $\psi$ is defined by \eqref{testfunction}. In the same way as we deduced \eqref{ineqA}, by H\"older's inequality we find

\begin{equation}\label{ineq2A}
\begin{aligned}
\int_0^\infty\!\!\int_{\mathbb R^N}  u^{d+1}\psi^{k - 1}|\psi_t| \, dx  dt
&\leq \Big(\int_0^\infty\!\!\int_{\mathbb R^N}   u^\tau \psi^k\Big)^{1/\sigma}
\Big(\int_0^\infty\!\!\int_{\mathbb R^N} \psi^{k-\sigma'}|\psi_t|^{\sigma'}\Big)^{1/\sigma'}\\[0.05in]
&= \Big(\int_0^{2R^\gamma}L(t) dt\Big)^{1/\sigma}
\Big(\int_0^\infty\!\!\int_{\mathbb R^N} \psi^{k-\sigma'}|\psi_t|^{\sigma'}\Big)^{1/\sigma'}\\[0.05in]
&\leq c R^{\gamma/(2\sigma)} \Big(\int_0^{2R^\gamma}L(t)^2 dt\Big)^{1/(2\sigma)}
\Big(\int_0^\infty\!\!\int_{\mathbb R^N} \psi^{k-\sigma'}|\psi_t|^{\sigma'}\Big)^{1/\sigma'}\\[0.05in]
&\leq C \Big(\int_0^{2R^\gamma}L(t)^2 dt\Big)^{1/(2\sigma)}R^{\frac{2N+\gamma}{2\sigma'}-\frac{\gamma}{2}},\\
\end{aligned}
\end{equation}
where $c>0$ is a constant and
$$\sigma=\frac{\tau}{d+1},\qquad \sigma'=\frac \tau{\tau-d-1}.$$
Note that since $\Gamma(d)>2$ we have $\sigma>1$.

By H\"older's inequality, similarly  to the estimate \eqref{ineqC}, we have
\begin{equation}\label{ineq2C}
\begin{aligned}
\int_0^\infty\!\!\int_{\mathbb R^N}  u^{d+m-1} & \psi^{k-m}|\nabla \psi |^m
\leq \Big(\int_0^\infty\!\!\int_{\mathbb R^N}   u^\tau \psi^k\Big)^{1/\theta}
\Big(\int_0^\infty\!\!\int_{\mathbb R^N}   \psi^{k-m\theta'}|\nabla \psi|^{m\theta'}\Big)^{1/\theta'}\\[0.05in]
&= \Big(\int_0^{2R^\gamma}L(t) dt\Big)^{1/\theta}
\Big(\int_0^\infty\!\!\int_{\mathbb R^N}   \psi^{k-m\theta'}|\nabla \psi|^{m\theta'}\Big)^{1/\theta'}\\[0.05in]
&\leq cR^{\gamma/(2\theta)} \Big(\int_0^{2R^\gamma}L(t)^2 dt\Big)^{1/(2\theta)}
\Big(\int_0^\infty\!\!\int_{\mathbb R^N}   \psi^{k-m\theta'}|\nabla \psi|^{m\theta'}\Big)^{1/\theta'}\\[0.05in]
&\leq C \Big(\int_0^{2R^\gamma}L(t)^2 dt\Big)^{1/(2\theta)}
R^{\frac{2N+\gamma}{2\theta'}-m+\frac{\gamma}{2}},
\\
\end{aligned}
\end{equation}
 where
\begin{equation}\label{eta2}
\theta=\frac{\tau}{d+m-1},\quad \theta'=\frac{\tau}{\tau-d-m+1},
\end{equation}
Note that from $\Gamma(d)>2$ we derive $\eta>1$.
Take now $\gamma=m$.

Next, we use \eqref{L1}, \eqref{ineq2A} and \eqref{ineq2C} in the estimate \eqref{ineqq2.0} of Lemma \ref{lpp1} to deduce
\begin{equation}\label{la2}
\int_0^{2R^\gamma}L(t)^2 dt \leq c\Big(\int_0^{2R^\gamma}L(t)^2 dt\Big)^{\frac{1}{2\sigma}}R^{\alpha_1}+c\Big(\int_0^{2R^\gamma}L(t)^2 dt\Big)^{\frac{1}{2\theta}}R^{\alpha_2},
\end{equation}
where, by $\eqref{eqd}$ and $\gamma=m$, we have
$$
\begin{aligned}
\alpha_1&=\frac{2N+m}{2\sigma'}-\frac{m-2\beta}{2}=\frac{(N+\beta)(p+q+d)-(2N+m)(d+1)}{p+q+d}<0,\\
\alpha_2&=\frac{2N+m}{2\theta'}-\frac{m-2\beta}{2}=\frac{(N+\beta)(p+q+d)-(2N+m)(d+m-1)}{p+q+d}<0.
\end{aligned}
$$

We further apply Young's inequality in the right hand-side of \eqref{la2} and obtain
$$
\int_0^{2R^\gamma}L(t)^2 dt \leq \frac{1}{4}\int_0^{2R^\gamma}L(t)^2 dt+C R^{\alpha_1 (2\sigma)'}+\frac{1}{4}\int_0^{2R^\gamma}L(t)^2 dt +C R^{\alpha_2 (2\theta)' },
$$
that is,
$$
\int_0^{2R^\gamma}L(t)^2 dt \leq 2C\Big( R^{\alpha_1 (2\sigma)'}+ R^{\alpha_2 (2\theta)' }\Big)\longrightarrow  0\quad\mbox{ as }R\to\infty.
$$
Proceeding as in the proof of Theorem \ref{thmain} this yields $L(t)\equiv0$ in $[0,\infty)$, which gives
$u\equiv 0$, contradiction.
\qed

\section*{Acknowledgments} RF is a member of the {\em Gruppo Nazionale per
l'Analisi Ma\-te\-ma\-ti\-ca, la Probabilit\`a e le loro Applicazioni}
(GNAMPA) of the {\em Istituto Nazionale di Alta Matematica} (INdAM)
and she was partly supported by  {\em Fondo Ricerca di
Base di Ateneo Esercizio} 2017-19 of the University of Perugia, named
{\em Problemi con non linearit\`a dipendenti dal gradiente}
and by INdAM-GNAMPA Project 2020 titled
{\em Equazioni alle derivate parziali : problemi e modelli}
(Prot\_U-UFMBAZ-2020).


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

The partially filled 5$f$ electronic shell exhibits itinerant-localized dual nature which remains a long-standing issue in condensed matter physics~\cite{LAReview}. The itinerant 5$f$ states embodied by the noticeable hybridization with conduction bands and the localization degree encoded in magnetic order enable the sensitivity of electronic structure to the external parameters, such as temperature, pressure and chemical doping. Since plutonium (Pu) situates on the edge between the itinerant and localized 5$f$ states of light and heavy actinides~\cite{RevModPhys.81.235}, 5$f$ states tend to involve in active chemical bonding and form affluent Pu-based compounds~\cite{Bauer2015Plutonium}, which exhibit novel quantum phenomena including unconventional superconductivity, magnetic order, nontrivial topology, and heavy-fermion behavior, to name a few~\cite{PhysRevLett.108.017001,PhysRevLett.111.176404,Chudo2013}.

PuSn$_3$ crystallizes in cubic AuCu$_{3}$ structure (space group $Pm$-3$m$) [see Fig.~\ref{fig:tstruct}(a)] with lattice constant 4.63 {\AA}~\cite{SARI1983301}. The temperature-independent paramagnetic order~\cite{Handbookferromag,PhysRevB.39.13115} is pretty anomalous because Pu-Pu distance is much larger than the Hill limit 3.40 {\AA}~\cite{Hilllimit} indicative of localized 5$f$ electrons. In addition, its low-temperature electrical resistivity displays a power-law relation~\cite{Brodsky1978}, which suggests the predominant distribution of conduction $s$ and $p$ states instead of 5$f$ states at the Fermi level. Hence the observed pseudogap at the Fermi level is attributed to the spin-orbit coupling, which enables the fully occupied 5$f_{5/2}$ states and partially distributed 5$f_{7/2}$ states. It is proposed that the $c$-$f$ hybridization becomes the dominant mechanism for Pu-5$f$ electron delocalization~\cite{PhysRevB.39.13115}. So far, photoemission spectroscopy, angle-resolved photoemission spectroscopy, x-ray adsorption spectroscopy and de Haas-van Alphen quantum oscillation experiments are still lacking which altogether provide subtle electronic structure of 5$f$ states.

\begin{figure}[ht]
\centering
\includegraphics[width=\columnwidth]{T_struct.pdf}
\caption{(Color online). (a) Crystal structure of PuSn$_3$. (b) Schematic picture of the first Brillouin zone of PuSn$_3$. Some high-symmetry $k$ points $X$ [0.5, 0.0, 0.0], $\Gamma$ [0.0, 0.0, 0.0], $M$ [0.5, 0.5, 0.0], and $R$ [0.5, 0.5, 0.5] are marked.
\label{fig:tstruct}}
\end{figure}

On the theoretical side, investigations concerning the structural stability, bulk modulus, band structure, electronic density of states and Fermi surface~\cite{PhysRevB.39.13115,BAIZAEE2005247,BAIZAEE2007287,PhysRevB.88.125106} attempt to elaborate narrow 5$f$ bands, Fermi surface topology and pseudogap in density of states at the Fermi level. Even though the Fermi surface and quantum oscillation are exhaustively interpreted, the missing narrow flat 5$f$ bands within the traditional density functional theory is ascribed to the underestimation of strong correlation among 5$f$ electrons~\cite{PhysRevB.88.125106}. Moreover, the unphysical density of states is mostly attributed to the exclusion of the spin-orbit coupling~\cite{BAIZAEE2007287}. Particularly, the underlying mechanism for the paramagnetic ground state has not been clearly understood. More importantly, the strongly correlated electronic states has not been fully quantified. Above all, the temperature-dependent electronic structure is rarely touched in previous literature. Consequently, it is crucial to study the temperature dependence of 5$f$ electronic structure and electronic correlation to disclose the relationship between the electronic structure and paramagnetic ground state of PuSn$_3$.

In the paper, we present the electronic structure of PuSn$_3$ dependence on temperature using the density functional theory in combination with the single-site dynamical mean-field theory. We endeavor to elucidate the itinerant-localized 5$f$ states and to uncover the strongly correlated 5$f$ states. We calculate the momentum-resolved spectral functions, density of states, Fermi surface, self-energy functions, valence state fluctuations and hybridization dynamics of PuSn$_3$. It is found that 5$f$ states become itinerant at low temperature accompanied by remarkable valence state fluctuations. Additionally, the onset of atomic multiplets and prominent $c$-$f$ hybridization at low temperature imply a temperature-driven 5$f$ itinerant-localized crossover. Especially, the change in Fermi surface topology implicates a possible Lifshitz transition induced by temperature. Finally, the orbital selective electronic correlation and hybridization dynamics are addressed.

The rest of this paper is organized as follows. In Sec.~\ref{sec:method}, the computational details are introduced. In Sec.~\ref{sec:results}, the electronic band structures, Fermi surface topology, total and partial 5$f$ density of states, 5$f$ self-energy functions, and probabilities of atomic eigenstates are presented. In Sec.~\ref{sec:dis}, we attempt to clarify some important topics about the 5$f$ itinerant-localized crossover and the hybridization gaps. Finally, Sec.~\ref{sec:summary} serves as a brief conclusion.

\section{Methods\label{sec:method}}
The well-established DFT + DMFT method combines realistic band structure calculation by DFT with non-perturbative many-body treatment of local interaction effects in DMFT~\cite{RevModPhys.68.13,RevModPhys.78.865}. The strong electronic correlation and large spin-orbit coupling are treated on the same footing.
Here we perform charge fully self-consistent calculations to explore the temperature-dependent electronic structures of PuSn$_3$ using DFT + DMFT method. The self-consistent implementation of this method is divided into DFT and DMFT parts, which are solved separately by using the \texttt{WIEN2K} code~\cite{wien2k} and the \texttt{EDMFTF} package~\cite{PhysRevB.81.195107}.

In the DFT calculation, the experimental crystal structure of PuSn$_3$~\cite{SARI1983301} was used and ignored the thermal expansion. The generalized gradient approximation was adopted to formulate the exchange-correlation functional~\cite{PhysRevLett.77.3865}. Additionally, the spin-orbit coupling was taken into account in a second-order variational manner. The muffin-tin radius for Pu and Sn were chosen as 2.7 au and 2.4 au, respectively. The $k$-points' mesh was $15 \times 15 \times 15$ and $R_{\text{MT}}K_{\text{MAX}} = 8.0$.

In the DMFT part, 5$f$ orbitals of plutonium were treated as correlated. The four-fermions' interaction matrix was parameterized using the Coulomb interaction $U = 5.0$~eV and the Hund's exchange $J_H=0.6$~eV via the Slater integrals~\cite{PhysRevB.59.9903}. The fully localized limit scheme was used to calculate the double-counting term for impurity self-energy function~\cite{jpcm:1997}. The vertex-corrected one-crossing approximation (OCA) impurity solver~\cite{PhysRevB.64.115111} was employed to solve the resulting multi-orbital Anderson impurity models.
 Note that we not only utilized the good quantum numbers $N$ (total occupancy) and $J$ (total angular momentum) to classify the atomic eigenstates, but also made a severe truncation ($N \in$ [3, 7]) for the local Hilbert space~\cite{PhysRevB.75.155113} to reduce the computational burden. The convergence criteria for charge and energy were $10^{-5}$ e and $10^{-5}$ Ry, respectively. It is worth noting that the direct output of OCA impurity solver is real axis self-energy $\Sigma (\omega)$ which was applied to calculate the momentum-resolved spectral functions $A(\mathbf{k},\omega)$ and density of states $A(\omega)$, as well as other physical observables.

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

\begin{figure*}[th]
\centering
\includegraphics[width=\textwidth]{T_akw.pdf}
\caption{(Color online). Momentum-resolved spectral functions $A(\mathbf{k},\omega)$ of PuSn$_3$ as a function of temperature under ambient pressure calculated by the DFT + DMFT method. The horizontal lines denote the Fermi level. (a) 580 K. (b) 290 K. (c) 116 K. (d) 29 K.
\label{fig:akw}}
\end{figure*}

\begin{figure*}[th]
\centering
\includegraphics[width=\textwidth]{T_dos.pdf}
\caption{(Color online). Electronic density of states of PuSn$_3$ obtained by the DFT + DMFT method.
Total density of states (thick solid lines) and partial 5$f$ density of states (color-filled regions) of 580 K (a), 290 K (b), 29 K (c). These peaks resulting from the quasiparticle multiplets are denoted with ``P1'', ``P2'', ``P3'', and ``P4''. The $j$-resolved 5$f$ partial density of states with $5f_{5/2}$ and $5f_{7/2}$ components represented by purple and green lines, respectively. 580 K (d), 290 K (e), 29 K (f). (g) The evolution of 5$f$ density of states against temperature in the vicinity of Fermi level. (h) The height of the central quasiparticle peak h$_{\rm QP}$(5$f_{5/2}$) and 5$f$ spectral weight at the Fermi level $A_{5f}$($\omega$ = 0) as a function of temperature.
\label{fig:dos}}
\end{figure*}

\begin{figure*}[th]
\centering
\includegraphics[width=\textwidth]{T_fs.pdf}
\caption{(Color online).
Three-dimensional Fermi surface and two-dimensional Fermi surface of PuSn$_3$ calculated by the DFT + DMFT method at 580 K (a, b, c) and 116 K (d, e, f). Two-dimensional Fermi surface are on the $k_x-k_y$ plane (with $k_z = \pi/2$). They are visualized with different colors.
\label{fig:FS}}
\end{figure*}

\begin{figure*}[th]
\centering
\includegraphics[width=\textwidth]{T_prob.pdf}
\caption{(Color online). Probabilities of $5f^{4}$ (blue), $5f^{5}$ (red), $5f^{6}$ (green) configurations (a) and $5f^{3}$ (purple), $5f^{7}$ (cyan) (b) for PuSn$_3$ derived by DFT + DMFT calculations. (c) 5$f$ occupancy as a function of temperature.
\label{fig:prob}}
\end{figure*}

\subsection{Momentum-resolved spectral functions}
It is illuminating to inspect the momentum-resolved spectral functions which encode interesting features of PuSn$_3$ [see Fig.~\ref{fig:akw}]. We examine the quasiparticle bands as a function of temperature to uncover the inherent nature of 5$f$ states. First of all, it is necessary to evaluate the reliability of our calculated band structures.
In comparison to the band structure of paramagnetic ground state in previous DFT calculations~\cite{PhysRevB.88.125106}, the low temperature momentum-resolved spectral functions [see Fig.~\ref{fig:akw}(d)] share similar characteristics.
For example, a striking electron pocket locates at $\Gamma$ point about --3 eV with the hole pockets at $X$ and $M$ points below the Fermi level. Generally, the band dispersion around $X$ and $M$ points are reasonably in line with the previous DFT calculations~\cite{PhysRevB.88.125106}. In special, the conduction bands along $\Gamma$ - $M$ high-symmetry line are very much alike.
Even though the flat bands parallel to the vicinity of Fermi level at low temperature [see Fig.~\ref{fig:akw}(d)] are absent in the DFT calculations, the basic features of energy bands are consistent with each other, evincing the validity of our results. The discrepancy in the existence of narrow 5$f$ bands around the Fermi level indicates the underestimation of strong correlation among 5$f$ electrons within traditional density functional theory, which demands a more rigorous and sophisticated quantum many-body approach.

Next we turn to the temperature effects on quasiparticle bands of PuSn$_3$. At high temperature [see Fig.~\ref{fig:akw}(a)], only conduction bands with noticeable dispersions cross the Fermi level, accompanied by the emergence of lower and upper Hubbard bands in the energy range of --3 eV $\sim$ --1 eV and 1 eV $\sim$ 3 eV, respectively.
When temperature decreases, 5$f$ electrons gradually get coherent and certain quasiparticle bands start to build up near the Fermi level. These bands are not evident at first because of their small quasiparticle weight and tiny band intensity.
Combined with the density of states in Fig.~\ref{fig:dos}, these flat narrow bands mainly root from 5$f$ states, signifying the itinerant tendency of 5$f$ states. As temperature further lowers, quasiparticle weight grows obvious and apparent narrow 5$f$ bands develop at --0.9 eV, --0.47 eV and the Fermi level, which hybridize with 5$p$ conduction bands of Sn atom to open evident hybridization gaps. It is noted that nearly dispersionless quasiparticle bands are split by spin-orbit coupling with energy gap about 0.43 eV. Meanwhile, the profiles of conduction bands above and below the Fermi level remain almostly unchanged with respect to temperature.
Accordingly, the significant quasiparticle bands, remarkable $c-f$ hybridization and salient spectral weight of 5$f$ states jointly imply intensifying itinerancy of 5$f$ bands with decreasing temperature. Thus, it is proposed that a potential localized to itinerant crossover of 5$f$ correlated electronic states is induced by temperature in PuSn$_3$.

\subsection{Density of states}
To elucidate the evolution of itinerant degree of 5$f$ correlated electronic states upon temperature, we explore the density of states and quasiparticle multiplets in a wide temperature range of 29 K $\sim$ 580 K. It is noticed that basic features exist in the whole temperature range [see Fig.~\ref{fig:dos}(a)-(c)]. Firstly, the overall profile of total density of states resemble each other including broad ``humps" from --3 eV $\sim$ --1 eV and 1 eV $\sim$ 3 eV, which are mainly assigned to the lower and upper Hubbard bands of plutonium's 5$f$ orbitals. Secondly, the fourteen-fold degenerated 5$f$ states are split into six-fold degenerated 5$f_{5/2}$ and eight-fold degenerated 5$f_{7/2}$ subbands [see Fig.~\ref{fig:dos}(e)-(f)] because of spin-orbit coupling~\cite{RevModPhys.81.235,PhysRevB.102.245111,PhysRevB.101.125123}. A sharp and narrow quasiparticle (QP) peak develops in the Fermi level, mostly belonging to 5$f_{5/2}$ orbital. Meanwhile, the two satellite peaks ``P1'' and ``P2'' at --0.9 eV and --0.47 eV with energy gap about 0.43 eV are ascribed to 5$f_{5/2}$ and 5$f_{7/2}$ orbitals, respectively.
Above the Fermi level, the reflected peaks ``P3'' and ``P4'' locate at 0.47 eV and 0.9 eV with respect to the central quasiparticle. So the five representative peaks are called ``quasiparticle multiplets''~\cite{2102.02034,2102.03085}. To trace the origin of these quasiparticle multiplets, it is deduced that ``P3'' and ``P4'' peaks are formed by a mix of 5$f_{5/2}$ and 5$f_{7/2}$ orbitals. Overall, it is mentioned that the atomic multiplets are induced by 5$f$ valence fluctuations, which leave fingerprints on the 5$f$ photoemission spectroscopy of PuSn$_3$.

As to the temperature-dependent coherence and itinerancy of 5$f$ electrons, it is instructive to examine the behavior of quasiparticle multiplets around the Fermi level. It reveals that the quasiparticle peak is too small to be seen in the Fermi level at high temperature, implicating the mostly localized and mainly incoherent 5$f$ states. When temperature slowly descends, the coherence of 5$f$ valence electron gradually builds up. Since the spectral weights of upper and lower Hubbard bands transfer to the Fermi level, the quasiparticle peak around the Fermi level rises up quickly. It is apparent that the height of 5$f_{5/2}$ quasiparticle peak [see Fig.~\ref{fig:dos}(g)] ascends swiftly and acutely, resulting in a sharp and intense quasiparticle peak at the Fermi level. Meanwhile, the quasipartilce weight at the Fermi level progressively magnifies, giving rise to coherent 5$f$ states. In a word, the significant hybridization between 5$f$ electrons of Pu atom and 5$p$ bands of Sn atom suggests the increasing itinerancy of 5$f$ states. The intensifying $c-f$ hybridization contributes to unveil the paramagnetic ground state. Hence the valence state fluctuations become predominant concurrently, which manifests the mixed-valence nature and hints the potential temperature-induced localized to itinerant crossover of 5$f$ states.

\subsection{Fermi surface topology}
Fermi surface topology is an effective physical quantity to capture the temperature-dependent electronic structure of PuSn$_3$. Figure~\ref{fig:FS} depicts the Fermi surface topology at two typical temperatures 580 K and 116 K, respectively. There are two doubly degenerated bands crossing the Fermi level (No. of bands: 18 and 19, 20 and 21), which are labeled by $\alpha$ and $\beta$, respectively. Note that the Fermi surface topology of $\alpha$ band takes an ellipsoid shape, which is in coincidence with those in previous DFT calculations~\cite{PhysRevB.88.125106}. Nevertheless, the Fermi surface topology of $\beta$ band resembles an anisotropic form, which acts different from the previous results~\cite{PhysRevB.88.125106}. The discrepancy of $\beta$ band might arise from the temperature effect on Fermi surface because the DFT calculations~\cite{PhysRevB.88.125106} are carried out at zero temperature and we perform the calculations at finite temperature.

As expected, the Fermi surface is sensitive to varying temperature. It is detected that the volume of Fermi surface for $\alpha$ band expands evidently [see Fig.~\ref{fig:FS} (a) and (d)] when temperature diminishes, providing sufficient evidence for the itinerant 5$f$ electrons at low temperature. Since the three-dimensional Fermi surface of $\beta$ band is hard to discern its inner structure, the two-dimensional Fermi surface of $\alpha$ and $\beta$ bands are visualized in Fig.~\ref{fig:FS} (c) and (f). Obviously, both $\alpha$ and $\beta$ bands intersect the $\Gamma$ - $X$ line and the distance between intersections of two bands almost remains unchanged with decreasing temperature, which is in conformity with the momentum-resolved spectral functions [see Fig.~\ref{fig:akw}]. On the contrary, it appears that the distance between the $\beta$ band intersections with Fermi level along the $M$ - $X$ line shrinks as temperature lowers. In particular, the anisotropic Fermi surface topology of $\beta$ band undergoes variation in surface envelope along the $\Gamma$ - $M$ line on the two-dimenional Fermi surface. In consequence, the Fermi surface topology actually changes dramatically with decreasing temperature, which demonstrates a probable Lifshitz transition for 5$f$ states and the advent of 5$f$ localized to itinerant crossover driven by temperature. Since the Fermi surface could be detected by subsequent dHvA quantum oscillation, the experimental results could clarify the underlying mechanism behind the paramagnetic ground state of PuSn$_3$ provided that no signature of Fermi surface nesting is observed.

\begin{figure*}[th]
\centering
\includegraphics[width=\textwidth]{T_sig.pdf}
\caption{(Color online). Real-frequency self-energy functions of PuSn$_3$ obtained by the DFT + DMFT method. (a)-(b) $Z|{\rm Im}\Sigma (\omega)|$ for the 5$f_{5/2}$ and 5$f_{7/2}$ states. $Z$ is the renormalization factor. (c) Electron effective masses for the 5$f_{5/2}$ and 5$f_{7/2}$ states as a function of temperature, where the green arrow denotes maximum electron effective mass for the 5$f_{5/2}$ states.
\label{fig:sig}}
\end{figure*}

\subsection{Valence state fluctuations}
It is worth mentioning that $\delta$-Pu displays obvious mixed-valence behavior with noninteger occupation number deviating from nominal value 5.0. The 5$f$ electron atomic eigenstates derived from the output of DMFT many-body states shed light on the valence state fluctuations and related mixed-valence behavior.
 Here $p_\Gamma$ is used to quantify the probability of 5$f$ electrons which stay in each atomic eigenstate $\Gamma$. Then the average 5$f$ valence electron is defined as $\langle n_{5f} \rangle = \sum_\Gamma p_\Gamma n_\Gamma$, where $n_\Gamma$ denotes the number of electrons in each atomic eigenstate $\Gamma$. Finally, the probability of 5$f^n$ electronic configuration can be expressed as $\langle w(5f^{n}) \rangle = \sum_\Gamma p_\Gamma \delta (n-n_\Gamma)$.

The calculated probabilities of 5$f^n$ electronic configuration for PuSn$_3$ are visualized in Fig.~\ref{fig:prob}(a) and (b). Apparently, the probability of 5$f^5$ electronic configuration is dominating, while the contributions of 5$f^3$ and 5$f^7$ electronic configurations are too small to be neglected. As is shown in Table~\ref{tab:prob}, at 580 K, the probability of 5$f^5$ electronic configuration accounts for 83.1\%, followed by the probabilities of 5$f^4$ and 5$f^6$ electronic configurations about 9.3\% and 7.4\%, respectively. Meanwhile, the occupation of 5$f$ electrons is approximately 4.98. When temperature gradually diminishes, the probability of 5$f^5$ electronic configuration declines slowly. At the same time, the probability of 5$f^4$ electronic configuration slightly decreases and the probability of 5$f^6$ electronic configuration appreciably augments. At 29 K, it is noticed that the probability of 5$f^5$ electronic configuration drops to 76.0\%, with 5$f^4$ and 5$f^6$ electronic configurations about 7.4\% and 16.3\%, respectively. At this stage, the occupation of 5$f$ electrons experiences a minor increase to 5.09.
In this scenario, 5$f$ valence electrons are prone to spend more time in the 5$f^6$ configuration, rendering the valence fluctuation and promoting quasiparticle multiplets. In low-temperature regime, it is speculated that PuSn$_3$ is a mixed-valence metal.
So it is inspiring to refer to isostructual compound PuIn$_3$~\cite{2102.03085}, whose electronic configuration probably suggests insensitivity to varying temperature. Analogous to cubic phase Pu$_3$Ga~\cite{2102.02034}, strong valence fluctuations play a key role in generating the quasiparticle multiplets at the Fermi level, and regulating the effective 5$f$ valence electrons.

Interestingly, the temperature-dependent patterns of 5$f^4$ and 5$f^6$ electronic configurations intersect about 387 K demonstrating the initiating itinerancy of 5$f$ states and appearance of electron coherence. Coincidently, the electron effective mass for 5$f_{5/2}$ states obtained from self-energy functions according to Eq
.~\ref{eqsigma} reaches the maximum value at 387 K [see Fig.~\ref{fig:sig}(c)]. Generally speaking, electronic correlation is usually strong in localized 5$f$ electrons and grows weaker in itinerant 5$f$ electrons. The coincidence between 5$f$ localized-itinerant nature and strong electronic correlation unravels the hidden connection of itinerant degree and electronic correlation, which requires further exploration.

\begin{table}[th]
\caption{Probabilities of $5f^{3}$, $5f^{4}$, $5f^{5}$, $5f^{6}$, and $5f^{7}$ for PuSn$_3$ at temperatures 580 K and 29 K, respectively.
 \label{tab:prob}}
\begin{ruledtabular}
\begin{tabular}{cccccc}
Temperatures & $5f^{3}$ & $5f^{4}$ & $5f^{5}$ & $5f^{6}$ & $5f^{7}$ \\
\hline
580 K    & 1.357$\times 10^{-3}$ & 0.093  & 0.831 & 0.074 & 4.934$\times 10^{-4}$ \\
29 K   & 9.921$\times 10^{-4}$ & 0.074  & 0.760 & 0.163 & 1.511$\times 10^{-3}$ \\
\end{tabular}
\end{ruledtabular}
\end{table}

\subsection{Self-energy functions}
In general, the electronic correlations are encapsulated in the electron self-energy functions~\cite{RevModPhys.68.13,RevModPhys.78.865}. Figure~\ref{fig:sig} illustrates the renormalized imaginary part of self-energy functions $Z|{\rm Im}\Sigma(\omega)|$ for 5$f_{5/2}$ and 5$f_{7/2}$ states, as well as the electron effective mass.
Here $Z$ means the quasiparticle weight or renormalization factor, which denotes the electronic correlation strength and can be obtained from the real part of self-energy functions via the following equation~\cite{RevModPhys.68.13}:
\begin{equation}
Z^{-1} = \frac{m^\star}{m_e} = 1 - \frac{\partial \text{Re} \Sigma(\omega)}{\partial \omega} \Big|_{\omega = 0}. \label{eqsigma}
\end{equation}

As is stated above, $Z|{\rm Im}\Sigma(0)|$ embodies the low-energy electron scattering rate~\cite{PhysRevB.99.125113}.
Both $Z|{\rm Im}\Sigma_{5f_{5/2}}(0)|$ and $Z|{\rm Im}\Sigma_{5f_{7/2}}(0)|$ approximate zero at low temperature, suggesting the itinerant 5$f$ states and the metallic feature. And they develop finite values at high temperatues indicating increscent low-energy electron scattering rate.
With the increment of temperature, $Z|{\rm Im}\Sigma_{5f_{5/2}}(\omega)|$ grows up rapidly and saturates in the energy range of [--0.5 eV, --0.25 eV].
Concurrently, $Z|{\rm Im}\Sigma_{5f_{7/2}}(\omega)|$ climbs up sharply in the high-energy regime ($|\omega|$ $>$ 0.5 eV). Therefore the enhancement of $Z|{\rm Im}\Sigma_{5f_{7/2}}(\omega)|$ becomes more remarkable than that of $Z|{\rm Im}\Sigma_{5f_{5/2}}(\omega)|$, which concretises the growing localization of 5$f_{7/2}$ state and explains its incoherent nature at high temperature. Now that self-energy functions of 5$f_{5/2}$ and 5$f_{7/2}$ states evince differentiated temperature-dependent patterns, it is proposed that 5$f$ electronic correlation are orbital dependent, which commonly exists in actinide compounds.

The evaluated electron effective masses $m^*$ for 5$f_{5/2}$ and 5$f_{7/2}$ states~\cite{RevModPhys.68.13} according to equation~\ref{eqsigma} are given in Fig.~\ref{fig:sig}(c). With the decrement of temperature, the electron effective mass for 5$f_{5/2}$ states surges up rapidly to reach the maximum value up to 40 $m_e$ at 387 K, followed by a monotonically decreasing trend below 387 K. The diminishing electron effective mass for 5$f_{5/2}$ states implies the weakening electronic correlations at low temperature, which is associated with the enhancive itinerancy of 5$f$ states.
Conversely, the electron effective mass for 5$f_{7/2}$ states increases steadily until 387 K. After that it dramatically saturates to a finite value at low temperature. The maximum electron effective mass for 5$f_{7/2}$ states at low temperature hints strong electronic correlation and localized 5$f$ electrons.
The distinct behaviors of electron effective masses for 5$f_{5/2}$ and 5$f_{7/2}$ states are associated with the intrinsic electron correlated characteristic, which confirms the orbital selective correlated states.

\subsection{Discussions\label{sec:dis}}
In this section, we discuss about the heavy-fermion state and temperature-dependent itinerant 5$f$ states to unveil the strongly correlated electronic structure of PuSn$_3$ and underlying mechanism of the paramagnetic ground state.

\textbf{Electronic heat capacity.}
To explore the heavy fermion behavior of PuSn$_3$, we evaluate the specific heat coefficient within the framework of Fermi-liquid theory. Taking into account the electronic degree of freedom and lattice vibrations, the heat capacity of solid contains two parts $C_v(T)=\gamma T+\beta T^3$. The linear part of heat capacity $\gamma T$ gives the linear specific heat coefficient $\gamma$, which is expressed by the following equation:
\begin{equation}
\gamma = \pi k^2_{B} \sum_{\alpha} \frac{A_{\alpha}(0)}{Z_{\alpha}},
\end{equation}
where $\alpha$ is the orbital index, $A_{\alpha}(0)$ is the spectral weight at the Fermi level, and $Z_{\alpha}$ is the orbital-resolved renormalization factor~\cite{PhysRevLett.101.056403,RevModPhys.68.13}. Since the effective mass enhancement is usually not quite substantial in Pu-based compounds, the specific heat coefficient $\gamma >$ 100 mJ/(mol$\times$K$^2$) is defined as a threshold for Pu-based heavy-fermion compound~\cite{Bauer2015Plutonium}.
Under the neglection of contributions from the conduction bands, the calculated electronic specific heat coefficient of 5$f$ states at 387 K is approximately 96 mJ/(mol$\times$K$^2$). Referring to typical Pu-based heavy-fermion compound PuIn$_3$ with specific heat coefficient 307 mJ/(mol$\times$K$^2$), it is conjectured that PuSn$_3$ might be a promising candidate of Pu-based heavy-fermion compound.

\textbf{Orbital selectivity.}
Here we analyze the temperature-driven itinerant to localized crossover and orbital dependent electronic correlations of 5$f$ states for PuSn$_3$. Primarily, flat narrow quasiparticle bands of 5$f$ states emerge at the Fermi level with amplifying quasiparticle weight and they hybridize with conduction bands to open obvious gaps. Then the itinerancy and coherence of 5$f$ electrons get strengthening at low temperature.
Due to spin-orbit coupling, 5$f$ states are split into 5$f_{5/2}$ and 5$f_{7/2}$ states, which become itinerant asynchronously.
Combined with the density of states, it is deduced that the quasiparticle peak rises up quickly at low temperature accompanied by predominant quasiparticle multiplets around the Fermi level. In addition, it is pointed out that quasiparticle multiplets are evoked by valence state fluctuations. It is observed that the central quasiparticle peak is mainly contributed by 5$f_{5/2}$ states which grows coherence at higher temperature than that of 5$f_{7/2}$ states.
Ultimately, the evolution pattern of 5$f^n$ electronic configuration demonstrates remarkable valence state fluctuations and distinct mixed-valence behavior, which is closely related to the enhancing itinerancy and developing coherence of 5$f$ states. Briefly, it is speculated that temperature-induced localized to itinerant crossover of 5$f$ states is orbital differentiation.
On the other hand, self-energy functions and electron effective masses of 5$f_{5/2}$ and 5$f_{7/2}$ states display differentiated temperature-dependent behaviors, which signifies the orbital dependent 5$f$ electronic correlation in reminiscence with Pu$_3$Ga~\cite{2102.02034}.

\begin{figure*}[t!]
\centering
\includegraphics[width=\textwidth]{T_hyb.pdf}
\caption{(Color online). Temperature dependent quasiparticle band structures (along the $\Gamma$ - $M$ direction) of cubic phase PuSn$_3$ acquired by the DFT + DMFT method. (a) $T$ = 580 K. (b) $T$ = 116 K. (c) $T$ = 29 K.
The white horizontal dashed lines denote the Fermi level. Here, we utilized a periodical Anderson model [see Eq.~\ref{eq:pam}] to fit the low-energy band structures. The colorful dashed lines are the fitting results. The renormalized $5f$ energy level $\epsilon_0$ is denoted by purple line, the unrenormalized band dispersion for conduction electrons $\epsilon(k)$ is plotted in black line. The upper and lower branches of hybridized bands are denoted with green and blue lines, respectively.
\label{fig:hyb}}
\end{figure*}

\textbf{Evolution of hybridization gap.}
As is already illustrated in Fig.~\ref{fig:akw}(a)-(d), strong hybridization exists between the 5$f$ bands and conduction bands at low temperature. The $c-f$ hybridization will open a hybridization gap for the conduction bands. Since this gap is adjacent to the Fermi level, the physical properties of PuSn$_{3}$ should be affected largely by it. Hence, it is considerable to determine the size of the hybridization gap and elucidate its temperature dependence relation. Phenomenologically, the low-energy hybridized bands can be well described by a simple mean-field hybridization band picture (i.e. the periodical Anderson model)~\cite{RevModPhys.92.011002}. According to this picture, the energy dispersions are expressed as:
\begin{equation}
\label{eq:pam}
E_{\pm}(k) = \frac{[\epsilon_0 + \epsilon(k)] \pm \sqrt{[\epsilon_0 - \epsilon(k)]^2 + 4|V_k|^2}}{2},
\end{equation}
where $\epsilon_0$ means the renormalized $5f$ energy level, $\epsilon(k)$ is the unrenormalized band dispersion for conduction electrons, and $V_k$ denotes the strength of hybridization. In the left side of this equation, the ``+'' and ``--'' symbols mean the upper and lower branches of hybridized bands, respectively. In Fig.~\ref{fig:hyb}(a) the band dispersion data at $T = 580$~K are shown. At such a high temperature the hybridization is negligible (i.e. $|V_k| = 0$ and $\epsilon_0 = 0$), so Eq.~(\ref{eq:pam}) is simplified to $E_{\pm}(k) = \epsilon(k)$. Thus, we used the data at $T = 580$~K to calibrate $\epsilon(k)$ (see the black dashed lines in Fig.~\ref{fig:hyb}). The band structures at $T = 116$~K and $29$~K are shown in Fig.~\ref{fig:hyb}(b) and (c), respectively. After fitting Eq.~(\ref{eq:pam}), we derive $\epsilon_0 = -26.5$~meV and $|V_k| = 42.5$~meV for $T = 116$~K, and $\epsilon_0 = -19$~meV and $|V_k| = 50$~meV for $T = 29$~K. Thus, the direct hybridization gaps ($\Delta \approx 2|V_k|$) are 85~meV and 100~meV at 116~K and 29~K, respectively. The results demonstrate that the 5$f$ energy level are pulled toward the Fermi level at low temperature, enhancing the $c-f$ hybridization strength, so as to enlarge the hybridization gap $\Delta$. Similarly, the quasiparticle multiplets also hybridize with the conduction bands. Therefore, multiple hybridization gaps will open at low temperature, which is quite complex and further work is ongoing.

\section{conclusion\label{sec:summary}}
In summary, we studied the detail electronic structure of PuSn$_3$ by employing a state-of-the-art first-principles many-body approach. The temperature dependence of itinerant to localized crossover of 5$f$ states and the correlated electronic states were addressed systematically. As temperature declines, 5$f$ electrons develop coherence and becomes itinerant, exhibiting outstanding quasiparticle multiplets and pronounced valence state fluctuations, accompanied by conspicuous spectal weight and noteworthy $c-f$ hybridization. Especially, the change in Fermi surface topology induced by temperature hints a Lifshitz transition which could be detected by quantum oscillation experiment. Accordingly, 5$f$ states experience a temperature-driven itinerant to localized crossover. Above all, 5$f$ states manifest orbital selective electronic correlation, expressing themselves as orbital dependent electron effective mass and renormalized bands. Our calculated results not only provide a comprehensive picture about how the 5$f$ correlated electronic states evolve with respect to temperature, but also gain deep insight into the complex electronic structure and mysterious paramagnetic state of PuSn$_3$. Further studies about the other Pu-based compounds are undertaken.

\begin{acknowledgments}
This work was supported by the National Natural Science Foundation of China (No.~11874329, No.~11934020, No.~22025602), and the Science Challenge Project of China (No.~TZ2016004).
\end{acknowledgments}