\section{Introduction}
This paper is devoted to studying the global dynamics of solutions to the energy subcritical defocusing semilinear wave equation
\begin{equation}
  \label{eq:NLW:semi:3d}
  \Box\phi=|\phi|^{p-1}\phi,\quad \phi(0, x)=\phi_0(x),\quad \pa_t\phi(0, x)=\phi_1(x)
\end{equation}
in $\mathbb{R}^{3+1}$ with $1<p<5$.

The existence of global classical solution has early been obtained by J\"{o}rgens in \cite{Jorgens61:energysub:NLW:lowerd}, with various kinds of extensions in \cite{Brenner79:globalregularity:NW}, \cite{Brenner81:globalregularity:d9}, \cite{Pecher76:NLW:global}, \cite{segal63:semigroup}, \cite{Vonwahl75:NW}, \cite{velo85:global:sol:NLW}, \cite{Velo89:globalsolution:NLW} and references therein.
 Strauss in \cite{Strauss:NLW:decay} investigated the asymptotics of the global solution for the superconformal case $3\leq p<5$. He showed that the solution scatters to linear wave in $H^1$ with compactly supported initial data (for asymptotic behaviours of this type, we refer to the author's companion paper \cite{yang:scattering:NLW} and references therein). This scattering result relied on the pointwise decay estimate $t^{\ep-1}$ for any $\epsilon >0$, which has later been improved to $t^{-1}$ for the strictly superconformal case $3<p<5$ and to $t^{-1}\ln t$ for the conformal case $p=3$ by Wahl in \cite{vonWahl72:decay:NLW:super}. The starting point of these asymptotic decay estimates is the conservation of approximate conformal energy derived by using the conformal Killing vector field as multiplier. The superconformal structure of the equation plays the role that the corresponding conformal energy (including the potential energy contributed by the nonlinearity) is controlled by the initial data as the spacetime integral arising from the nonlinearity is nonnegative, which in particular implies the uniform bound of the $L^2$ norm of the solution and the time decay of the potential energy. These a priori bounds force the solution to decay in time by viewing the nonlinearity as inhomogeneous term of the linear wave equation.
Another geometric point of view to see this conformal structure is the method of conformal compactification (see \cite{ChristodoulouYangM}, \cite{ChDNull}), based on which together with the representation formula for linear wave equation,
Bieli-Szpak obtained shaper decay estimate for the solution with compactly supported initial data in \cite{Roger:3DNLW:symmetry}, \cite{Bieli:3DNLW}.

To go beyond the superconformal case, Pecher in \cite{Pecher82:decay:3d} observed that the potential energy still decays in time but with a weaker decay rate. For this case, the spacetime integral arising from the nonlinearity mentioned above changes sign to be negative. This term could be controlled by using Gronwall's inequality with the price that the conformal energy grows in time with a rate linearly depending on the coefficient. Since the conformal energy contains weights in $t$, the potential energy still decays in time when $p$ is not too small (sufficiently close to $3$ so that the coefficient is small). This weaker energy decay estimate is sufficiently strong to conclude the pointwise decay estimate for the solution when $p>\frac{1+\sqrt{13}}{2}$. As a consequence, the solution scatters in energy space for $p>2.7005$.

 The aim of this paper are two folds: firstly we obtain pointwise decay estimate for the solution with data in some weighted energy space which is weaker than the conformal energy space required in previous works. We prove that the solution decays as quickly as linear waves (with the same initial data) for all $p>\frac{1+\sqrt{17}}{2}$, covering additional part of the subconformal range. Secondly for even smaller $p$ with lower bound $2$, we show that the solution decays at least $t^{-\frac{1}{3}}$. This decay estimate immediately leads to the scattering result in energy space for $p>2.3542$, hence refining Pecher's pointwise decay estimates and scattering result in $\mathbb{R}^{3+1}$.

More precisely, for some fixed constant $1<\ga_0<2$ define the weighted energy norm of the initial data
\begin{align*}
  \mathcal{E}_{k,\ga_0}[\phi]=\sum\limits_{l\leq k}\int_{\mathbb{R}^3}(1+|x|)^{\ga_0+2l}(|\nabla^{l+1}\phi_0|^2+|\nabla^l \phi_1|^2)+(1+|x|)^{\ga_0}|\phi_0|^{p+1}dx.
\end{align*}
Then we have
\begin{Thm}
\label{thm:main}
 Consider the Cauchy problem to the energy subcritical defocusing semilinear wave equation \eqref{eq:NLW:semi:3d}. For initial data $(\phi_0, \phi_1)$ bounded in $\mathcal{E}_{1, \ga_0}[\phi]$ for some constant $1<\ga_0<2$, the solution is global in time and  satisfies the following decay estimates:
\begin{itemize}
\item For the case when
\[
\frac{1+\sqrt{17}}{2}<p<5, \quad \max\{\frac{4}{p-1}-1, 1\}<\ga_0<\min\{p-1, 2\},
\]
then
\begin{equation*}
|\phi(t, x)|\leq C (1+\mathcal{E}_{1, \ga_0}[\phi] )^{\frac{p-1}{2}}(1+t+|x|)^{-1}(1+||x|-t|)^{-\frac{\ga_0-1}{2}};
\end{equation*}
\item Otherwise if $2<p\leq \frac{1+\sqrt{17}}{2}$ and $1<\ga_0<p-1$, then
\begin{equation*}
|\phi(t, x)|\leq C \sqrt{\mathcal{E}_{1, \ga_0}[\phi] }  (1+t+|x|)^{-\frac{3+(p-2)^2}{(p+1)(5-p)}\ga_0}(1+||x|-t|)^{-\frac{\ga_0}{p+1}}
\end{equation*}
for some constant $C$ depending on $\ga_0$, $p$ and the zeroth order weighted energy $\mathcal{E}_{0, \ga_0}[\phi] $.
\end{itemize}
 \end{Thm}
 As a consequence of the above pointwise decay estimate, we extend Pecher's scattering result to a larger range of $p$. Recall the linear operator $\mathbf{L}(t)$ defined in \cite{yang:scattering:NLW}
 \begin{align*}
   \Box\mathbf{L}(t)(f, g)=0,\quad \mathbf{L}(0)(f, g)=f(x),\quad \pa_t \mathbf{L}(0)(f, g)=g(x).
 \end{align*}
\begin{Cor}
 \label{cor:scattering:3D}
 For $p>p_*$ (defined in the last section and $p_*<2.3542$) and initial data bounded in $\mathcal{E}_{1, p-1}[\phi]$, the solution $\phi$ of \eqref{eq:NLW:semi:3d} is uniformly bounded in the following mixed spacetime norm
 \begin{align*}
   \|\phi\|_{L_t^p L_x^{2p}}<\infty.
 \end{align*}
 Consequently the solution scatters in energy space, that is, there exists pairs $(\phi_0^{\pm}(x), \phi_1^{\pm}(x))$ such that
 \begin{align*}
      \lim\limits_{t\rightarrow\pm\infty}\|\phi(t, x)-\mathbf{L}(t)(\phi_0^{\pm}(x), \phi_1^{\pm}(x))\|_{\dot{H}_x^1}+\| \pa_t\phi(t,x)-\pa_t \mathbf{L}(t)(\phi_0^{\pm}(x), \phi_1^{\pm}(x))\|_{ L_x^2}=0.
\end{align*}
\end{Cor}
We give several remarks.
\begin{Remark}
 One can also derive the pointwise decay estimates for the derivatives of the solution by assuming the boundedness of the second order weighted energy of the initial data.
\end{Remark}

\begin{Remark}
 The precise decay estimate obtained by Pecher in \cite{Pecher82:decay:3d} is the following
\[
|\phi|\leq C t^{\frac{6+2p-2p^2}{3+p}+\ep},\quad \frac{1+\sqrt{13}}{2}<p\leq 3
\]
with initial data bounded in $\mathcal{E}_{1, 2}[\phi]$. Theorem \ref{thm:main} improves this decay estimate with weaker assumption on the initial data.
\end{Remark}

\begin{Remark}
  Note that the solution to the linear wave equation
  \[
  \Box\phi^{lin}=0,\quad \phi^{lin}(0, x)=\phi_0(x),\quad \pa_t\phi^{lin}(0, x)=\phi_1
  \]
  with data $(\phi_0, \phi_1)$ bounded in $\mathcal{E}_{1, \ga_0}[\phi]$ for some $1<\ga_0<1$ has the following pointwise decay property
  \begin{align*}
    |\phi^{lin}(t, x)|\leq C \sqrt{\mathcal{E}_{1, \ga_0}}(1+t+|x|)^{-1}(1+|t-|x||)^{-\frac{\ga_0-1}{2}}
  \end{align*}
  for some universal constant $C$. Thus when $\frac{1+\sqrt{17}}{2}<p<5$, for arbitrary large data $(\phi_0, \phi_1)$ bounded in $\mathcal{E}_{1, \ga_0}[\phi]$, the solution to the nonlinear equation \eqref{eq:NLW:semi:3d} decays as quickly as the solution to the linear equation with the same initial data. This pointwise decay property is consistent with the scattering result obtained in the author's companion paper \cite{yang:scattering:NLW}, in which it has shown that the solution to \eqref{eq:NLW:semi:3d} scatters in critical Sobolev space $\dot{H}^{\frac{3}{2}-\frac{2}{p-1}}$ and the energy space $\dot{H}^{1}$ when $\frac{1+\sqrt{17}}{2}<p<5$.
\end{Remark}
\begin{Remark}
  Our scattering result in energy space applies to power even below the Strauss exponent $p_c=1+\sqrt{2}$, for which small data global solution and scattering hold for the pure power semilinear wave equation with power above $p_c$ (see for example \cite{Pecher88:scattering:sharpp:3D}) while finite time blow up can occur with power below $p_c$ (see John's work in \cite{John79:blowup:NLW:3d}).
\end{Remark}

As mentioned above, the existing approach (see for example \cite{Pecher82:decay:3d}, \cite{Velo87:decay:NLW}) to study the asymptotic behavior of solutions to \eqref{eq:NLW:semi:3d} relied on the following time decay of the potential energy
\begin{align}
\label{eq:timedecay:3D:Pecher}
  \int_{\mathbb{R}^3}|\phi|^{p+1}dx\leq C (1+t)^{\max\{4-2p, -2\}},\quad 1<p<5,
\end{align}
which is based on the following energy estimate
\begin{align*}
  \int_{\mathbb{R}^{3}}t^2 |\phi|^{p+1}(t, x)dx+\int_0^{t}\int_{\mathbb{R}^3}(2p-6)s|\phi|^{p+1}(s, x)dxds\leq C\mathcal{E}_{0, 2}[\phi]
\end{align*}
 obtained by using the conformal Killing vector field $t^2\pa_t+r^2 \pa_r$ ($r=|x|$) as multiplier. Here the constant $C$ depends only on $p$. With this a priori decay estimate for the solution, a type of $L^q$ estimate for linear wave equation (prototype of Strichartz estimate, see for example \cite{Brenner75:Lp:LW}) yields the pointwise decay estimate for the solution. This approach only makes use of the time decay of the solution. However, it is well known that linear waves have improved decay away from light cone, which can be quantified by decay in $u=t-|x|$. Our improvement comes from thoroughly utilizing such $u$ decay of linear waves.

 The method we used to explore this $u$ decay is the vector field method originally introduced by Dafermos-Rodnianski in \cite{newapp}. The new ingredient is the $r$-weighted energy estimate derived by using the vector field $r^{\ga}(\pa_t+\pa_r)$ as multiplier with $0\leq \ga\leq 2$.
  Applying this to equation \eqref{eq:NLW:semi:3d}, we obtain that
  \begin{align*}
    \iint_{\mathbb{R}^{3+1}}\frac{p-1-\ga}{p+1}r^{\ga-1}|\phi|^{p+1}dxdt\leq C \mathcal{E}_{0, \ga}[\phi].
  \end{align*}
  See details in \cite{yang:scattering:NLW}. To obtain a useful estimate for the solution, we require that $0<\ga<p-1$. On the other hand, combined with an integrated local energy estimate obtained by using the vector field $f(r)\pa_r$ as multiplier and the classical energy conservation, the energy flux through the outgoing null hypersurface $\H_u$ (constant $u$ hyersurface) decays in terms of $u$. In particular,
  \begin{align*}
    \int_{\H_u} |\phi|^{p+1}d\sigma \leq C(1+|u|)^{-\ga}\mathcal{E}_{0, \ga}[\phi].
  \end{align*}
Integrating in terms of $u$, we then get that
\begin{align*}
   \iint_{\mathbb{R}^{3+1}} (1+|u|)^{\ga-1-\ep}|\phi|^{p+1}dxdt\leq C\mathcal{E}_{0, \ga}[\phi],\quad \forall \ep>0
\end{align*}
by assuming that $\ga>1$ (this forces that $p>2$). This together with the above $r$-weighted energy estimate leads to the spacetime bound
\begin{align*}
  \iint_{\mathbb{R}^{3+1}}(1+t+|x|)^{\ga-1-\ep}|\phi|^{p+1}dxdt\leq C\mathcal{E}_{0, \ga}[\phi],
\end{align*}
which is one of the main results obtained in \cite{yang:scattering:NLW} as restated precisely in the following Proposition \ref{prop:spacetime:bd}. For the subconformal case $p<3$, since $\ga$ can be as large as $p-1$, in terms of time decay, this spacetime bound is stronger than \eqref{eq:timedecay:3D:Pecher} as $p-2>2p-5$. Our improvement on the asymptotic decay properties of the solution heavily relies on this uniform spacetime bound.

To show the pointwise decay estimate of solution to \eqref{eq:NLW:semi:3d}, we start by obtaining a uniform weighted energy flux bound through the backward light cone $\N^{-}(q)$ emanating from the point $q=(t_0, x_0)$. Consider the vector field
\begin{align*}
  X=u_+^{\ga}(\pa_t-\pa_r)+v_+^{\ga}(\pa_t+\pa_r),\quad v_+=\sqrt{1+(t+|x|)^2},\quad u_+=\sqrt{1+u^2}.
\end{align*}
The case when $\ga=2$ corresponds to the conformal Killing vector field while the case when $1<\ga<2$ has been widely used (see for examples
\cite{Dafermos17:C0Kerr}, \cite{LindbladMKG}). Applying this vector field as multiplier to the region bounded by the backward light cone $\N^{-}(q)$, we obtain that
\begin{align*}
  \int_{\mathcal{N}^{-}(q)}\big((1+\frac{x\cdot(x-x_0)}{|x||x-x_0|})v_+^{\ga}+u_+^{\ga}\big)|\phi|^{p+1}
\leq C \mathcal{E}_{0, \ga_0}[\phi],\quad \ga<\ga_0
\end{align*}
 for which the above uniform spacetime bound plays the role that it controls the spacetime integral without a definite sign (see details in Proposition \ref{prop:EF:cone:NW:3d}).
 Once we have this uniform potential energy bound, we apply the representation formula to demonstrate the pointwise decay for the solution. The nonlinear term can be estimated by interpolation between the $L^\infty$ estimate of the solution and the above potential energy bound.
When $p>\frac{1+\sqrt{17}}{2}$ is sufficiently large, it turns out that the coefficient of the $L^\infty$ norm of the solution is integrable from $0$ to $t_0$. Thus Gronwall's inequality leads to the decay properties of the solution. On the other hand when $2<p\leq \frac{1+\sqrt{17}}{2}$ (the lower bound for $p$ comes from the fact that we need $\ga>1$), we split the integral of the nonlinear term into the region close to the tip point $q$ and the region far away. The argument for the region close to $q$ is the same as the case when $p>\frac{1+\sqrt{17}}{2}$ due to the fact that the coefficient is still integrable on a small interval. The integral on the region far away can be bounded directly by using the above uniform potential energy bound, which however losses decay.

The above argument only works in the exterior region $\{t+1\leq |x|\}$. Similar argument with minor modifications also applies to the interior region $\{t+1>|x|\}$ after conformal transformation (for simplicity we only consider the solution in the future $t\geq 0$). Pick up the hyperboloid $\mathbb{H}$ passing through the 2-sphere $\{t=0, |x|=1\}$. The region enclosed by $\mathbb{H}$ contains the interior region and is conformally equivalent to the compact backward cone $\{\tilde{t}+|\tilde{x}|\leq \frac{5}{6}\}$. The study of the asymptotic behavior of solution to \eqref{eq:NLW:semi:3d} in the interior region is then reduced to control the growth of solution to a class of semilinear wave equation on this compact cone with initial data determined by the original solution on the hyperboloid $\mathbb{H}$, which has already been understood. The argument to control the solution on the compact region is similar to that used for studying the solution in the exterior region.

The plan of the paper is as follows: In Section 2, we define some notations. In Section 3, we use vector field method to derive a uniform weighted energy estimate for the solution through backward light cones, based on which we obtain the pointwise decay estimate for the solution in the exterior region. In addition, we show quantitative and necessary properties for the solution on the hyperboloid $\mathbb{H}$, used as initial data for the solution in the interior region. In section 4, we study a class of semilinear wave equation on a compact backward cone. The approach is similar but this section is independent of others. In the last section, we conduct the conformal transformation and apply the result of Section 4 to conclude the pointwise decay estimate for the solution in the interior region.

\section{Preliminaries and notations}
\label{sec:notation}

We use the standard polar local coordinate system $(t, r,
\om)$ of Minkowski space as well as the null coordinates $u=\frac{t-r}{2}$, $v=\frac{t+r}{2}$, in which $\om$ is the coordinate of the unit sphere.
Introduce a null frame $\{L, \Lb, e_1, e_2\}$ with
\[
L=\pa_v=\pa_t+\pa_r,\quad \Lb=\pa_u=\pa_t-\pa_r
\]
and $\{e_1, e_2\}$ an orthonormal basis of the sphere with
constant radius $r$. At any fixed point $(t, x)$, we may choose $e_1$, $e_2$ such that
\begin{equation}
 \label{eq:nullderiv}
 \begin{split}
&\nabla_{e_i}L=r^{-1}e_i,\quad \nabla_{e_i}\Lb=-r^{-1}e_i, \quad \nabla_{e_1}e_2=\nabla_{e_2}e_1=0,\quad \nabla_{e_i}e_i=-r^{-1}\pa_r,
 \end{split}
\end{equation}
where $\nabla$ is the covariant derivative in Minkowski space. Defines the functions
\[
u_+=\sqrt{1+u^2},\quad v_+=\sqrt{1+v^2}.
\]
Through out this paper, the exterior region will be referred as $\{(t, x)| u=\frac{t-|x|}{2}\leq -1, t\geq 0\}$ while the interior region will be $\{(t, x)| u\geq -1, t\geq 0\}$.
Let $\H_u$ be the outgoing null hypersurface $\{t-|x|=2u, |x|\geq 2\}$ and $\Hb_v$ be the incoming null hypersurface $\{t+|x|=2v, |x|\geq 2\}$. In the exterior region, we may also use the truncated ones $\H_u^{v_1}$, $\Hb_{v}^{u_1}$ defined as follows
\begin{align*}
  \H_u^{v_1}=\H_u\cap\{-u\leq v\leq v_1\}, \quad
  \Hb_v^{u_1}=\Hb_v\cap\{-v\leq u\leq u_1\}
\end{align*}
and the domain $\mathcal{D}_{u_1}^{ u_2}$  bounded by $\H_{u_1}^{-u_2}$, $\Hb_{-u_2}^{u_1}$ and the initial hypersurface for all $u_2<u_1\leq -1$.

Additional to the above null hypersurfaces, we will also use the hyperboloid
\begin{equation}
\label{eq:def4Hyperboloid}
\Hy:=\left\{(t, x)|(t^*)^2-|x|^2=2R^* t^*\right\}, \quad t^*=t+3, \quad R^*=\frac{5}{6},
\end{equation}
which splits into the future part $\Hy^{+}=\Hy\cap \{t\geq 0\}$ and the past part $\Hy^{-}=\Hy\cap \{t<0\}$. We may note here that the interior region defined above lies inside this hyperboloid.

For any $q=(t_0, x_0)\in \mathbb{R}\times \mathbb{R}^3$ and $r>0$, denote $\B_q(r)$ as the 3-dimensional ball at time $t_0$ with radius $r$ centered at point $q$, that is,
\begin{align*}
  \B_q( r)=\{(t,x)| t=t_0, |x-x_0|\leq r\}.
\end{align*}
 The boundary of $\B_q( r)$ is the 2-sphere $\S_q(r)$. On the initial hypersurface $\{t=0\}$, we use $\B_{r_1}^{r_2}$ to denote the annulus $\{t=0, r_1\leq |x|\leq r_2\}$.

For any fixed point $q=(t_0, x_0)$, let $\mathcal{N}^{-}(q)$ be the past null cone of the point $q$ (for simplicity we are only concerned for the solution in the future $\{t\geq 0\}$) and $\mathcal{J}^{-}(q)$ to be the past of the point $q$, that is, the region bounded by $\mathcal{N}^{-}(q)$ and the initial hypersurface.
Additional to the standard coordinates $(t, x)$ as well as the associated polar coordinates, let $(\tilde{t}, \tilde{x})$ be the new coordinates centered at the point $q=(t_0, x_0)$
\[
\tilde{t}=t-t_0,\quad \tilde{x}=x-x_0,\quad \tilde{r}=|\tilde{x}|,\quad \tilde{\om}=\frac{\tilde{x}}{|\tilde{x}|},\quad \tilde{u}=\f12 (\tilde{t}-\tilde{r}),\quad \tilde{v}=\f12(\tilde{t}+\tilde{r}).
\]
We also have the associated null frame $\{\tilde{L}, \tilde{\Lb}, \tilde{e}_1, \tilde{e}_2\}$ verifying the same relation \eqref{eq:nullderiv}. Under this new coordinates, the past null cone $\mathcal{N}^{-}(q)$ can be characterized by $\{\tilde{v}=0\}\cap\{0\leq t\leq t_0\}$. Through out this paper, the coordinates $(\tilde{t}, \tilde{x})$ are always referred to be the translated ones centered at the point  $q=(t_0, x_0)$ unless it is clearly emphasized.

For simplicity, for integrals in this paper, we will omit the volume form unless it is specified. More precisely we will use
\begin{align*}
  \int_{\mathcal{D}}f,\quad \int_{\H} f, \quad \int_{\Hb}f, \quad \int_{\{t=constant\}} f
\end{align*}
to be short for
\begin{align*}
  \int_{\mathcal{D}}f dxdt, \quad \int_{\H} f 2r^{2}dv d\om , \quad \int_{\Hb}f 2r^{2}dud\om, \quad \int_{\{t=constant\}} f dx
\end{align*}
respectively. Here $\om$ are the standard coordinates of unit sphere.

Finally we make a convention through out this paper to avoid too many constants that $A\les B$ means that there exists a constant $C$, depending possibly on $p$, $\ga_0$ the weighted energy $\mathcal{E}_{0, \ga_0}[\phi]$ such that $A\leq CB$.

\section{A uniform weighted energy flux bound}

In this section, we establish a uniform weighted energy flux bound on any backward light cone in terms of the zeroth order initial energy, based on the spacetime bound for the solution derived in the author's companion paper \cite{yang:scattering:NLW}, from which we recall the following:
\begin{Prop}
  \label{prop:spacetime:bd}
  For all $2<p\leq 5 $ and $1<\ga_0<\min\{2,  p-1\}$, the solution $\phi$ of \eqref{eq:NLW:semi:3d} is uniformly bounded in the following sense
\begin{align}
  \label{eq:spacetime:bd}
  \iint_{\mathbb{R}^{3+1}} v_+^{\ga_0-\ep-1}|\phi|^{p+1}dxdt \leq C \mathcal{E}_{0, \ga_0}[\phi]
\end{align}
for some constant $C$ relying only on $p$, $\ep$ and $\ga_0$.
\end{Prop}
\begin{proof}
  See the main theorem in \cite{yang:scattering:NLW}.
\end{proof}
Using this spacetime bound, we establish the following uniform weighted energy flux bound.
\begin{Prop}
\label{prop:EF:cone:NW:3d}
Let $q=(t_0, x_0)$ be any point in $\mathbb{R}^{3+1}$ with $t_0\geq 0$. Then for solution $\phi$ of the nonlinear wave equation \eqref{eq:NLW:semi:3d} and for all $1< \ga<\ga_0< \min\{2, p-1\}$, we have the following uniform bound
\begin{equation}
\label{eq:Eflux:ex:EF}
\begin{split}
&\int_{\mathcal{N}^{-}(q)}\big((1+\tau)v_+^{\ga}+u_+^{\ga}\big)|\phi|^{p+1}
\leq C \mathcal{E}_{0, \ga_0}[\phi]
\end{split}
\end{equation}
for some constant $C$ depending only on $p$, $\ga_0$ and $\ga$ and independent of the point $q$. Here $\tau=\om\cdot \tilde{\om}$, $r_0=|x_0|$ and the tilde components are measured under the coordinates $(\tilde{t}, \tilde{x})$ centered at the point $q=(t_0, x_0)$.
\end{Prop}
\begin{proof}
Define the energy momentum tensor for the scalar field $\phi$
\begin{align*}
  T[\phi]_{\mu\nu}=\pa_{\mu}\phi\pa_{\nu}\phi-\f12 m_{\mu\nu}(\pa^\ga \phi \pa_\ga\phi+\frac{2}{p+1} |\phi|^{p+1}),
\end{align*}
where $m_{\mu\nu}$ is the flat Minkowski metric.
For any vector field $X$ and any function $\chi$, define the current
\begin{equation*}
J^{X,  \chi}_\mu[\phi]=T[\phi]_{\mu\nu}X^\nu -
\f12\pa_{\mu}\chi \cdot|\phi|^2 + \f12 \chi\pa_{\mu}|\phi|^2.
\end{equation*}
By using Stokes' formula, we have the energy identity
\begin{equation}
\label{eq:energy:id}
\iint_{\pa\mathcal{D}}i_{  J^{X, \chi}[\phi]} d\vol =\iint_{\mathcal{D}}  T[\phi]^{\mu\nu}\pi^X_{\mu\nu}+
\chi \pa_\mu\phi\pa^\mu\phi -\f12\Box\chi\cdot|\phi|^2 +\chi \phi\Box\phi+X(\phi)(\Box\phi-|\phi|^{p}\phi) d\vol
\end{equation}
for any domain $\mathcal{D}$ in $\mathbb{R}^{3+1}$. Here $\pi^X=\f12 \cL_X m$  is the deformation tensor of the metric $m$ along the vector field $X$ and $i_Z d\vol$ is the contraction of the vector field with the volume form $d\vol$.

For the weighted energy flux estimate \eqref{eq:Eflux:ex:EF}, we choose the vector field $X$ as follows:
\[
X=v_+^{\gamma} L+u_+^\gamma \Lb.
\]
Take the region $\cD$ to be $\mathcal{J}^{-}(q)$ which is bounded by the backward light cone $\N^{-}(q)$ and the initial hypersurface. For the above chosen vector field $X$, we can compute that
\[
\nabla_{L}X=\gamma v_+^{\gamma-2}v L,\quad \nabla_{\Lb}X= \gamma u_+^{\gamma-2}u \Lb,\quad \nabla_{e_i}X=r^{-1}(v_+^\gamma-u_+^{\gamma}) e_i.
\]
Then the non-vanishing components of the deformation tensor $\pi_{\mu\nu}^X$ are
\[
\pi^X_{L\Lb}=-\gamma \left(v_+^{\gamma-2}v+u_+^{\gamma-2}u\right),\quad \pi^X_{e_i e_i}=r^{-1}(v_+^{\ga}-u_+^{\gamma}).
\]
Therefore we can compute that
\begin{align*}
&T[\phi]^{\mu\nu}\pi^X_{\mu\nu}=2T[\phi]^{L\Lb}\pi^X_{L\Lb}+T[\phi]^{e_ie_i}\pi^X_{e_ie_i}\\
&=-\f12\ga (v_+^{\gamma-2}v+u_+^{\gamma-2}u)(|\nabb\phi|^2+\frac{2}{p+1}|\phi|^{p+1})+r^{-1}(v_+^{\ga}-u_+^{\gamma})(|\nabb\phi|^2-\pa^\mu \phi \pa_\mu\phi-\frac{2}{p+1}|\phi|^{p+1})\\
&=\left(-\f12\ga(v_+^{\gamma-2}v+u_+^{\gamma-2}u)+r^{-1}(v_+^{\ga}-u_+^{\gamma})\right)|\nabb\phi|^2-r^{-1}(v_+^{\ga}-u_+^{\gamma})\pa^\mu\phi \pa_\mu\phi\\
&\quad +\left(-\f12\ga(v_+^{\gamma-2}v+u_+^{\gamma-2}u)-r^{-1}(v_+^{\ga}-u_+^{\gamma})\right)\frac{2}{p+1}|\phi|^{p+1}.
\end{align*}
Now choose the function $\chi$ as follows:
\[
\chi=r^{-1}(v_+^{\ga}-u_+^{\gamma}).
\]
For such spherically symmetric function $\chi$ (with respect to the coordinates $(t, x)$), we can compute that
\begin{align*}
\Box \chi=-r^{-1}L\Lb (r\chi)=-2r^{-1}L\Lb(v_+^\ga-u_+^\ga)=0,\quad r>0.
\end{align*}
At $r=0$ it grows at most $r^{\gamma-3}$. Therefore we can write that
\begin{align*}
&T[\phi]^{\mu\nu}\pi^X_{\mu\nu}+
\chi \pa_\mu\phi \pa^\mu\phi -\f12\Box\chi\cdot|\phi|^2+\chi\phi\Box\phi  \\
&=\left(\chi-\f12\ga(v_+^{\gamma-2}v+u_+^{\gamma-2}u)\right)|\nabb\phi|^2 + \left(\frac{p-1}{p+1}\chi-\frac{(v_+^{\ga-2}v+u_+^{\ga-2}u)\ga}{p+1}\right)|\phi|^{p+1}.
\end{align*}
Denote $f(s)=(1+s^2)^{\frac{\ga}{2}}$. Then $\chi=\frac{f(v)-f(u)}{v-u}$. It can be checked directly that the derivative $f'(s)$ is concave. In particular we conclude that
\begin{align*}
  \chi=\frac{f(v)-f(u)}{v-u}\geq \f12 (f'(v)+f'(u))=\f12 \ga(v_+^{\ga-2}v+u_+^{\ga-2}u).
\end{align*}
Therefore the coefficient of $|\nabb\phi|^2$ is nonnegative. On the other hand, the coefficient of $|\phi|^{p+1}$ can be trivially bounded by
\begin{align*}
  |\frac{p-1}{p+1}\chi-\frac{(v_+^{\ga-2}v+u_+^{\ga-2}u)\ga}{p+1}|\leq Cv_+^{\ga-1}
\end{align*}
for some constant $C$ depending only on $p$ and $\ga$. We remark here that this coefficient is also nonnegative for the super-conformal case when $p\geq 3$. We use Proposition \ref{prop:spacetime:bd} to control this potential term for the sub-conformal case when the sign is indefinite.

We next compute the boundary integrals on the left hand side of the energy identity \eqref{eq:energy:id}, which consists of the integral on the initial hypersurface $\B_{(0, x_0)}(t_0)$ and the integral on the backward light cone $\mathcal{N}^{-}(q)$. Let's first compute the boundary integral on the initial hypersurface, under the coordinates system $(t, x)$. As the initial hypersurface $\B_{(0, x_0)}(t_0)$ has the volume form $dx$, the contraction reads
\begin{align*}
i_{J^{X, \chi}[\phi ]}d\vol&= T[\phi ]_{0 L}X^L+T[\phi ]_{\Lb 0}X^{\Lb}- \f12 \pa_t\chi |\phi|^2+ \f12\chi \pa_t|\phi|^2\\
&= \f12 v_+^\ga( |L\phi|^2+|\nabb\phi|^2+\frac{2}{p+1}|\phi|^{p+1})-\f12\pa_t\chi \cdot |\phi|^2+\f12\chi \pa_t|\phi|^2 \\
&\quad +\f12 u_+^\ga(|{\Lb}\phi|^2+|\nabb\phi|^2+\frac{2}{p+1}|\phi|^{p+1})\\
&=\f12(u_+^\ga+v_+^\ga)( |\nabb\phi|^2+\frac{2}{p+1}|\phi|^{p+1})+\f12 v_+^\ga r^{-2}|L(r\phi)|^2\\
&\quad +\f12 u_+^\ga r^{-2}|{\Lb}(r\phi)|^2-\div(\om r^{-1}|\phi|^2(u_+^\ga+v_+^\ga)).
\end{align*}
Here $\om=\frac{x}{|x|}$ can be viewed as a vector on $\mathbb{R}^{3}$ and the divergence is taken over the initial hypersurface $\B_{(0, x_0)}(t_0)$. The integral of the divergence term and be computed by using integration by parts. Under the coordinates $\tilde{x}=x-x_0$ on the initial hypersurface, we have
\begin{align*}
\int_{\B_{(0, x_0)}(t_0)} \div (\om r^{-1}|\phi|^2(u_+^\ga+v_+^\ga))dx&=\int_{\B_{(0, x_0)}(t_0)} \div (\om r^{-1}|\phi|^2(u_+^\ga+v_+^\ga))d\tilde{x}\\
&=\int_{\S_{(0, x_0)}(t_0)} \tilde{r}^2 \tilde{\om} \cdot\om r^{-1}|\phi|^2(u_+^\ga+v_+^\ga)d\tilde{\om}.
\end{align*}
In particular we derive that
\begin{equation}
\label{eq:PWE:ex:bxt0}
\begin{split}
 &\int_{\B_{(0, x_0)}(t_0)} i_{J^{X,\chi}[\phi]}d\vol + \int_{\S_{(0, x_0)}(t_0)} \tilde{r}^2 \tilde{\om} \cdot\om r^{-1}|\phi|^2(u_+^\ga+v_+^\ga)d\tilde{\om}\\
&=\f12 \int_{\B_{(0, x_0)}(t_0)}     v_+^{\ga} r^{-2}|L(r\phi)|^2 + (u_+^{\ga}+v_+^{\ga})(|\nabb\phi|^2+\frac{2}{p+1}|\phi|^{p+1})+u_+^{\ga}r^{-2}|\Lb(r\phi)|^2 dx \\
 &\leq C \mathcal{E}_{0, \ga}[\phi]
 \end{split}
\end{equation}
for some constant $C$ depending only on $\ga$.

For the boundary integral on the backward light cone $\mathcal{N}^{-}(q)$, we shift to the coordinates centered at the point $q=(t_0, x_0)$.
 Recall the volume form
\[
d\vol=dxdt=d\tilde{x}d\tilde{t}=2\tilde{r}^2 d\tilde{v}d\tilde{u}d\tilde{\om}.
\]
Since the backward light cone $\mathcal{N}^{-}(q)$  can be characterized by $\{\tilde{v}=0\}$ under these new coordinates $(\tilde{t}, \tilde{x})$, we therefore have
\begin{align*}
-i_{J^{X, \chi}[\phi]}d\vol=J_{\tilde{\Lb}}^{X, \chi}[\phi]\tilde{r}^2d\tilde{u}d\tilde{\om}= ( T[\phi]_{\tilde{\Lb}\nu}X^\nu -
\f12(\tilde{\Lb}\chi) |\phi|^2 + \f12 \chi\cdot\tilde{\Lb}|\phi|^2 ) \tilde{r}^2d\tilde{u}d\tilde{\om}.
\end{align*}
For the main quadratic terms, we first can compute that
\begin{align*}
T[\phi]_{\tilde{\Lb}\nu}X^\nu =T[\phi]_{\tilde{\Lb}\tilde{\Lb}}X^{\tilde{\Lb}}+T[\phi]_{\tilde{\Lb}\tilde{L}}X^{\tilde{L}}+T[\phi]_{\tilde{\Lb}\tilde{e}_i}X^{\tilde{e}_i}.
\end{align*}
Since the vector field $X$ is given under the coordinates $(t, x)$, we need to write it under the new null frame $\{\tilde{L}, \tilde{\Lb}, \tilde{e}_1, \tilde{e}_2\}$ centered at the point $q$. Note that
\begin{align*}
\pa_r=\om \cdot \nabla=\om \cdot \tilde{\nabla}=\om\cdot \tilde{\om}\pa_{\tilde{r}}+ \om\cdot (\tilde{\nabla}-\tilde{\om}\pa_{\tilde{r}}).
\end{align*}
Here $\om=\frac{x}{|x|}$, $\nabla=(\pa_{x^1}, \pa_{x^2}, \pa_{x^3})$. Thus we can write that
\begin{align*}
X &=(v_+^\ga+u_+^\ga)\pa_{t}+(v_+^\ga-u_+^\ga)\pa_r\\
&=(v_+^\ga+u_+^\ga)\pa_{\tilde{t}}+(v_+^\ga-u_+^\ga)(\om\cdot \tilde{\om}\pa_{\tilde{r}}+ \om\cdot \tilde{\nabb})\\
&=\f12 \left(u_+^\ga+v_+^\ga+(v_+^\ga-u_+^\ga)\om\cdot \tilde{\om}\right)\tilde{L}+\f12 \left(u_+^\ga+v_+^\ga-(v_+^\ga-u_+^\ga)\om\cdot \tilde{\om}\right)\tilde{\Lb}+(v_+^\ga-u_+^\ga)\om\cdot \tilde{\nabb}.
\end{align*}
Here $\tilde{\nabb}=\tilde{\nabla}-\pa_{\tilde{r}}$. Thus we can compute the quadratic terms
\begin{align*}
T[\phi]_{\tilde{\Lb}\nu}X^\nu
=&\left((1-\tau)v_+^\ga+(1+\tau)u_+^\ga\right)|{\tilde{\Lb}}\phi|^2 +\left((1+\tau)v_+^\ga+(1-\tau)u_+^\ga\right)(|\tilde{\nabb}\phi|^2+\frac{2}{p+1}|\phi|^{p+1})\\
&+2 (v_+^\ga-u_+^\ga) ({\tilde{\Lb}}\phi) (\om\cdot \tilde{\nabb})\phi.
\end{align*}
 Here recall that $\tau=\om\cdot \tilde{\om}$.
 It turns out that these terms are nonnegative but we need to estimate them together with the lower order terms arising from the function $\chi$.
We compute that
\begin{align*}
&\tilde{\Lb}(r)=-\pa_{\tilde{r}}(r)=-\tilde{\om}_i\pa_i(r)=-\tilde{\om}\cdot \om =-\tau,\\
&\tilde{\nabb}(r)=(\tilde{\nabla}-\tilde{\om}\pa_{\tilde{r}})(r)=\om-\tilde{\om}\tau.
\end{align*}
Therefore we can write
\begin{align*}
&-\f12 r^2 (\tilde{\Lb}{\chi})|\phi|^2+\f12 r^2\chi \tilde{\Lb}|\phi|^2=(r\chi)( {\tilde{\Lb}}(r\phi)+\tau \phi) \phi-\f12(r\tilde{\Lb}(r\chi)+\tau r\chi)|\phi|^2,\\
&r^2|\tilde{\Lb}\phi|^2=|{\tilde{\Lb}}(r\phi)-\tilde{\Lb}(r)\phi|^2=|{\tilde{\Lb}}(r\phi)|^2+\tau^2|\phi|^2+2 {\tilde{\Lb}}(r\phi) \tau\phi,\\
& r^2|\tilde{\nabb}\phi|^
=|\tilde{\nabb}(r\phi)|^2+(1-\tau^2)|\phi|^2-2(\om-\tilde{\om}\tau)\cdot\tilde{\nabb}(r\phi) \phi,\\
& r^2 ({\tilde{\Lb}}\phi)  (\om\cdot \tilde{\nabb})\phi={\tilde{\Lb}}(r\phi) (\om \cdot \tilde{\nabb})(r\phi)-\tau(1-\tau^2)|\phi|^2+\phi \tau(\om\cdot \tilde{\nabb})(r\phi) -(1-\tau^2){\tilde{\Lb}}(r\phi)\phi.
\end{align*}
Notice that
\[
\om\cdot \tilde{\nabb}=\om\cdot (\tilde{\om}\times \tilde{\nabla})=\om\times \tilde{\om}\cdot \tilde{\nabb}.
\]
Since $v_+\geq u_+$, we therefore can show that the quadratic terms are nonnegative
\begin{align*}
&\left((1-\tau)v_+^\ga+(1+\tau)u_+^\ga\right)|\tilde{\Lb}(r\phi)|^2+\left((1+\tau)v_+^\ga+(1-\tau)u_+^\ga\right)|\tilde{\nabb}(r\phi)|^2+2 (v_+^\ga-u_+^\ga){\tilde{\Lb}}(r\phi) (\om \cdot \tilde{\nabb})(r\phi)\\
&\geq 2\sqrt{(v_+^{\ga}+u_+^{\ga})^2-\tau^2(v_+^{\ga}-u_+^{\ga})^2}|\tilde{\nabb}(r\phi)||\tilde{\Lb}(r\phi)|-2 (v_+^\ga-u_+^\ga)\sqrt{1-\tau^2}|{\tilde{\Lb}}(r\phi)| |\tilde{\nabb}(r\phi)|\\
& \geq 0.
\end{align*}
For the other lower order terms, we compute  that
\begin{align*}
&\left((1-\tau)v_+^\ga+(1+\tau)u_+^\ga\right)(\tau^2|\phi|^2+2{\tilde{\Lb}}(r\phi)\tau\phi )+(r\chi)( {\tilde{\Lb}}(r\phi)+\tau \phi) \phi-\f12(r\tilde{\Lb}(r\chi)+\tau r\chi)|\phi|^2\\
&+\left((1+\tau)v_+^\ga+(1-\tau)u_+^\ga\right)\big((1-\tau^2)|\phi|^2-2(\om-\tilde{\om}\tau)\tilde{\nabb}(r\phi) \phi\big)\\
&+(v_+^\ga-u_+^\ga)\big(-2\tau(1-\tau^2)|\phi|^2+2\tau \phi (\om\cdot \tilde{\nabb})(r\phi)-2\phi (1-\tau^2){\tilde{\Lb}}(r\phi) \big)\\
&=(-\f12 r\tilde{\Lb}(r\chi)+v_+^\ga+u_+^\ga)|\phi|^
+2(v_+^\ga+u_+^\ga)(\tau {\tilde{\Lb}}-\om\cdot \tilde{\nabb})(r\phi) \phi\\
&=-r^2\tilde{r}^{-1} \tilde{\Om}_{ij}(r^{-3}(v_+^\ga+u_+^\ga) \om_j\tilde{\om}_i |r\phi|^2)+\tilde{r}^{-2}r^2\tilde{\Lb}(r^{-1}\tau\tilde{r}^2(v_+^\ga+u_+^\ga) |\phi|^2)\\
&+(-\f12 r\tilde{\Lb}(r\chi)+v_+^\ga+u_+^\ga)|\phi|^2-\tilde{r}^{-2}r^2\tilde{\Lb}(r^{-3}\tau\tilde{r}^2(v_+^\ga+u_+^\ga)) |r\phi|^2+r^2 \tilde{r}^{-1}\tilde{\Om}_{ij}(r^{-3}(v_+^\ga+u_+^\ga) \om_j\tilde{\om}_i) |r\phi|^2.
\end{align*}
We can compute that
\begin{align*}
&\tilde{r}^{-1}\tilde{\Om}_{ij}(r^{-3}\om_j\tilde{\om}_i)=-2r^{-4}(1-2\tau^2)-2\tau \tilde{r}^{-1}r^{-3},\\
&\tilde{r}^{-2}r^4\tilde{\Lb}(r^{-3}\tilde{r}^2\tau)=4\tau^2-1-2r\tilde{r}^{-1}\tau.
\end{align*}
Thus the coefficients of $|\phi|^2$ in the last line in the previous equation verify
\begin{align*}
&(-\f12 r\tilde{\Lb}(r\chi)+v_+^\ga+u_+^\ga)-\tilde{r}^{-2}r^4\tilde{\Lb}(r^{-3}\tau\tilde{r}^2(v_+^\ga+u_+^\ga)) +r^4 \tilde{r}^{-1} \tilde{\Om}_{ij}(r^{-3}(v_+^\ga+u_+^\ga) \om_j\tilde{\om}_i) \\
&=-r(\pa_t-\tilde{\om}\cdot \nabla)(v_+^\ga-u_+^\ga)-r\tau (\pa_t-\tilde{\om}\cdot \nabla)(u_+^\ga+v_+^\ga)+r(\pa_r-\tau \tilde{\om}\cdot \nabla)(u_+^\ga+v_+^\ga)\\
&\quad +(u_+^\ga+v_+^\ga)\left(1-(4\tau^2-1-2r\tilde{r}^{-1}\tau)-2(1-2\tau^2)+2\tau \tilde{r}^{-1}r\right)\\
&=r(\pa_t+\pa_r)u_+^\ga+r(\pa_r-\pa_t)v_+^\ga-\tau r(\pa_t+\pa_r)u_+^\ga-\tau r(\pa_t-\pa_r)v_+^\ga=0.
\end{align*}
The above computations imply that the lower order terms can be written as a divergence form and hence can be estimated by using integration by parts:
\begin{align*}
  &\int_{\mathcal{N}^{-}(q)}\big(-r^2\tilde{r}^{-1} \tilde{\Om}_{ij}(r^{-3}(v_+^\ga+u_+^\ga) \om_j\tilde{\om}_i |r\phi|^2)+\tilde{r}^{-2}r^2\tilde{\Lb}(r^{-1}\tau\tilde{r}^2(v_+^\ga+u_+^\ga) |\phi|^2)\big)r^{-2}\tilde{r}^2 d\tilde{u}d\tilde{\om}\\
  &= \int_{\S_{(0, x_0)}(t_0)}r^{-1}\tau \tilde{r}^2 (u_+^\ga+v_+^\ga)|\phi|^2d\tilde{\om}.
\end{align*}
This term is an integral on the sphere on the initial hypersurface and cancels the one arising from the boundary integral on $\B_{(0, x_0)}(t_0)$.
Keeping the potential part and discarding the quadratic terms which are nonnegative, we in therefore derive that
\begin{equation*}
\begin{split}
&\frac{2}{p+1}\int_{\mathcal{N}^{-}(q)}((1+\tau)v_+^\ga+(1-\tau)u_+^\ga ) |\phi|^{p+1} \tilde{r}^2d\tilde{u}d\tilde{\om}\\
&\leq -\int_{\mathcal{N}^{-}(q)}i_{J^{X,\chi}[\phi]}d\vol+ \int_{\S_{(0, x_0)}(t_0)}r^{-1}\tau \tilde{r}^2 (u_+^\ga+v_+^\ga)|\phi|^2d\tilde{\om}.
\end{split}
\end{equation*}
Combining this estimate with \eqref{eq:PWE:ex:bxt0} and by using the uniform spacetime bound of Proposition \ref{prop:spacetime:bd}, we then derive that
\begin{align*}
  \frac{2}{p+1}\int_{\mathcal{N}^{-}(q)}((1+\tau)v_+^\ga+(1-\tau)u_+^\ga ) |\phi|^{p+1}
 & \leq \iint_{\mathcal{J}^{-}(q)}\left|\frac{p-1}{p+1}\chi-\frac{(v_+^{\ga-2}v+u_+^{\ga-2}u)\ga}{p+1}\right||\phi|^{p+1}\\
  &\quad +\int_{\pa\mathcal{J}^{-}(q)}i_{J^{X, \chi}[\phi]}d\vol\\
  &\leq C \mathcal{E}_{0, \ga+\ep}[\phi]
\end{align*}
for some constant $C$ depending only on $\ga$, $p$ and $0<\ep<p-1-\ga$. The proposition then follows by letting $0<\ep<\ga_0-\ga$.
\end{proof}

\section{The pointwise decay of the solution in the exterior region}
In this section, we make use of the weighted energy flux bound derived in the previous section to investigate the asymptotic behaviour of the solution in the exterior region $\{t+2\leq |x|\}$.

We need the following integration lemma.
\begin{Lem}
\label{lem:integration:ex:ab}
  Assume $1<\ga<2$ and $\a$, $\b$ nonnegative such that $\b+\a\ga>2$. Fix $q=(t_0, x_0)$ in the exterior region. For the $2$-sphere $\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})$ on the backward light cone $\mathcal{N}^{-}(q)$, we have
  \begin{equation}
    \label{eq:integration:ex:ab}
    \begin{split}
    &\int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})} ((1+\tau)r^{\ga}+(r_0-t_0)^\ga)^{-\a}r^{-\b}d\tilde{\om}\\
    &\leq C (r_0-\tilde{r})^{2-\b-\ga+\ep} r_0^{-2}\left((r_0-\tilde{r})^{(1-\a)\ga}+(r_0-t_0)^{(1-\a)\ga}\right)
    \end{split}
  \end{equation}
  for some constant $C$ depending only on $\ep$, $\ga$, $\a$ and $\b$.
  Here $\tau=\om\cdot \tilde{\om}$, $r_0=|x_0|$ and $0\leq \tilde{r}\leq t_0<r_0$. And the small positive constant $\ep$ only appears for the case when $\a=1$.
\end{Lem}
\begin{proof}
Denote $s=-\om_0\cdot \tilde{\om}$ with $\om_0=r_0^{-1}x_0$. Note that
  \begin{align*}
    &r^2=|x_0+\tilde{x}|^2=\tilde{r}^2+r_0^2+2r_0\tilde{r}\om_0\cdot \tilde{\om}=(\tilde{r}-r_0s)^2+(1-s^2)r_0^2,\\
    &(1+\tau)r=r+r\om\cdot \tilde{\om}=r+(\tilde{x}+x_0)\cdot \tilde{\om}=r+\tilde{r}-r_0s.
  \end{align*}
We can write the integral as
  \begin{align*}
    &\int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})} ((1+\tau)r^{\ga}+(r_0-t_0)^\ga)^{-\a}r^{-\b}d\tilde{\om} =2\pi\int_{-1}^1 r^{-\b}(r^{\ga-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga})^{-\a}ds.
  \end{align*}
  Without loss of generality, we may assume that $t_0\geq \frac{9}{10}r_0$. Otherwise it trivially has
   \begin{align*}
    \int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})} ((1+\tau)r^{\ga}+(r_0-t_0)^\ga)^{-\a}r^{-\b}d\tilde{\om}
    \leq 4\pi 10^{\b+\a\ga} r_0^{-\b-\a\ga}
  \end{align*}
   since $r\geq r_0-t_0\geq \frac{1}{10}r_0$.

For the integral on $s\leq 0$, we trivially bound that
  \begin{align*}
    r\geq r_0,\quad r+\tilde{r}-r_0s\geq r+\tilde{r}\geq r_0.
  \end{align*}
  Thus
  \begin{align*}
    \int_{-1}^0 r^{-\b}(r^{\ga-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga})^{-\a}ds\leq r_0^{-\b-\a\ga}.
  \end{align*}
  Define $s_0=1-(1-\tilde{r}r_0^{-1})^2$. On the interval $[0, s_0]$, note that
  \begin{align*}
    \sqrt{1-s} \ r_0\geq r_0-\tilde{r}.
  \end{align*}
  Therefore, we can show that
  \begin{align*}
    \tilde{r}-r_0s \leq r_0(1-s)\leq r_0\sqrt{1-s},\quad r_0s-\tilde{r}\leq r_0-\tilde{r}\leq r_0\sqrt{1-s}
  \end{align*}
  as $\tilde{r}\leq t_0<r_0$. This in particular implies that
  \begin{align*}
    \sqrt{1-s}\ r_0\leq r\leq \sqrt{2(1-s)}\ r_0,\quad \sqrt{(\tilde{r}-r_0s)^2+(1-s^2)r_0^2}+\tilde{r}-r_0s\geq \frac{1}{3}\sqrt{1-s}r_0.
  \end{align*}
  Here the second inequality follows from the inequality
  \[
  \sqrt{a^2+b^2}+b\geq (\sqrt{2}-1)|a|,\quad \forall |b|\leq |a|
  \]
  together with the bound $|\tilde{r}-r_0s|\leq \sqrt{1-s}r_0\leq \sqrt{1-s^2}r_0$.

  Therefore on the interval $[0, s_0]$, we can estimate that
  \begin{align*}
    \int_{0}^{s_0} r^{-\b}(r^{\ga-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga})^{-\b} ds
    &\leq 3^{\a} \int_{0}^{s_0}(\sqrt{1-s} r_0)^{-\b-\a\ga} ds\\
    &\leq \frac{2 \times 3^{\a} }{\b+\a\ga-2}r_0^{-2}(r_0-\tilde{r})^{2-\b-\a\ga}.
  \end{align*}
  Here we used the assumption $\b+\a\ga> 2$.

  Finally on the interval $[s_0, 1]$, notice that
  \begin{align*}
    2r\geq r_0s-\tilde{r}+\sqrt{1-s}r_0=r_0-\tilde{r}+(\sqrt{1-s}-(1-s))r_0\geq r_0-\tilde{r}.
  \end{align*}
  Moreover
  \begin{align*}
    r+\tilde{r}-r_0s=\frac{(1-s^2)r_0^2}{r+r_0s-\tilde{r}}\geq\frac{(1-s)r_0^2}{4(r_0-\tilde{r})}.
  \end{align*}
  Therefore we can estimate that
  \begin{align*}
    &\int_{s_0}^{1} r^{-\b}(r^{\ga-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga})^{-\a} ds\\
    &\leq 2^{\b} \int_{s_0}^{1}\big((r_0-t_0)^{\ga}+2^{\ga-3}(r_0-\tilde{r})^{\ga-2}(1-s)r_0^2\big)^{-\a}(r_0-\tilde{r})^{-\b}ds\\
    &= 2^{\b+3-\ga}(\a-1)^{-1}(r_0-\tilde{r})^{2-\b-\ga} r_0^{-2}\left((r_0-t_0)^{(1-\a)\ga}-((r_0-t_0)^{\ga}+2^{\ga-3}(r_0-\tilde{r})^{\ga})^{1-\a}\right)\\
    &\leq C_\ep (r_0-\tilde{r})^{2-\b-\ga+\ep} r_0^{-2}\left((r_0-\tilde{r})^{(1-\a)\ga}+(r_0-t_0)^{(1-\a)\ga}\right)
  \end{align*}
  for some constant $C_{\ep}$ depending only on $\ep$, $\a$, $\b$ and $\ga$.
  The loss of $\ep$ comes from the case when $\a=1$. Since
  \begin{align*}
    r_0^{-\a\ga-\b}\leq (r_0-\tilde{r})^{2-
    \b-\ga+\ep} r_0^{-2}\left((r_0-\tilde{r})^{(1-\a)\ga}+(r_0-t_0)^{(1-\a)\ga}\right)
  \end{align*}
  due to the assumption $\b+\a\ga> 2$, we thus conclude the Lemma.
\end{proof}

\iffalse
The above lemma is sufficient to deal with the case when $p>\frac{1+\sqrt{17}}{2}$ and $(p-1)(\ga_0+1)>4$. The following two specific lemmas are for the case when $2<p\leq \frac{1+\sqrt{17}}{2}$ and $\ga_0\leq p-1$, which in particular implies that $(p-1)(\ga_0+1)\leq 4$.
\begin{Lem}
  Fix $q=(t_0, x_0)$ in the exterior region. For the $2$-sphere $\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})$ on the backward light cone $\mathcal{N}^{-}(q)$, we have the following estimate
  \begin{equation}
    \label{eq:integration:Lemma:Gronwall:p}
    \begin{split}
    &\int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})} ((1+\tau)r^{\ga_0}+(r_0-t_0)^{\ga_0})^{-\frac{p-1}{2}} d\tilde{\om} \leq C r_0^{\frac{1-p}{2}\ga_0}
    \end{split}
  \end{equation}
  for some constant $C$ depending only on $p$ and $\ga_0$.
  Here $\tau=\om\cdot \tilde{\om}$, $r_0=|x_0|$ and $0\leq \tilde{r}\leq \f12 t_0<r_0$.
\end{Lem}
\begin{proof}
Denote $s=-\om_0\cdot \tilde{\om}$ with $\om_0=r_0^{-1}x_0$. Note that
  \begin{align*}
    &r^2=|x_0+\tilde{x}|^2=\tilde{r}^2+r_0^2+2r_0\tilde{r}\om_0\cdot \tilde{\om}=(\tilde{r}-r_0s)^2+(1-s^2)r_0^2,\\
    &(1+\tau)r=r+r\om\cdot \tilde{\om}=r+(\tilde{x}+x_0)\cdot \tilde{\om}=r+\tilde{r}-r_0s.
  \end{align*}
  Without loss of generality, we may assume that $t_0\geq \frac{9}{10}r_0$. Otherwise it trivially has
   \begin{align*}
    \int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})} ((1+\tau)r^{\ga_0}+(r_0-t_0)^{\ga_0})^{-\frac{p-1}{2}} d\tilde{\om}
    \les r_0^{-\frac{p-1}{2}\ga_0}
  \end{align*}
   since $r\geq r_0-t_0\geq \frac{1}{10}r_0$.

    In spherical coordinates, we can write the integral on the sphere as
  \begin{align*}
    &\int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})} ((1+\tau)r^{\ga_0}+(r_0-t_0)^{\ga_0})^{-\frac{p-1}{2}} d\tilde{\om}=2\pi\int_{-1}^1 (r^{\ga_0-1}(r+\tilde{r}-r_0 s)+(r_0-t_0)^{\ga_0})^{-\frac{p-1}{2}}ds.
  \end{align*}
  Since we restricted to the case $\tilde{r}\leq \f12 t_0\leq \f12 r_0$, it always holds that
  \begin{align*}
    r^2=(\tilde{r}-r_0s)^2+(1-s^2)r_0^2\geq \frac{1}{9}r_0^2.
  \end{align*}
  For the integral on $[-1, r_0^{-1}\tilde{r} ]$, we trivially bound that
  \begin{align*}
    \int_{-1}^{r_0^{-1}\tilde{r}} (r^{\ga_0-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga_0})^{-\frac{p-1}{2}}ds\les r_0^{-\frac{p-1}{2}\ga_0}.
  \end{align*}
  On $[r_0^{-1}\tilde{r}, 1]$, we estimate that
  \begin{align*}
    r+\tilde{r}-r_0s=\frac{(1-s^2)r_0^2}{r+r_0s-\tilde{r}}\geq \frac{(1-s)r_0^2}{2r}\geq \frac{(1-s)r_0}{10}.
  \end{align*}
  Therefore we can bound that
  \begin{align*}
    &\int_{r_0^{-1}\tilde{r}}^1 (r^{\ga_0-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga_0})^{-\frac{p-1}{2}}ds\\
    &\les \int_{r_0^{-1}\tilde{r}}^{1}(r_0^{\ga_0}(1-s)+(r_0-t_0)^{\ga_0})^{-\frac{p-1}{2}}ds\\
    &\les r_0^{-\ga_0} \left((r_0^{\ga_0-1}(r_0-\tilde{r})+(r_0-t_0)^{\ga_0})^{\frac{3-p}{2}}-(r_0-t_0)^{\ga_0 \frac{3-p}{2}}
    \right)\\
    &\les r_0^{\frac{1-p}{2}\ga_0}.
  \end{align*}
  Here we used the fact that $\tilde{r}\leq \f12 t_0\leq \f12 r_0$ and $p<3$.
\end{proof}

Still we restricted to the situation when $2<p\leq \frac{1+\sqrt{17}}{2}$ and $(p-1)(\ga_0+1)\leq 4$ in the exterior region.
\begin{Lem}
Assume that $p\ga>2$.
  Fix $q=(t_0, x_0)$ in the exterior region. For the $2$-sphere $\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})$ on the backward light cone $\mathcal{N}^{-}(q)$, we have the following estimate
  \begin{equation}
    \label{eq:integration:Lemma:direct:p}
    \begin{split}
    &\int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})} ((1+\tau)r^{\ga}+(r_0-t_0)^\ga)^{-p} d\tilde{\om}\leq C (r_0-\tilde{r})^{2-\ga} r_0^{-2} (r_0-t_0)^{(1-p)\ga }
    \end{split}
  \end{equation}
  for some constant $C$ depending only on $p$ and $\ga$.
  Here $\tau=\om\cdot \tilde{\om}$, $r_0=|x_0|$ and $0\leq \tilde{r}\leq t_0<r_0$.
\end{Lem}
\begin{proof}
Denote $s=-\om_0\cdot \tilde{\om}$ with $\om_0=r_0^{-1}x_0$. Note that
  \begin{align*}
    &r^2=|x_0+\tilde{x}|^2=\tilde{r}^2+r_0^2+2r_0\tilde{r}\om_0\cdot \tilde{\om}=(\tilde{r}-r_0s)^2+(1-s^2)r_0^2,\\
    &(1+\tau)r=r+r\om\cdot \tilde{\om}=r+(\tilde{x}+x_0)\cdot \tilde{\om}=r+\tilde{r}-r_0s.
  \end{align*}
  We can express the integral as follows
  \begin{align*}
    &\int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})} ((1+\tau)r^{\ga}+(r_0-t_0)^\ga)^{-p}d\tilde{\om}=2\pi\int_{-1}^1 (r^{\ga-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga})^{-p}ds.
  \end{align*}
  For similar reason as the previous lemma, the inequality holds trivially when $t_0\leq \frac{9}{10}r_0$ or restricted to the interval $s\leq 0$. It suffices to assume that $t_0\geq \frac{9}{10}r_0$ and $s\geq 0$.

  Define $s_0=1-(1-\tilde{r}r_0^{-1})^2$. On the interval $[0, s_0]$, note that
  \begin{align*}
    \sqrt{1-s} \ r_0\geq r_0-\tilde{r}.
  \end{align*}
  Therefore, we can show that
  \begin{align*}
    \tilde{r}-r_0s \leq r_0(1-s)\leq r_0\sqrt{1-s},\quad r_0s-\tilde{r}\leq r_0-\tilde{r}\leq r_0\sqrt{1-s}
  \end{align*}
  as $\tilde{r}\leq t_0<r_0$. This in particular implies that
  \begin{align*}
    \sqrt{1-s}\ r_0\leq r\leq \sqrt{2(1-s)}\ r_0,\quad \sqrt{(\tilde{r}-r_0s)^2+(1-s^2)r_0^2}+\tilde{r}-r_0s\geq \frac{1}{3}\sqrt{1-s}r_0.
  \end{align*}
  Therefore on the interval $[0, s_0]$, we can estimate that
  \begin{align*}
    \int_{0}^{s_0} (r^{\ga-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga})^{-p} ds
    &\les \int_{0}^{s_0}(\sqrt{1-s} r_0)^{-p\ga} ds\\
    &\les  r_0^{-2} (r_0-\tilde{r})^{2-p\ga}.
  \end{align*}
  Here we used the assumption $p \ga> 2$.

  Finally on the interval $[s_0, 1]$, notice that
  \begin{align*}
    2r\geq r_0s-\tilde{r}+\sqrt{1-s}r_0=r_0-\tilde{r}+(\sqrt{1-s}-(1-s))r_0\geq r_0-\tilde{r}.
  \end{align*}
  Moreover
  \begin{align*}
    r+\tilde{r}-r_0s=\frac{(1-s^2)r_0^2}{r+r_0s-\tilde{r}}\geq\frac{(1-s)r_0^2}{4(r_0-\tilde{r})}.
  \end{align*}
  Therefore we can estimate that
  \begin{align*}
    &\int_{s_0}^{1} (r^{\ga-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga})^{-p} ds\\
    &\les \int_{s_0}^{1}\big((r_0-t_0)^{\ga}+(r_0-\tilde{r})^{\ga-2}(1-s)r_0^2\big)^{-p}ds\\
    &\les  (r_0-\tilde{r})^{2-\ga} r_0^{-2} (r_0-t_0)^{(1-p)\ga }.
  \end{align*}

   \begin{align*}
    &\int_{s_0}^{1} (r^{\ga-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga})^{-p} ds\\
    &\les \int_{s_0}^{1}\big((r_0-t_0+(r_0-\tilde{r})^{-1}r_0\tilde{r}(1-s))^{\ga}+(r_0-\tilde{r})^{\ga-2}(1-s)r_0^2\big)^{-p}ds\\
    &\les  (r_0-\tilde{r})^{2-\ga} r_0^{-2} (r_0-t_0)^{(1-p)\ga }.
  \end{align*}
  Let's consider for the particular case when $\tilde{r}=\f12 t_0$.
   \begin{align*}
    &\int_{s_0}^{1} (r^{\ga-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga})^{-p} ds\\
    &\les \int_{0}^{(1-r_0^{-1}\tilde{r})^2}\big((r_0-t_0+ r_0 s)^{\ga}+ r_0^{\ga} s\big)^{-p}ds.
  \end{align*}
  Let $s_*$ be the smallest solution of
  \begin{align*}
    (r_0-t_0+ r_0  s)^{\ga}= 4^{\ga}r_0^{\ga} s,\quad 0\leq s\leq (1-r_0^{-1}\tilde{r})^2.
  \end{align*}
  Then
    \begin{align*}
    &\int_{0}^{s_*} (r^{\ga-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga})^{-p} ds\\
    &\geq \int_{0}^{s_*} (r_0-t_0+ r_0 s)^{-p\ga} ds\\
    &\geq  u_0^{1-p\ga}(u_0+r_0s_*)^{-1}s_*=4^{-1}u_0^{1-p\ga}r_0^{-1}s_*^{1-\frac{1}{\ga}} .
  \end{align*}
  Notice that
  \begin{align*}
    r_0 s_*^{\frac{1}{\ga}}(4-s_*^{1-\frac{1}{\ga}})=u_0.
  \end{align*}
  Without loss of generality, one may assume that $t_0\geq \frac{9}{10}r_0$. Thus in particular, one has $r_0^{-1}u_0\leq \frac{1}{10}$. Therefore it can be checked that
  \begin{align*}
    s_*\leq 2 r_0^{-\ga}u_0^{\ga}\leq 2 \cdot 10^{-\ga}\leq \frac{1}{5}.
  \end{align*}
  Here we keep in mind that $\ga>1$.
  Hence from the above equation for $s_*$, we conclude that
  \begin{align*}
   2r_0^{-\ga}u_0^{\ga}\geq s_*\geq  10^{-1} r_0^{-\ga} u_0^{\ga}.
  \end{align*}
  Therefore
  \begin{align*}
    &\int_{0}^{s_*} (r^{\ga-1}(r+\tilde{r}-r_0s)+(r_0-t_0)^{\ga})^{-p} ds\\
    &\geq   u_0^{1-p\ga}r_0^{-1}s_*^{1-\frac{1}{\ga}}\geq u_0^{1-p\ga}r_0^{-1}(r_0^{-\ga}u_0^{\ga})^{1-\frac{1}{\ga}}=u_0^{\ga-p\ga}  r_0^{-\ga}.
    \end{align*}
    No improvement.
\end{proof}
It turns out that this lemma is implied by Lemma \ref{lem:integration:ex:ab} with $\a=p$, $\b=0$.

\fi

We are now ready to prove the following decay estimates for the solution in the exterior region.
\begin{Prop}
\label{prop:NLW:def:3d:ex}
In the exterior region $\{2+t\leq |x|\}$,
the solution $\phi$ to the equation \eqref{eq:NLW:semi:3d} satisfies the following $L^\infty$ decay estimates:
\begin{itemize}
\item when $p$ and $\ga_0$ verify the relation
\[
\frac{1+\sqrt{17}}{2}<p<5, \quad \max\{\frac{4}{p-1}-1, 1\}<\ga_0<\min\{p-1, 2\},
\]
then
\begin{equation}
\label{eq:phi:pt:Br:largep}
|\phi(t_0, x_0)|\leq C (1+t_0+|x_0|)^{-1}(1+|x_0|-t_0)^{-\frac{\ga_0-1}{2}}\sqrt{\mathcal{E}_{1, \ga_0}[\phi] };
\end{equation}
\item when $2<p\leq \frac{1+\sqrt{17}}{2}$ and $1<\ga_0<p-1$, then
\begin{equation}
\label{eq:phi:pt:Br:smallp}
|\phi(t_0, x_0)|\leq C  |x_0|^{-\frac{3+(p-2)^2}{(p+1)(5-p)}\ga_0}(|x_0|-t_0)^{-\frac{(p-1)\ga_0}{p+1}} \sqrt{\mathcal{E}_{1, \ga_0}[\phi] }
\end{equation}
for some constant $C$ depending on $\ga_0$, $p$ and the zeroth order weighted energy $\mathcal{E}_{0, \ga_0}[\phi] $.
\end{itemize}
\end{Prop}
\begin{proof}
The proof for this decay estimate relies on the representation formula for linear wave equations. The nonlinearity will be controlled by using the weighted energy estimates in Proposition \ref{prop:EF:cone:NW:3d}. Note that for $q=(t_0, x_0)$ in the exterior region, we have
\begin{equation}
\label{eq:rep4phi:ex}
\begin{split}
4\pi\phi(t_0, x_0)&=\int_{\tilde{\om}}t_0  \phi_1(x_0+t_0\tilde{\om})d\tilde{\om}+\pa_{t_0}\big(\int_{\tilde{\om}}t_0  \phi_0(x_0+t_0\tilde{\om})d\tilde{\om}   \big)-\int_{\mathcal{N}^{-}(q)}|\phi|^{p-1} \phi \ \tilde{r} d\tilde{r}d\tilde{\om}.
\end{split}
\end{equation}
For the linear evolution part, one can use the standard vector field method to show that
\begin{align*}
  |\int_{\tilde{\om}}t_0  \phi_1(x_0+t_0\tilde{\om})d\tilde{\om}+\pa_{t_0}\big(\int_{\tilde{\om}}t_0  \phi_0(x_0+t_0\tilde{\om})d\tilde{\om}   \big)|

  &\les r_0^{-1} (r_0-t_0)^{-\frac{\ga_0-1}{2}}\sqrt{\mathcal{E}_{1, \ga_0}[\phi] }
\end{align*}
for $\ga_0>1$ and $r_0=|x_0|\geq t_0+2$.

For the case when $ \frac{1+\sqrt{17}}{2}<p<5$, by using the weighted energy estimate \eqref{eq:Eflux:ex:EF} and the bound \eqref{eq:integration:ex:ab} with $\a=\frac{p-1}{2}$, $\b=\frac{p+1}{2}$, we can estimate that
\begin{align*}
  &|\int_{\mathcal{N}^{-}(q)}\Box \phi \ \tilde{r} d\tilde{r}d\tilde{\om}|\\
  &\leq \left(\int_{\mathcal{N}^{-}(q)}((1+\tau)r^{\ga_0}+(r_0-t_0)^{\ga_0})|\phi|^{p+1}\ \tilde{r}^{2} d\tilde{r}d\tilde{\om}\right)^{\frac{p-1}{p+1}}\\
  &\quad \cdot \left(\int_{\mathcal{N}^{-}(q)}((1+\tau)r^{\ga_0}+(r_0-t_0)^{\ga_0})^{-\frac{p-1}{2}}|\phi|^{\frac{p+1}{2}}\ \tilde{r}^{\frac{3-p}{2}} d\tilde{r}d\tilde{\om}\right)^{\frac{2}{p+1}}\\

  &\les \left(\int_{0}^{t_0}(\sup\limits_{x}|r\phi|^{\frac{p+1}{2}}) (r_0-\tilde{r})^{\frac{3-p}{2}-\ga_0+\ep} r_0^{-2}\left((r_0-\tilde{r})^{\frac{3-p}{2}\ga_0}+(r_0-t_0)^{\frac{3-p}{2}\ga_0}\right)\tilde{r}^{\frac{3-p}{2}} d\tilde{r} \right)^{\frac{2}{p+1}}.
\end{align*}
When $p\geq 3$, we estimate that
\begin{align*}
   (r_0-\tilde{r})^{\frac{3-p}{2}-\ga_0+\ep} r_0^{\frac{p-3}{2}}\left((r_0-\tilde{r})^{\frac{3-p}{2}\ga_0}+(r_0-t_0)^{\frac{3-p}{2}\ga_0}\right)\leq  (2+t_0-\tilde{r})^{\frac{3-p}{2}-\ga_0+\ep} (2+t_0)^{\frac{p-3}{2}}.
\end{align*}
When $p<3$, we can choose $\ep$ sufficiently small such that $$\frac{3-p}{2}-\ga_0+\ep+\frac{3-p}{2}\ga_0\leq 0$$ due to the assumption $(p-1)(\ga_0+1)>4$ for this case. Thus we can bound that
\begin{align*}
   (r_0-\tilde{r})^{\frac{3-p}{2}-\ga_0+\ep} r_0^{\frac{p-3}{2}}\left((r_0-\tilde{r})^{\frac{3-p}{2}\ga_0}+(r_0-t_0)^{\frac{3-p}{2}\ga_0}\right)\leq  (2+t_0-\tilde{r})^{\frac{3-p}{2}+\frac{1-p}{2}\ga_0+\ep} (2+t_0)^{\frac{p-3}{2}}.
\end{align*}
We therefore derive that
\begin{align*}
  |\phi(t_0, x_0)|&\les \sqrt{\mathcal{E}_{1, \ga_0}[\phi](\B_2^{\infty})}r_0^{-1}(r_0-t_0)^{\frac{1-\ga_0}{2}}\\
  &+\left(r_0^{-\frac{p+1}{2}}\int_{0}^{t_0}(\sup\limits_{x}|r\phi|^{\frac{p+1}{2}})(2+t_0-\tilde{r})^{\frac{3-p}{2}-\ga_0+\ep+\frac{\ga_0}{2}\max\{0, 3-p\}} (2+t_0)^{\frac{p-3}{2}} \tilde{r}^{\frac{3-p}{2}} d\tilde{r} \right)^{\frac{2}{p+1}}.
\end{align*}
Now define the function
\begin{equation*}
\mathcal{M}(t_0)=\sup\limits_{|x|\geq t_0+2} |u_+^{\frac{\ga_0-1}{2}}r\phi|^{\frac{p+1}{2}}
\end{equation*}
and
\begin{equation*}
  f(t_0, \tilde{r})=(2+t_0-\tilde{r})^{\frac{3-p}{2}-\ga_0+\ep+\frac{\ga_0}{2}\max\{0, 3-p\}} (2+t_0)^{\frac{p-3}{2}} \tilde{r}^{\frac{3-p}{2}}.
\end{equation*}
Here $u_+=1+\f12|t-r|$. Since on the cone $\mathcal{N}^{-}(q)$, $u_+\geq \f12 (r_0-t_0)$,
then the previous inequality leads to
\begin{align*}
  \mathcal{M}(t_0)\les (\mathcal{E}_{1, \ga_0}[\phi] )^{\frac{p+1}{4}}+\int_{0}^{t_0}\mathcal{M}(t_0-\tilde{r})f(t_0, \tilde{r}) d\tilde{r}.
\end{align*}
When $p\geq 3$, by the assumption on $p$ and $\ga_0$ of the main Theorem \ref{thm:main}, we in particular have $\ga_0>1$. Choose $\ep$ sufficiently small such that $0<\ep<\ga_0-1$. Then since $3\leq p<5$, we can bound that
\begin{align*}
  \int_0^{t_0}f(t_0,\tilde{r})d\tilde{r} &=\int_0^{t_0}(2+t_0-\tilde{r})^{\frac{3-p}{2}-\ga_0+\ep} (2+t_0)^{\frac{p-3}{2}} \tilde{r}^{\frac{3-p}{2}}d\tilde{r}\\
  &\les \int_0^{\f12 t_0}(2+t_0)^{-\ga_0+\ep}  \tilde{r}^{\frac{3-p}{2}}d\tilde{r}+\int_{\f12 t_0}^{t_0}(2+t_0-\tilde{r})^{\frac{3-p}{2}-\ga_0+\ep} d\tilde{r}\\
  &\les (2+t_0)^{-\ga_0+\ep}t_0^{\frac{5-p}{2}}+2^{\frac{5-p}{2}-\ga_0+\ep}+(2+t_0)^{\frac{5-p}{2}-\ga_0+\ep}\\
  &\les 1
\end{align*}
as $\frac{5-p}{2}+\ep<\ga_0$. The implicit constant relies only on $p$, $\ga_0$ and $\ep$. For the case when $p<3$, the assumption implies that $(\ga_0+1)(p-1)>4$. Let $\ep>0$ verify the relation $$\frac{5-p}{2}-\frac{(p-1)\ga_0}{2}+\ep<0.$$ Then
\begin{align*}
  \int_0^{t_0}f(t_0,\tilde{r})d\tilde{r} &=\int_0^{t_0}(2+t_0-\tilde{r})^{\frac{3-p}{2}-\ga_0+\frac{(3-p)\ga_0}{2}+\ep} (2+t_0)^{\frac{p-3}{2}} \tilde{r}^{\frac{3-p}{2}}d\tilde{r}\\
  &\les \int_0^{\f12 t_0}(2+t_0)^{-\frac{\ga_0(p-1)}{2}+\ep}  \tilde{r}^{\frac{3-p}{2}}d\tilde{r}+\int_{\f12 t_0}^{t_0}(2+t_0-\tilde{r})^{\frac{3-p}{2}-\frac{p-1}{2}\ga_0+\ep} d\tilde{r}\\
  &\les (2+t_0)^{-\frac{p-1}{2}\ga_0+\ep}t_0^{\frac{5-p}{2}}+2^{\frac{5-p}{2}-\frac{p-1}{2}\ga_0+\ep}+(2+t_0)^{\frac{5-p}{2}-\frac{p-1}{2}\ga_0+\ep}\\
  &\les 1.
\end{align*}
In any case the function $f(t_0, \tilde{r})$ is integrable on the interval $[0, t_0]$, which is independent of $t_0$. Thus by using Gronwall's inequality, we conclude that
\begin{align*}
   \mathcal{M}(t_0)\les (\mathcal{E}_{1, \ga_0}[\phi] )^{\frac{p+1}{4}}.
\end{align*}
The decay estimate for $\phi$ then follows from the definition of $\mathcal{M}(t_0)$ for the case when $p>\frac{1+\sqrt{17}}{2}$.

\bigskip

For the smaller power $p$, we instead have weaker decay estimates for the solution. Since in this case we assumed $1<\ga_0<p-1\leq \frac{\sqrt{17}-1}{2}$, in particular we have $$p\ga_0>2, \quad (p-1)(\ga_0+1)>2,\quad \frac{(p-1)\ga_0}{5-p}<\frac{(p-1)^2}{5-p}<1.$$ To estimate the nonlinear term, we split the integral on the backward light cone $\mathcal{N}^{-}(q)$ into two parts: the first part is restricted to time interval $[\f12 t_0^{\frac{(p-1)\ga_0}{5-p}}, t_0]$, on which we use the same argument as above and the second part can be directly estimated by the weighted energy flux. More precisely we have
\begin{align*}
  &|\int_{\mathcal{N}^{-}(q)\cap\{\tilde{r}\geq \f12 t_0^{\frac{p-1}{5-p}\ga_0}\}}\Box \phi \ \tilde{r} d\tilde{r}d\tilde{\om}|\\
  &\les \left(\int_{\mathcal{N}^{-}(q)}((1+\tau)r^{\ga_0}+(r_0-t_0)^{\ga_0})|\phi|^{p+1}\ \tilde{r}^{2} d\tilde{r}d\tilde{\om}\right)^{\frac{p}{p+1}}\\
  &\quad \cdot \left(\int_{\mathcal{N}^{-}(q)\cap\{\tilde{r}\geq \f12 t_0^{\frac{p-1}{5-p}\ga_0}\}}((1+\tau)r^{\ga_0}+(r_0-t_0)^{\ga_0})^{-p} \tilde{r}^{1-p} d\tilde{r}d\tilde{\om}\right)^{\frac{1}{p+1}}\\
  &\les  \left(\int_{\f12 t_0 ^{\frac{p-1}{5-p}\ga_0}}^{t_0}  (r_0-\tilde{r})^{2-\ga_0}   r_0^{-2}(r_0-t_0)^{(1-p)\ga_0} \tilde{r}^{1-p} d\tilde{r} \right)^{\frac{1}{p+1}}\\
  &\les r_0^{-\frac{3+(p-2)^2}{(p+1)(5-p)}\ga_0}(r_0-t_0)^{\frac{(1-p)\ga_0}{p+1}}.
\end{align*}
Here we used estimate \eqref{eq:integration:ex:ab} without loss of $\ep$ as $p<3$ (see the proof for Lemma \ref{lem:integration:ex:ab}). Thus we derive that
\begin{align*}
  |\int_{\mathcal{N}^{-}(q) }\Box \phi \ \tilde{r} d\tilde{r}d\tilde{\om}|
  \les  &r_0^{-\frac{3+(p-2)^2}{(p+1)(5-p)}\ga_0}(r_0-t_0)^{\frac{(1-p)\ga_0}{p+1}}\\
 &+ \left(\int_{0}^{ \f12 t_0^{\frac{(p-1)\ga_0}{5-p}}}(\sup\limits_{x}|r\phi|^{\frac{p+1}{2}}) (r_0-\tilde{r})^{1-\frac{(p-1)(\ga_0+1)}{2}} r_0^{-2} \tilde{r}^{\frac{3-p}{2}} d\tilde{r} \right)^{\frac{2}{p+1}}.

\end{align*}
Again here we do not loss the $\ep $ decay of the bound \eqref{eq:integration:ex:ab} as $p\leq \frac{1+\sqrt{17}}{2}<3$. Since $\tilde{r}\leq \f12 t_0^{\frac{(p-1)\ga_0}{5-p}}\leq \f12 t_0$, we in particular have that $r_0\les r$, $r_0\les r_0-\tilde{r}$. Define
\begin{align*}
  \mathcal{M}_1(t)=\sup\limits_{|x|\geq t+2}( |\phi(t, x)| |x|^{\frac{3+(p-2)^2}{(p+1)(5-p)}\ga_0 }(|x|-t)^{\frac{(p-1)\ga_0}{p+1}})^{\frac{p+1}{2}}.
\end{align*}
We then derive that
\begin{align*}
  \mathcal{M}_1(t_0)&\les 1+\int_0^{\f12 t_0^{\frac{(p-1)\ga_0}{5-p}}}\mathcal{M}_1(t_0-\tilde{r}) r_0^{\frac{p+1}{2}+1-\frac{(p-1)(\ga_0+1)}{2}-2}\tilde{r}^{\frac{3-p}{2}}d\tilde{r}\\
  & \les  1+\int_0^{\f12 t_0^{\frac{(p-1)\ga_0}{5-p}}}\mathcal{M}_1(t_0-\tilde{r}) t_0^{ -\frac{(p-1)\ga_0 }{2} } \tilde{r}^{\frac{3-p}{2}}d\tilde{r}.
\end{align*}
Here we note that
\begin{align*}
  \frac{3+(p-2)^2}{(p+1)(5-p)}\ga_0\leq 1,\quad \frac{3+(p-2)^2}{(p+1)(5-p)}\ga_0+\frac{(p-1)\ga_0}{p+1}\leq \frac{\ga_0+1}{2}.
\end{align*}
By using Gronwall's inequality, we conclude that
\begin{align*}
  \mathcal{M}_1(t_0)\les 1.
\end{align*}
The pointwise decay estimate for $\phi$ for the case when $2<p\leq \frac{1+\sqrt{17}}{2}$ then follows.
\end{proof}

In order to study the asymptotic behaviour of the solution in the interior region $\{t+2\geq |x|\}$, we use the method of conformal compactification, which requires to understand the solution on the hyperboloid $\Hy$ defined in the Section \ref{sec:notation}.
\begin{Prop}
\label{prop:energyflux:H:ex}
Assume that $p$ and $\ga_0$ verifies the relation
\begin{align*}
  2<p<5,\quad 1<\ga_0<\min\{2, p-1\}.
\end{align*}
 Then we have the following weighted energy flux bound through the future part of the hyperboloid
\begin{equation}
\label{eq:EE:hyB}
\begin{split}
E[\phi](\Hy)+\int_{\mathbb{H}^+}r^{\ga_0}|L(r\phi)|^2+|{\Lb} \phi|^2+r^2| L \phi|^2+|\nabb (r\phi)|^2+\frac{2r^2}{p+1}|\phi|^{p+1} dtd\om\leq C \mathcal{E}_{0, \ga_0}[\phi]
\end{split}
\end{equation}
for some universal constant $C$.

In addition, for the large $p$ case
\[
\frac{1+\sqrt{17}}{2}<p<5, \quad \max\{\frac{4}{p-1}-1, 1\}<\ga_0<\min\{p-1, 2\},
\]
we also have the energy bound for the first order derivatives
\begin{equation}
\label{eq:EE:hyB:1derivative}
\begin{split}
E[Z\phi](\Hy)+\int_{\mathbb{H}^+}r^{\ga_0}|L Z(r\phi)|^2+|{\Lb} Z \phi|^2+r^2| L Z\phi|^2+|\nabb Z (r\phi)|^2  dtd\om\leq C \mathcal{E}_{1, \ga_0}[\phi]^{p-1}
\end{split}
\end{equation}
for all $Z\in \Gamma=\{\pa_{\mu}, \Om_{\mu\nu}=x^{\mu}\pa_{\nu}-x^{\nu}\pa_\mu\}$ and some constant $C$ depending on $\mathcal{E}_{0, \ga_0}[\phi]$, $p$ and $\ga_0$. Here the hyperboloid $\Hy$ is parameterized by $(t, \om)$ and $E[\phi](\Hy)$ denotes the energy flux of $\phi$ through $\Hy$.
\end{Prop}
\begin{proof}
 The proof goes in the same manner by applying the energy identity \eqref{eq:energy:id} to the same vector fields $X$, $Y$ and function $\chi$ in the proof of the main theorem in the author's companion paper \cite{yang:scattering:NLW} but for the domain $\mathcal{D}$ being the subregion of the exterior region bounded by $\Hy^+$, the initial hypersurface and the incoming null hypersurface $\Hb_{v_0}$. The bulk integral and the boundary integral on the initial hypersurface as well as the incoming null hypersurface $\Hb_{v_0}$ could be found in section 4 in \cite{yang:scattering:NLW}.
 It remains to compute the boundary integral on the hyperboloid $\Hy^+$.

Define the functions
\[
\tau_1=\f12\sqrt{(t^*-R^*)^2+r^2},\quad \tau_0=\frac{t^*-R^*}{r}=\sqrt{1+r^{-2}(R^*)^2}.
\]
Here recall that$ R^*=\frac{5}{6}$, $t^*=t+3$. In particular we have
\[
d\tau_1=\f12\tau_1^{-1}((t^*-R^*)dt+rdr),\quad \pa_{\tau_1}=\frac{\tau_1}{t^*-R^*}\pa_t+\frac{\tau_1}{r}\pa_r.
\]
 Then the hyperboloid $\mathbb{H}^+$ can be parameterized by $(\tau_1, \om)$ or $(t, \om)$.
 We therefore can compute that
\begin{align*}
-2\int_{\mathbb{H}^+\cap\{v\leq v_0\}}i_{J^{X, Y, \chi}[\phi]}d\vol&=2\int_{\mathbb{H}^+ \cap\{v\leq v_0\}} (J^{X, Y, \chi}[\phi])^u (dt+dr)r^2d\om+(J^{X, Y, \chi}[\phi])^v (dr-dt) r^2 d\om\\
&=\int_{\mathbb{H}_+ \cap\{v\leq v_0\} } (1+\tau_0)r^{\ga_0}|L(r\phi)|^2dtd\om - \int_{\mathbb{H}^+ \cap\{v\leq v_0\}}\pa_\tau (r^{\ga_1+1}|\phi|^2)d\tau d\om \\
&\quad +\int_{\mathbb{H}^+ \cap\{v\leq v_0\} }(\tau_0-1)r^{\ga_0}(|\nabb(r\phi)|+\frac{2r^2}{p+1}|\phi|^{p+1})\quad  dt d\om.
\end{align*}
Here for the particular choice of $X=r^{\ga_0} L$, $Y=\f12 \ga_0 r^{\ga_0-2}|\phi|^2 L$ and $\chi=r^{\ga_0-1}$, we can compute that
\begin{align*}
-2r^2(J^{X, Y, \chi}[\phi])^u&= r^{\ga_0}|L(r\phi)|^2-\f12 L(r^{\ga_0+1}|\phi|^2),\\
2r^2(J^{X, Y, \chi}[\phi])^v&= -r^{\ga_0}(|\nabla(r\phi)|^2+\frac{2r^2}{p+1}|\phi|^{p+1})-\f12 \Lb(r^{\ga_0+1}|\phi|^2).
\end{align*}
By using integration by parts, we have the identity
\begin{align*}
  -\int_{\mathbb{H}^+ \cap\{v\leq v_0\}}\pa_\tau (r^{\ga_1+1}|\phi|^2)d\tau d\om+\int_{\B_2^{\infty}\cap\{v\leq v_0\}} \pa_r(r^{\ga_0+1}|\phi|^2)r^{-2}dx+\int_{\Hb_{v_0}}\Lb(r^{\ga_0+1}|\phi|^2)du d\om=0.
\end{align*}
We therefore conclude from the energy identity \eqref{eq:energy:id} that
\begin{align*}
  \int_{\Hy^+}r^{\ga_0}|L(r\phi)|^2 dtd\om\leq \int_{\B_{2}^{\infty}}  r^{\ga_0}(  r^{-2}|L(r\phi)|^2 +|\nabb\phi|^2+\frac{2}{p+1}|\phi|^{p+1}) dx\leq C \mathcal{E}_{0, \ga_0}[\phi]
\end{align*}
for some universal constant $C$. For more details, we refer the interested reader to \cite{yang:scattering:NLW}.

Now to prove the estimate \eqref{eq:EE:hyB}, we conduct the classical energy estimate derived by using the vector field $\pa_t$ as multiplier. It suffices to compute the energy flux through the hyperboloid $\Hy$ for solution $\phi$ of \eqref{eq:NLW:semi:3d}. For this we compute that
\begin{align*}
E[\phi](\Hy)&=-2\int_{\Hy}i_{J^{\pa_t, 0, 0}[\phi]}d\vol\\
&=-2\int_{\mathbb{H} } (J^{\pa_t, 0, 0}[\phi])^0 r^2 dr d\om-(J^{\pa_t, 0, 0}[\phi])^r r^2 dt  d\om\\
 &=\int_{\mathbb{H} }(\tau_0 (|\pa\phi|^2+\frac{2}{p+1}|\phi|^{p+1} ) +2\pa_t\phi \pa_r\phi ) \quad r^2 dt  d\om\\
 &=\int_{\mathbb{H} }(\tau_0 (|\nabb\phi|^2+\frac{2}{p+1}|\phi|^{p+1} )+\frac{\tau_0-1}{2}|\Lb\phi|^2+\frac{1+\tau_0}{2}|L\phi|^2 ) \quad r^2 dt  d\om.
\end{align*}
Since $1\leq \tau_0\leq 2$ on $\Hy^+$ and $\tau_0-1=(R^*)^2(1+\tau_0)^{-1}r^{-2}$, the energy conservation then leads
\begin{align*}
  \int_{\Hy^+} (|L\phi|^2+|\nabb\phi|^2+\frac{2}{p+1}|\phi|^{p+1}+r^{-2}|\Lb\phi|^2)r^2 dtd\om\leq C \mathcal{E}_{0, 0}[\phi]
\end{align*}
for some universal constant $C$.
This together with the above weighted energy bound for $|L(r\phi)|$ implies the inequality \eqref{eq:EE:hyB}.

As for the energy estimate \eqref{eq:EE:hyB:1derivative} for the derivatives,
consider the equation for $Z\phi$
\[
\Box Z\phi=Z(|\phi|^{p}\phi)
\]
with nonlinearity $Z(|\phi|^p\phi)$. The associated energy momentum tensor for $Z\phi$ is
\begin{align*}
  T[Z\phi]_{\mu\nu}=\pa_{\mu}Z\phi \pa_\nu Z\phi-\f12 m_{\mu\nu}\pa^\ga Z\phi \pa_\ga Z\phi.
\end{align*}
The energy identity \eqref{eq:energy:id} still holds but without the potential part $|\phi|^{p}\phi$. The above computations for $\phi$ then lead to
\begin{equation}
\label{eq:bd4LZrphi}
\begin{split}
  &\int_{\mathbb{H}^+}r^{\ga_0}|L Z(r\phi)|^2+|{\Lb} Z^k \phi|^2+r^2| L Z^k\phi|^2+|\nabb Z^k (r\phi)|^2 dtd\om\\
  &\les \mathcal{E}_{1, \ga_0}[\phi] +\iint_{\mathcal{D}}|X (Z\phi)+\chi Z\phi \Box Z\phi|+ |\Box Z\phi||\pa_t Z\phi|dxdt\\
  &\les \mathcal{E}_{1, \ga_0}[\phi] +\iint_{\mathcal{D}}(r^{\ga_0-1}|L(r Z\phi)|+|\pa_t Z\phi|)|Z\phi| |\phi|^{p-1} dx dt.
  \end{split}
\end{equation}
Here recall that the region $\mathcal{D}$ is bounded by the hyperboloid $\Hy^+$ and the initial hypersurface and by our convention, the implicit constant relies only on $\mathcal{E}_{0, \ga_0}[\phi]$, $\ga_0$ and $p$.

To bound the bulk integral on the right hand side of the above inequality, we instead apply the above $r$-weighted energy estimate and the classical energy estimate to the domain $\cD_{u_1}^{u_2}$ bounded by the outgoing null hypersurface $\H_{u_1}$, the incoming null hypersurface $\Hb_{-u_2}$ and the initial hypersurface. The $r$-weighted energy estimate with the same choice of the vector fields $X$, $Y$ and the function $\chi$ shows that
\begin{align*}
  \int_{\H_{u_1}}r^{\ga_0} |L(rZ\phi)|^2 dvd\om &\les \mathcal{E}_{1, \ga_0}[\phi]  +\iint_{\cD_{u_1}^{u_2}}r^{\ga_0-1}|L(rZ\phi)||Z\phi||\phi|^{p-1} dxdt\\
  &\les \mathcal{E}_{1, \ga_0}[\phi]  +\iint_{\cD_{u_1}^{u_2}} r^{\ga_0-2}|L(rZ\phi)|^2u_+^{-1-\ep}+ u_+^{1+\ep} r^{\ga_0}|Z\phi|^2|\phi|^{2p-2} dxdt.
\end{align*}
The integral of the first term on the right hand side can be absorbed by using Gronwall's inequality. For the integral of the nonlinearity, we rely on the pointwise decay estimate of $\phi$ obtained in Proposition \ref{prop:NLW:def:3d:ex} as well as the energy estimate for $\phi$ in Proposition \ref{prop:EF:cone:NW:3d}. First we note that in the exterior region
\begin{align*}
  |Z\phi|^2\les |L(r\phi)|^2+ u_+^2|\Lb \phi|^2+ r^2 |\nabb\phi|^2+|\phi|^2,\quad \forall t+2\leq |x|.
\end{align*}
The integral of $|L(r\phi)|^2+r^2|\nabb\phi|^2$ can be bounded by using the weighted energy estimate \eqref{eq:Eflux:ex:EF}. For $|\Lb\phi|^2$, recall the energy estimate for $\phi$
\begin{align*}
  \int_{\Hb_{-u_2}^{u_1}}|\Lb\phi|^2 \les (u_1)_+^{-\ga_0}\mathcal{E}_{0, \ga_0}[\phi] ,\quad \forall u_2<u_1\leq -1.
\end{align*}
To bound the integral of $|\phi|^2$, we rely on the $r$-weighted energy estimate for $\phi$ through the outgoing null hypersurface $\H_u$
\begin{align*}
  \int_{\H_u}r^{-1-\ep}|\phi|^2\les \int_{\om}|r\phi(-u, u, \om)|^2 d\om+u_+^{1-\ga_0}\int_{\H_u}r^{\ga_0}|L(r\phi)|^2 dvd \om \les u_+^{1-\ga_0}\mathcal{E}_{0, \ga_0}[\phi],\quad \forall u\leq -1.
\end{align*}
By our assumption on $p$, $\ga_0$, we in particular have the lower bound for $p$ and choose $\ep$ such that
\[
p>\frac{1+\sqrt{17}}{2}>\frac{5}{2}, \quad (\ga_0+1)(p-1)>4+3\ep.
\]
Therefore
 we can show that
\begin{align*}
  &\iint_{\cD_{u_1}^{u_2}}   u_+^{1+\ep} r^{\ga_0}|Z\phi|^2|\phi|^{2p-2} dxdt\\
  &\les \mathcal{E}_{1, \ga_0}[\phi]^{p-1} \iint_{\cD_{u_1}^{u_2}}   u_+^{1+\ep-(\ga_0-1)(p-1)} r^{\ga_0-2p+2}(|L(r\phi)|^2+ u_+^2|\Lb \phi|^2+ r^2 |\nabb\phi|^2+|\phi|^2)\\

  &\les \mathcal{E}_{1, \ga_0}[\phi]^{p-1} \iint_{\cD_{u_1}^{u_2}}   u_+^{\ga_0-2-\ep} r^{-1-\ep}(u_+^2|\Lb \phi|^2+|\phi|^2)+r^{\ga_0-3}(|L(r\phi)|^2 + r^2 |\nabb\phi|^2 )\\
  &\les \mathcal{E}_{1, \ga_0}[\phi]^{p-1} \mathcal{E}_{0, \ga_0}[\phi].
\end{align*}
This in particular implies that
\begin{align*}
  \int_{\H_{u_1}}r^{\ga_0} |L(rZ\phi)|^2 dvd\om  +\iint_{\cD_{u_1}^{u_2}}r^{\ga_0-1}|L(rZ\phi)||Z\phi||\phi|^{p-1} dxdt \les \mathcal{E}_{1, \ga_0}[\phi]^{p-1}.
\end{align*}
Here without loss of generality we may assume that $\mathcal{E}_{1, \ga_0}[\phi]\geq 1.$ Based on these computations and by using energy estimate for $Z\phi$, we further can show that
\begin{align*}
 \int_{\H_{u_1}^{u_2}}|L Z\phi|^2+\int_{\Hb_{-u_2}^{u_1}}|\Lb Z\phi|^2 &\les \mathcal{E}_{1, 0}[\phi]  +\iint_{\cD_{u_1}^{u_2}}|\pa_t Z\phi||Z\phi||\phi|^{p-1} dxdt\\
 &\les (u_1)_+^{-\ga_0}\mathcal{E}_{1, \ga_0}[\phi]  +\iint_{\cD_{u_1}^{u_2}}r_+^{-1-\ep}|\pa_t Z\phi|^2+r^{1+\ep}|Z\phi|^2|\phi|^{2p-2} dxdt.
\end{align*}
The integral of the first term can be absorbed by using Gronwall's inequality while the second term has been estimated above by choosing $\ep$ such that $1+\ep<\ga_0$. We hence conclude that
\begin{align*}
  \iint_{\cD_{u_1}^{u_2}}|\pa_t Z\phi||Z\phi||\phi|^{p-1} \les (u_1)_+^{-\ga_0}\mathcal{E}_{1, \ga_0}[\phi]+(u_1)_+^{-\ga_0}\mathcal{E}_{1, \ga_0}[\phi]^{p-1}\les (u_1)_+^{-\ga_0}\mathcal{E}_{1, \ga_0}[\phi]^{p-1}.
\end{align*}
The weighted energy estimate \eqref{eq:EE:hyB:1derivative} then follows in view of \eqref{eq:bd4LZrphi}.
\end{proof}

The above proposition will play the role that the solution in the interior region has uniform bounded energy flux. The method for proving the decay estimates is similar to the above argument for deriving the decay estimates for the solution in the exterior region after conformal transformation. The nonlinearity will be controlled by the weighted energy flux through backward light cone. To use Gronwall's inequality, one needs first bound the linear evolution with prescribed data, for which, in the interior region, will be the data on the hyperboloid $\Hy$. For large $p$ when the solution decays sufficiently fast, one can use the standard energy estimates to control the linear evolution, which however, fails for the smaller $p$ case. We instead rely on the representation formula together with the uniform weighted energy flux bound through backward light cones.

Define the region inclosed by the hyperboloid $\Hy$
\begin{align*}
\mathbf{D}:=\left\{(t, x)|(t^*)^2-|x|^2\geq (R^*)^{-1} t^*\right\},\quad \mathbf{D}^+=\mathbf{D}\cap\{t\geq 0\}, \quad R^*=\frac{5}{6},\quad t^*=t+3.
\end{align*}
Let $\phi^{lin}_{H}$ be the linear evolution in $\mathbf{D}$, that is,
\begin{align*}
  \Box \phi^{lin}_{H}=|\phi|^{p-1}\phi (1-\mathbf{1}_{\mathbf{D}^+}),\quad \phi^{lin}_H(0, x)=\phi_0,\quad \pa_t \phi_H^{lin}(0, x)=\phi_1,
\end{align*}
where $\mathbf{1}_{\mathbf{D}^+}$ stands for the characteristic function of the set $\mathbf{D}^+$. We see that $\phi^{lin}_H$ coincides with $\phi$ in the region $(\mathbb{R}^{3}/\mathbf{D})\cap\{t\geq 0\}$.

We have the following estimate for $\phi^{lin}_H$ inside $\mathbf{D}$.
\begin{Prop}
\label{prop:NLW:def:3d:ex:ID:H}
Let $p$ and $\ga_0$ verify the same assumptions as in Proposition \ref{prop:NLW:def:3d:ex}. Then inside the hyperboloid $\mathbf{D}$, for large $p>\frac{1+\sqrt{17}}{2}$, we have
\begin{equation}
\label{eq:phi:pt:Br:largep:lin}
|\phi^{lin}_H(t_0, x_0)|\leq C (2+t_0+|x_0|)^{-1}(2+||x_0|-t_0|)^{-\frac{\ga_0-1}{2}}\mathcal{E}_{1, \ga_0}[\phi]^{\frac{p-1}{2}},
\end{equation}
while for the case $2<p\leq \frac{1+\sqrt{17}}{2}$ and $1<\ga_0<p-1$, we have
\begin{equation}
\label{eq:phi:pt:Br:smallp:lin}
|\phi^{lin}_H (t_0, x_0)|\leq C  (2+t_0+|x_0|)^{-\frac{3+(p-2)^2}{(p+1)(5-p)}\ga_0}(1+||x_0|-t_0|)^{-\frac{\ga_0}{p+1}} \sqrt{\mathcal{E}_{1, \ga_0}[\phi] }
\end{equation}
 for some constant $C$ depending on $\ga_0$, $p$ and the zeroth order weighted energy $\mathcal{E}_{0, \ga_0}[\phi] $.
\end{Prop}
\begin{proof}
The larger $p$ case of estimate \eqref{eq:phi:pt:Br:largep:lin} follows directly by using the standard energy method, in view of the weighted energy bounds \eqref{eq:EE:hyB}, \eqref{eq:EE:hyB:1derivative} from the previous Proposition for the initial data for $\phi^{lin}_H$ on $\Hy$. Details for this decay estimate for linear waves could be found, for example, in \cite{yang1}.

 For the small $2<p\leq \frac{1+\sqrt{17}}{2}$ case, which requires that $\ga_0<p-1$,  the above energy method fails. Denote
\[
u_0=1+|t_0-|x_0||,\quad v_0=2+t_0+|x_0|.
\]
 Recall that for $q=(t_0, x_0)$, we have
\begin{equation*}
\begin{split}
4\pi\phi^{lin}_H(t_0, x_0)&=\int_{\tilde{\om}}t_0  \phi_1(x_0+t_0\tilde{\om})d\tilde{\om}+\pa_{t_0}\big(\int_{\tilde{\om}}t_0  \phi_0(x_0+t_0\tilde{\om})d\tilde{\om}   \big)-\int_{\mathcal{N}^{-}(q)/\mathbf{D}}|\phi|^{p-1} \phi \ \tilde{r} d\tilde{r}d\tilde{\om}.
\end{split}
\end{equation*}
Decay estimates for the linear evolution part can be carried out by using standard vector field method
\begin{align*}
  |\int_{\tilde{\om}}t_0  \phi_1(x_0+t_0\tilde{\om})d\tilde{\om}|+|\pa_{t_0}\big(\int_{\tilde{\om}}t_0  \phi_0(x_0+t_0\tilde{\om})d\tilde{\om}   \big)|\leq C v_0^{-1}u_0^{-\frac{\ga_0-1}{2}}\sqrt{\mathcal{E}_{1, \ga_0}[\phi]}.
\end{align*}
We now need to control the contribution of the nonlinear part from the exterior region. The case when $t_0\leq 20$ is trivial since in this case $|t_0|+|x_0|\leq 20$ (confined in $\mathbf{D}$). Minor modification of the argument for estimating the nonlinear terms in Proposition \ref{prop:NLW:def:3d:ex} also applies to the case $|t_0-|x_0||\leq 2$ (that is Lemma \ref{lem:integration:ex:ab} holds for $|t_0-|x_0||\leq 10$). Alternatively by moving the origin around, the decay estimates of Proposition \ref{prop:NLW:def:3d:ex} are also valid for $q=(t_0, x_0)$ with $|t_0-|x_0||\leq 10$. Hence in the sequel, it suffices to consider the case when $t_0\geq 20$ and $t_0>|x_0|+10$.

First we can estimate that
\begin{align*}
  |\int_{\mathcal{N}^{-}(q)/\mathbf{D}}|\phi|^{p-1} \phi \ \tilde{r} d\tilde{r}d\tilde{\om}|&\leq \left(\int_{\mathcal{N}^{-}(q)}|\phi|^{p+1}((1+\tau)v_+^{\ga}+u_+^{\ga})\tilde{r}^2 d\tilde{r}d\tilde{\om}\right)^{\frac{p}{p+1}}\\
  & \quad \cdot \left(\int_{\mathcal{N}^{-}(q)/\mathbf{D}}((1+\tau)v_+^{\ga}+u_+^{\ga})^{-p}\tilde{r}^{1-p} d\tilde{r}d\tilde{\om}\right)^{\frac{1}{p+1}}\\
  &\leq (\mathcal{E}_{0, \ga_0}[\phi])^{\frac{p}{p+1}} \left(\int_{\mathcal{N}^{-}(q)/\mathbf{D}}((1+\tau)v_+^{\ga}+u_+^{\ga})^{-p}\tilde{r}^{1-p} d\tilde{r}d\tilde{\om}\right)^{\frac{1}{p+1}}
\end{align*}
for all $1<\ga<\ga_0$.
Denote $s=-\om_0\cdot \tilde{\om}$ with $\om_0=r_0^{-1}x_0$. Recall that
  \begin{align*}
    &r^2=|x|^2=(\tilde{r}-r_0s)^2+(1-s^2)r_0^2,\quad r\tau=\tilde{r}-r_0s.
  \end{align*}
On $\mathcal{N}^{-}(q)/\mathbf{D}$, $\tilde{r}$ and $s$ have to verify the relation
\begin{align*}
  r^2+4\geq t^2,\quad t=t_0-\tilde{r},
\end{align*}
that is,
\begin{align*}
 s\leq s_*, \quad  2\tilde{r}(t_0-r_0s_*)+4= (t_0-r_0)(t_0+r_0).
\end{align*}
As $-1\leq s\leq 1$, to make the set $\mathcal{N}^{-}(q)/\mathbf{D}$ non-empty, it in particular requires that $$\tilde{r}\geq \frac{t_0-r_0}{2}-2\geq \frac{1}{5}u_0.$$
 Here keep in mind that we have assumed that $t_0\geq 20$, $t_0-r_0\geq 10$. For the case when $\tilde{r}\geq \frac{t_0+r_0}{2}-2$, it can be showed that
\[
s_*\geq 1,\quad r\tau=\tilde{r}-r_0s\geq 0,\quad r\geq \frac{1}{2}(\tilde{r}-r_0s+\sqrt{1-s^2}r_0)\geq \frac{1}{4}(\tilde{r}-r_0+\sqrt{1-s}r_0).
\]
Therefore we can estimate that
\begin{align*}
 & \int_{\mathcal{N}^{-}(q)/\mathbf{D}\cap\{\tilde{r}\geq\frac{t_0+r_0}{2}-2\}}((1+\tau)v_+^{\ga}+u_+^{\ga})^{-p}\tilde{r}^{1-p} d\tilde{r}d\tilde{\om}\\
 &\les \int_{-1}^{1}\int_{\frac{t_0+r_0}{2}-2}^{t_0} (\tilde{r}-r_0+\sqrt{1-s}r_0 )^{-p\ga}\tilde{r}^{1-p}d\tilde{r}ds\\
  &\les t_0^{1-p} u_0 \int_{-1}^{1}(u_0+\sqrt{1-s}r_0)^{-p\ga}ds.

\end{align*}
For the case when $r_0\leq \f12 t_0$, we trivially have that
\begin{align*}
\int_{-1}^{1}(u_0+\sqrt{1-s}r_0)^{-p\ga}ds\les t_0^{-p\ga}
\end{align*}
as $u_0=t_0-r_0+1\geq \frac{1}{4}t_0$. When $r_0\geq \f12 t_0$, we show that
\begin{align*}
\int_{-1}^{1}(u_0+\sqrt{1-s}r_0)^{-p\ga}ds &\les \int_{-1}^{1-r_0^{-2}u_0^2} (\sqrt{1-s}r_0)^{-p\ga}ds+\int_{1-r_0^{-2}u_0^2}^{1} u_0^{-p\ga}ds\\
&\les r_0^{-2}u_0^{2-p\ga}+ r_0^{-p\ga}(r_0^{-2}u_0^2)^{1-\frac{p\ga}{2}}\\
&\les t_0^{-2}u_0^{2-p\ga}.
\end{align*}
Here notice that $p\ga>2$. We thus conclude that
\begin{align*}
 \int_{\mathcal{N}^{-}(q)/\mathbf{D}\cap\{\tilde{r}\geq\frac{t_0+r_0}{2}-2\}}((1+\tau)v_+^{\ga}+u_+^{\ga})^{-p}\tilde{r}^{1-p} d\tilde{r}d\tilde{\om} \les   t_0^{-1-p} u_0^{3-p\ga}.
\end{align*}
Next we consider the case when $\frac{t_0-r_0}{2}-2\leq \tilde{r}\leq \frac{t_0+r_0}{2}-2$. By the definition of $s_*$, we have
\begin{align*}
 & \int_{\mathcal{N}^{-}(q)/\mathbf{D}\cap\{\frac{t_0-r_0}{2}-2\leq \tilde{r}\leq\frac{t_0+r_0}{2}-2\}}((1+\tau)v_+^{\ga}+u_+^{\ga})^{-p}\tilde{r}^{1-p} d\tilde{r}d\tilde{\om}\\
 &\les \int_{\frac{t_0-r_0}{2}-2}^{\frac{t_0+r_0}{2}-2} \int_{-1}^{s_*} ((r+\tilde{r}-r_0s) r^{\ga-1}+u_+^{\ga})^{-p}\tilde{r}^{1-p} d\tilde{r}ds.
\end{align*}
For the case when $s_*\leq r_0^{-1}\tilde{r}$, that is,
\begin{align*}
s_*=\frac{2\tilde{r}t_0+4-t_0^2+r_0^2}{2\tilde{r}r_0}\leq r_0^{-1}\tilde{r}\Longleftrightarrow (\tilde{r}-\f12 t_0)^2\geq 2+\f12 r_0^2-\frac{1}{4}t_0^2,
\end{align*}
and for the situation $2+\f12 r_0^2\leq \frac{1}{4}t_0^2$, then
\begin{align*}
  r\geq \frac{1}{4}(\tilde{r}-r_0s+\sqrt{1-s}r_0).
\end{align*}
Let's distinguish for two cases: when $r_0$ is small compared to $t_0$, that is, $r_0\leq \frac{1}{10}t_0$, then
\begin{align*}
  r\geq \frac{1}{4}(\tilde{r}-r_0)\geq \frac{1}{100}t_0.
\end{align*}
Otherwise we have $\frac{1}{10}t_0\leq r_0\leq\frac{\sqrt{2}}{2}\sqrt{t_0^2-8} $ and then we can show that
\begin{align*}
  \tilde{r}-r_0s+\sqrt{1-s}r_0&\geq \tilde{r}-r_0+\sqrt{\frac{(t_0-r_0)(t_0+r_0-2\tilde{r})-4}{2\tilde{r}r_0}}r_0\\
  &\geq \tilde{r}-r_0+\frac{1}{20}\sqrt{r_0(\frac{t_0+r_0}{2}-\tilde{r})}\\
  &\geq \frac{1}{100}t_0,\quad \forall s\leq s_*,\quad \frac{t_0-r_0}{2}-2\leq \tilde{r}\leq \frac{t_0+r_0}{2}-2.
\end{align*}
Thus for the case when $8+2r_0^2\leq t_0^2$, we always have
\begin{align*}
 \int_{\frac{t_0-r_0}{2}-2}^{\frac{t_0+r_0}{2}-2} \int_{-1}^{s_*} ((r+\tilde{r}-r_0s) r^{\ga-1}+u_+^{\ga})^{-p}\tilde{r}^{1-p} d\tilde{r}d\tilde{\om}\les t_0^{2-p-p\ga}.
\end{align*}
Now it remains to consider the case when $2+\f12 r_0^2>\frac{1}{4}t_0^2$. For the integral on $r_0\leq \tilde{r}\leq \frac{t_0+r_0}{2}-2$, similarly we can estimate that
\begin{align*}
& \int_{r_0}^{\frac{t_0+r_0}{2}-2} \int_{-1}^{s_*} ((r+\tilde{r}-r_0s) r^{\ga-1}+u_+^{\ga})^{-p}\tilde{r}^{1-p} d\tilde{r}ds\\
&\les  \int_{r_0}^{\frac{t_0+r_0}{2}-2} \int_{-1}^{s_*} (\tilde{r}-r_0+\sqrt{1-s}r_0)^{-p\ga}t_0^{1-p} d\tilde{r}ds\\
&\les  \int_{r_0}^{\frac{t_0+r_0}{2}-2}  (1-s_*)^{1-\f12 p\ga}t_0^{1-p-p\ga} d\tilde{r}ds\\
&\les t_0^{-1-p}u_0^{1-\frac{p\ga}{2}}
\end{align*}
as $p\ga\leq p(p-1)\leq 4$.

Now we need to estimate the integral on $[\frac{t_0-r_0}{2}-2, r_0]$. And we first consider the case when $\tilde{r}\leq \f12 r_0$.
 Since $2+\f12 r_0^2\geq \frac{1}{4}t_0^2$, $t_0\geq 20$ and $s\leq s_*$, in particular we have
 \begin{align*}
 &r\geq t\geq t_0-\tilde{r}\geq t_0-\f12 r_0\geq \frac{1}{10}t_0,\\
 & r+\tilde{r}-r_0s=\frac{(1-s^2)r_0^2}{r+r_0s-\tilde{r}}\geq \frac{1}{100}(1-s)t_0,\\
 &1-s_
 =\frac{(t_0-r_0)(t_0+r_0-2\tilde{r})-4}{2\tilde{r}r_0}\geq \frac{u_0}{10 \tilde{r}}.
\end{align*}
The second inequality holds trivially when $\tilde{r}-r_0s\geq 0$. Otherwise we use the bound of $r\leq 2r_0$.
Therefore we can estimate that
\begin{align*}
  &\int^{\f12 r_0}_{\frac{t_0-r_0}{2}-2}\int_{-1}^{s_{*}}((r+\tilde{r}-r_0s) r^{\ga-1}+u_+^{\ga})^{-p}\tilde{r}^{1-p} d\tilde{r}ds\\
  &\les \int^{\f12 r_0}_{\frac{t_0-r_0}{2}-2}\int_{-1}^{s_{*}}(1-s)^{-p}t_0^{-p\ga} \tilde{r}^{1-p} d\tilde{r}ds\\
  &\les \int^{\f12 r_0}_{\frac{t_0-r_0}{2}-2} (u_0\tilde{r}^{-1})^{1-p}t_0^{-p\ga} \tilde{r}^{1-p} d\tilde{r}\\
  &\les u_0^{1-p}t_0^{1-p\ga}.
\end{align*}
Finally it remains to consider the integral on $[\f12 r_0, r_0]$ with $2+\f12 r_0^2>\frac{1}{4}t_0^2$.
Denote $$s_{**}=\min\{s_*, r_0^{-1}\tilde{r}\}.$$
For the integral restricted on $-1\leq s\leq s_{**}$, note that
\begin{align*}
   r^2=(\tilde{r}-r_0s)^2+(1-s^2)r_0^2 &\geq \f12 (\tilde{r}-r_0s_*)^2+\f12 (1-s_*^2)r_0^2+\f12 (1-s)r_0^2\\
   &\geq \f12 (t_0-r_0)^2+\f12 (1-s)r_0^2.
\end{align*}
Here recall that $s_*$ is defined such that $r=t=t_0-\tilde{r}$.
Therefore we can show that
\begin{align*}
  &\int_{\f12 r_0}^{r_0}\int_{-1}^{s_{**}}((r+\tilde{r}-r_0s) r^{\ga-1}+u_+^{\ga})^{-p}\tilde{r}^{1-p} d\tilde{r}ds\\
  &\les \int^{ r_0}_{\f12 r_0}\int_{-1}^{s_{**}}(u_0^2+(1-s)r_0^2)^{-\f12 p\ga} t_0^{1-p} d\tilde{r}ds\\
  &\les \int^{ r_0}_{\f12 r_0} (r_0^{-2}u_0^2+1-r_0^{-1}\tilde{r})^{1-\frac{p\ga}{2}} t_0^{1-p-p\ga} d\tilde{r}\\
  &\les t_0^{2-p-p\ga} (1+(r_0^{-2}u_0^2)^{2-\frac{p\ga}{2}})\\
  &\les t_0^{2-p-p\ga}
\end{align*}
as $p\ga<4$ and $u_0< r_0$.

Lastly for the integral on $s_{**}=r_0^{-1}\tilde{r}<s\leq s_*$, which in particular requires that
\begin{align*}
  \tilde{r}\leq r_*=\f12 t_0+\sqrt{2+\f12 r_0^2-\frac{1}{4}t_0^2}<r_0.
\end{align*}
Moreover
\begin{align*}
  r+\tilde{r}-r_0s=\frac{(1-s^2)r_0^2}{r+r_0s-\tilde{r}}\geq \frac{(1-s)r_0^2}{ r},\quad r\leq 2r_0,\quad \forall \tilde{r}\leq r_0.
\end{align*}
Therefore
\begin{align*}
  (r+\tilde{r}-r_0s)r^{\ga-1}\geq (1-s)r_0^{2}r^{\ga-2}\geq 2^{\ga-2}(1-s)r_0^{\ga}\geq 2^{-2}(1-s)t_0^{\ga}.
\end{align*}
This leads to
\begin{align*}
  &\int_{\f12 r_0}^{r_*}\int_{r_0^{-1}\tilde{r}}^{s_{*}}((r+\tilde{r}-r_0s) r^{\ga-1}+u_+^{\ga})^{-p}\tilde{r}^{1-p} d\tilde{r}ds\\
  &\les \int^{ r_*}_{\f12 r_0}\int_{r_0^{-1}\tilde{r}}^{s_{*}}(1+(1-s)t_0^{\ga})^{-p} t_0^{1-p} d\tilde{r}ds\\

  &\les t_0^{1-p}\int^{ r_*}_{\f12 r_0}t_0^{-\ga}(1+(1-s_*)t_0^{\ga})^{1-p}  d\tilde{r} \\
  &\les t_0^{1-p}\int^{ r_*}_{\f12 r_0}t_0^{-\ga}(1+t_0^{\ga-2}u_0 (t_0+r_0-2\tilde{r}))^{1-p}  d\tilde{r} \\
  &\les  t_0^{3-p-2\ga}u_0^{-1}

\end{align*}
Since $2<p<\frac{1+\sqrt{17}}{2}$, $1<\ga<p-1$ and $u_0<t_0$, gathering all the above estimates, we have shown that
\begin{align*}
  |\int_{\mathcal{N}^{-}(q)/\mathbf{D}}|\phi|^{p-1} \phi \ \tilde{r} d\tilde{r}d\tilde{\om}|
  \les (\mathcal{E}_{0, \ga_0}[\phi])^{\frac{p}{p+1}} t_0^{\frac{3-p-2\ga}{p+1}}u_0^{-\frac{1}{p+1}}.
\end{align*}
Now we compute that
\begin{align*}
  3-p-2\ga+\frac{p^2-4p+7}{5-p}\g
  =-\frac{3-p}{5-p}(2p-4+(p+1)(\ga-1))< \frac{9-p^2}{5-p} (1-\ga)<1-\ga_0
\end{align*}
by choosing $\ga$ sufficiently close to $\ga_0$. This demonstrates that
\begin{align*}
  |\phi^{lin}_H(t_0, x_0)|\les t_0^{\frac{3-p-2\ga}{p+1}}u_0^{-\frac{1}{p+1}}\les t_0^{-\frac{(p-2)^2+3}{(p+1)(5-p)}\ga_0}u_0^{-\frac{\ga_0}{p+1}}.
\end{align*}
This proves \eqref{eq:phi:pt:Br:smallp:lin} and hence we finished the proof for the Proposition.

\iffalse
If $u_0^2\leq t_0^{2-\ga}$, then
\begin{align*}
  t_0^{-\frac{3-p}{5-p}(2p-4+(p+1)(\ga-1))}u_0^{-1}&\leq u_0^{-1-\frac{2}{2-\ga}\frac{3-p}{5-p}(2p-4+(p+1)(\ga-1))}\\
  &=u_0^{-1+\frac{2(p+1)(3-p)}{(5-p)}-\frac{2}{2-\ga}\frac{(3p-3)(3-p)}{(5-p)}}\\
  &\leq u_0^{-1+\frac{2(p+1)(3-p)}{(5-p)}-\frac{2(3p-3)(3-p)}{(5-p)}}\\
  &\leq u_0^{-1-\frac{4(p-2)(3-p)}{(5-p)}}
\end{align*}
computations
\begin{align*}
  (3-p)(p+3)-5+p=4-p^2+p>0,\quad 2<p<\frac{1+\sqrt{17}}{2}.
\end{align*}

\begin{align*}
  &\int^{\f12 r_0}_{\frac{t_0-r_0}{2}-2}\int_{-1}^{s_{**}}((r+\tilde{r}-r_0s) r^{\ga-1}+u_+^{\ga})^{-p}\tilde{r}^{1-p} d\tilde{r}ds\\
  &\les \int^{\f12 r_0}_{\frac{t_0-r_0}{2}-2}\int_{-1}^{s_{**}}(\sqrt{1-r_0^{-1}\tilde{r}}r_0)^{-p}\tilde{r}^{1-p} d\tilde{r}ds\les t_0^{-p\ga}u_0^{2-p}.
\end{align*}
As  and $\tilde{r}\leq \f12 r_0$, we have
\fi

\end{proof}

\section{Semilinear wave equation on a truncated backward light cone}
\label{sec:comp}

We study the solution to a class of semilinear wave equations on a compact region with smooth initial data $(\phi_0, \phi_1)$ which may blow up on the boundary. This is motivated by studying the asymptotic behaviour of solutions to subcritical defocusing nonlinear wave equation in the interior region. However the content in this section is independent and may be of independent interest.

Let $R>1$ be a constant and $\B_R$ be the ball with radius $R$ in $\mathbb{R}^3$. Denote $\cJ^+(\B_R)$ be the future maximal Cauchy development, that is, $(t, x)\in\mathbb{R}\times \mathbb{R}^{3}$ belongs $\cJ^+(\B_R)$ if and only if $x+t\om\in \B_R$ for all $\om\in \mathbb{S}^2$.
Consider the Cauchy problem to the following nonlinear wave equation
\begin{equation}
  \label{eq:NLW:3D:conf}
  \begin{cases}
    \Box \phi= \La^{3-p}|\phi|^{p-1}\phi,\\
    \phi(0, x)=\phi_0,\quad \pa_t\phi(0, x)=\phi_1, \quad x\in \B_R
  \end{cases}
\end{equation}
on $\mathcal{J}^{+}(\B_R)$ with $\La=((R-t)^2-|x|^2)^{-1}$.

For any fixed point $q=(t_0, x_0)\in \mathcal{J}^+(\B_R)$, recall that $\mathcal{N}^{-}(q)$ is the past null cone of the point $q$ in $\mathcal{J}^{+}(\B_R)$ (as pointed out before, we are only concerned with the solution in the future) and $\mathcal{J}^{-}(q)$ is the past of the point $q$, that is, the region bounded by $\mathcal{N}^{-}(q)$ and $\B_R$. As defined in Section \ref{sec:notation}, the tilde coordinates, quantities are referred to those ones with coordinates centered at the given point $q=(t_0, x_0)$.

At the fixed point $q=(t_0, x_0)$, define the following functions
\[
u_*=R-t+r, \quad v_*=R-t-r,\quad \tau=\frac{x\cdot (x-x_0)}{|x||x-x_0|}=\om\cdot \tilde{\om}.
\]
In particular $\La=u_*^{-1}v_*^{-1}$.

 Assume that the initial data $(\phi_0, \phi_1)$ are bounded in  the following weighted energy norm
\[
\tilde{\mathcal{E}}_{0, \gamma}=\int_{\B_R}(R-|x|)^{\gamma}|L\phi|^2+|{\Lb}\phi|^2+|\nabb\phi_0|^2+|\phi_0|^2+(R-|x|)^{p-3+\ga}|\phi_0|^{p+1} dx
\]
for some constant $0<\ga<1$. Define
\begin{align*}
  \mathcal{I}=\iint_{\cJ^+(\B_R)} \La^{3-p}|\phi|^{p+1}v_{*}^{\ga-1} d\vol.
\end{align*}
First we establish a weighted energy flux bound for the potential.
\begin{Prop}
\label{prop:EF:cone:gamma}
For any point $q=(t_0, x_0)\in \mathcal{J}^{+}(\B_R)$, the solution verifies the uniform bound
\begin{equation}
\label{eq:comp:v:EF}
\begin{split}
&\int_{\mathcal{N}^{-}(q)} ( v_*^\ga+(1-\tau)u_*^{\ga})   \La^{3-p} |\phi|^{p+1} \leq C (\tilde{\cE}_{0, \ga}+\mathcal{I})
\end{split}
\end{equation}
for some constant $C$ depending only on $R$, $p$ and $\ga$.
\end{Prop}
\begin{proof}
The proof is similar to that for Proposition \ref{prop:EF:cone:NW:3d}. For solution $\phi$ of the equation \eqref{eq:NLW:3D:conf}, define the associated energy momentum tensor
\begin{align*}
  T[\phi]_{\mu\nu}=\pa_{\mu}\phi\pa_{\nu}\phi-\f12 m_{\mu\nu}(\pa^\ga \phi \pa_\ga\phi+\frac{2}{p+1}\La^{3-p}|\phi|^{p+1}).
\end{align*}
Here $m_{\mu\nu}$ is the flat Minkowski metric.
Then
\begin{align*}
  \pa^\mu T[\phi]_{\mu\nu}=&(\Box\phi-\La^{3-p} |\phi|^{p-1}\phi)\pa_\nu\phi+\frac{p-3}{p+1}\La^{2-p}\pa_\nu\La |\phi|^{p+1}.
\end{align*}
Recall that current $J^{X, \chi}[\phi]$ defined for any vector field $X$ and any function $\chi$
\begin{equation*}
J^{X, \chi}_\mu[\phi]=T[\phi]_{\mu\nu}X^\nu -
\f12\pa_{\mu}\chi \cdot|\phi|^2 + \f12 \chi\pa_{\mu}|\phi|^2.
\end{equation*}
For solution $\phi$ of equation \eqref{eq:NLW:3D:conf}, we derive the following energy identity
\begin{equation}
\label{eq:energy:id:conf}
\iint_{\mathcal{D}}\pa^\mu  J^{X,\chi}_\mu[\phi] d\vol =\iint_{\mathcal{D}}\frac{p-3}{p+1}\La^{2-p}X(\La) |\phi|^{p+1}+ T[\phi]^{\mu\nu}\pi^X_{\mu\nu}+
\chi \pa_\mu\phi\pa^\mu\phi -\f12\Box\chi\cdot|\phi|^2+\chi\phi\Box\phi d\vol
\end{equation}
for any domain $\mathcal{D}$ in $\mathcal{J}^+(\B_{R})$.

Apply the above energy identity to the domain $\mathcal{J}^{-}(q)$ for some fixed point $q=(t_0, x_0)$ with vector field $X$
\[
X=(R-t-r)^{\gamma} L+(R-t+r)^\gamma \Lb=v_*^\ga L+u_*^\ga \Lb.
\]
Since the bulk integral on the right hand side is an integral over a spacetime region, we can compute it under the null frame $(L, \Lb, e_1, e_2)$. We first can compute that
\[
\nabla_{L}X=-2\gamma v_*^{\gamma-1},\quad \nabla_{\Lb}X=-2\gamma u_*^{\gamma}\Lb,\quad \nabla_{e_i}X=r^{-1}(v_*^\gamma-u_*^{\gamma}) e_i.
\]
In particular, the non-vanishing components of the deformation tensor $\pi_{\mu\nu}^X$ are
\[
\pi^X_{L\Lb}=2\gamma \left(v_*^{\gamma-1}+u_*^{\gamma-1}\right),\quad \pi^X_{e_i e_i}=r^{-1}(v_*^{\ga}-u_*^{\gamma}).
\]
Therefore we have
\begin{align*}
&T[\phi]^{\mu\nu}\pi^X_{\mu\nu}=2T[\phi]^{L\Lb}\pi^X_{L\Lb}+T[\phi]^{e_ie_i}\pi^X_{e_ie_i}\\
&=\ga (v_*^{\gamma-1}+u_*^{\gamma-1})(|\nabb\phi|^2+\frac{2}{p+1}\La^{3-p}|\phi|^{p+1})+r^{-1}(v_*^{\ga}-u_*^{\gamma})(|\nabb\phi|^2-\pa^\mu \phi \pa_\mu\phi-\frac{2}{p+1}\La^{3-p}|\phi|^{p+1})\\
&=\left(\ga(v_*^{\gamma-1}+u_*^{\gamma-1})+r^{-1}(v_*^{\ga}-u_*^{\gamma})\right)|\nabb\phi|^2-r^{-1}(v_*^{\ga}-u_*^{\gamma})\pa^\mu\phi \pa_\mu\phi\\
&\quad +\left(\ga(v_*^{\gamma-1}+u_*^{\gamma-1})+r^{-1}(u_*^{\ga}-v_*^{\gamma})\right)\frac{2}{p+1}\La^{3-p}|\phi|^{p+1}.
\end{align*}
Now take the function $\chi$ to be
\[
\chi=r^{-1}(v_*^{\ga}-u_*^{\gamma}).
\]
We may note that
\begin{align*}
\Box \chi=-r^{-1}L\Lb (r\chi)=-2r^{-1}L\Lb(v_*^\ga-u_*^\ga)=0,\quad r>0.
\end{align*}
Moreover we can compute that
\begin{align*}
  X(\La)=(u_*^\ga \Lb+v_*^\ga L)(u_*v_*)^{-1}=-(u_*v_*)^{-2}(-2 u_*^{\ga}v_*-2v_*^{\ga}u_*)=2\La^{2}(u_*^{\ga}v_*+v_*^{\ga}u_*),
\end{align*}
Therefore for solution $\phi$ to \eqref{eq:NLW:3D:conf}, we have
\begin{align*}
&T[\phi]^{\mu\nu}\pi^X_{\mu\nu}+
\chi \pa_\mu\phi \pa^\mu\phi -\f12\Box\chi\cdot|\phi|^2 +\chi\phi\Box\phi +\frac{p-3}{p+1}\La^{2-p}X(\La) |\phi|^{p+1}\\
&=\left(\ga(v_*^{\gamma-1}+u_*^{\gamma-1})+r^{-1}(v_*^{\ga}-u_*^{\gamma})\right)|\nabb\phi|^2+r^{-1}(v_*^{\ga}-u_*^{\gamma}) \La^{3-p}|\phi|^{p+1}\\
&\quad +\left(\ga(v_*^{\gamma-1}+u_*^{\gamma-1})+r^{-1}(u_*^{\ga}-v_*^{\gamma})\right)\frac{2}{p+1}\La^{3-p}|\phi|^{p+1}+2\frac{p-3}{p+1}\La^{4-p}(u_*v_*^{\ga}+u_*^{\ga}v_*) |\phi|^{p+1}\\
&=\left(\ga(v_*^{\gamma-1}+u_*^{\gamma-1})-r^{-1}(u_*^{\ga}-v_*^{\gamma})\right)(|\nabb\phi|^2+\frac{2}{p+1}\La^{3-p}|\phi|^{p+1})\\
&\quad + \La^{3-p}|\phi|^{p+1} \frac{p-3}{p+1}(2(v_*^{\gamma-1}+u_*^{\gamma-1})-r^{-1}(u_*^{\ga}-v_*^{\gamma}))
\end{align*}
Now note that $v_*\geq 0$, $u_*\geq 0$ and $2r=u_*-v_*$. Define $f(u)=u^{\ga}$ for $u>0$. As $0< \ga<1$, we conclude that the function $f'(u)=\ga u^{\ga-1}$ is convex. Therefore
\begin{align*}
\frac{f(u_*)-f(v_*)}{u_*-v_*}=\int_{0}^{1}f'(su_*+(1-s)v_*)ds \leq \int_{0}^{1}sf'(u_*)+(1-s)f'(v_*)ds=\frac{f'(u_*)+f'(v_*)}{2},
\end{align*}
which implies that
\[
\ga(v_*^{\gamma-1}+u_*^{\gamma-1})+r^{-1}(v_*^{\ga}-u_*^{\gamma})\geq 0.
\]
Since $0< \ga<1$, we conclude that
\[
2(v_*^{\gamma-1}+u_*^{\gamma-1})-r^{-1}(u_*^{\ga}-v_*^{\gamma})\geq 0.
\]
Therefore the bulk integral on the right hand side of the energy identity \eqref{eq:energy:id:conf} is nonnegative for the super-conformal case $p\geq 3$. Otherwise, we make use of the priori bound $\mathcal{I}$.  This leads to the following energy estimate
\begin{equation}
\label{eq:positive:ga:comp}
\begin{split}
\iint_{\mathcal{J}^{-}(q)}\pa^\mu J^{X, \chi}_\mu[\phi]d\vol&=\int_{\mathcal{N}^{-}(q)}i_{ J^{X, \chi}[\phi]}d\vol+\int_{\mathcal{J}^{-}(q)\cap{\B_R}}i_{ J^{X, \chi}[\phi]}d\vol\\
&\geq - \frac{4|p-3|}{p+1}\iint_{\cJ^{-}(q)}\La^{3-p}|\phi|^{p+1}  v_*^{\gamma-1} d\vol\\
&\geq -\frac{4|p-3|}{p+1}\mathcal{I}
\end{split}
\end{equation}
by using Stokes' formula. Here we note that $0<\ga<1$ and $u_*\geq v_*$. The boundary integral on the initial hypersurface $\mathcal{J}^{-}(q)\cap{\B_R}$ can be bounded by the initial data. The above inequality then gives control on the weighted energy flux through the backward light cone $\mathcal{N}^{-}(q)$. To find the explicit form of this weighted energy flux, we shift to the coordinates centered at the point $q=(t_0, x_0)$.
 Recall from the proof for Proposition \ref{prop:EF:cone:NW:3d} that
\begin{align*}
-i_{J^{X,\chi}[\phi]}d\vol=J_{\tilde{\Lb}}^{X,\chi}[\phi]\tilde{r}^2d\tilde{u}d\tilde{\om}= ( T[\phi]_{\tilde{\Lb}\nu}X^\nu -
\f12(\tilde{\Lb}\chi) |\phi|^2 + \f12 \chi\cdot\tilde{\Lb}|\phi|^2 ) \tilde{r}^2d\tilde{u}d\tilde{\om}.
\end{align*}
For the main quadratic terms, we first can compute that
\begin{align*}
T[\phi]_{\tilde{\Lb}\nu}X^\nu =T[\phi]_{\tilde{\Lb}\tilde{\Lb}}X^{\tilde{\Lb}}+T[\phi]_{\tilde{\Lb}\tilde{L}}X^{\tilde{L}}+T[\phi]_{\tilde{\Lb}\tilde{e}_i}X^{\tilde{e}_i}.
\end{align*}
 We expand the vector field $X$ under the new null frame $\{\tilde{L}, \tilde{\Lb}, \tilde{e}_1, \tilde{e}_2\}$ centered at the point $q$. Recall those computations in the proof for Proposition \ref{prop:EF:cone:NW:3d}, we can write that
\begin{align*}
X
&=\f12 \left(u_*^\ga+v_*^\ga+(v_*^\ga-u_*^\ga)\om\cdot \tilde{\om}\right)\tilde{L}+\f12 \left(u_*^\ga+v_*^\ga-(v_*^\ga-u_*^\ga)\om\cdot \tilde{\om}\right)\tilde{\Lb}+(v_*^\ga-u_*^\ga)\om\cdot \tilde{\nabb}.
\end{align*}
Here $\tilde{\nabb}=\tilde{\nabla}-\pa_{\tilde{r}}$. Denote $\tau=\om\cdot \tilde{\om}$. Then we can compute the quadratic terms
\begin{align*}
T[\phi]_{\tilde{\Lb}\nu}X^\nu
=&\left((1-\tau)v_*^\ga+(1+\tau)u_*^\ga\right)|{\tilde{\Lb}}\phi|^2 +\left((1+\tau)v_*^\ga+(1-\tau)u_*^\ga\right)(|\tilde{\nabb}\phi|^2+\frac{\La^{3-p}}{p+1}|\phi|^{p+1})\\
&+2 (v_*^\ga-u_*^\ga) ({\tilde{\Lb}}\phi) (\om\cdot \tilde{\nabb})\phi.
\end{align*}
Similar to the proof of Proposition \ref{prop:EF:cone:NW:3d}, we write the above quantity in terms of $r\phi$ and show that the quadratic terms are nonnegative.
Indeed since $u_*^\ga\geq v_*^\ga$, we therefore can bound that
\begin{align*}
&\left((1-\tau)v_*^\ga+(1+\tau)u_*^\ga\right)|\tilde{\Lb}(r\phi)|^2+\left((1+\tau)v_*^\ga+(1-\tau)u_*^\ga\right)|\tilde{\nabb}(r\phi)|^2+2 (v_*^\ga-u_*^\ga){\tilde{\Lb}}(r\phi) (\om \cdot \tilde{\nabb})(r\phi)\\
&\geq \left((1-\tau)v_*^\ga+(1+\tau)u_*^\ga\right)|\tilde{\Lb}(r\phi)|^2+\left((1+\tau)v_*^\ga+(1-\tau)u_*^\ga\right)|\tilde{\nabb}(r\phi)|^2\\
  &\quad -2 (u_*^\ga-v_*^\ga)\sqrt{1-\tau^2}|{\tilde{\Lb}}(r\phi)| |\tilde{\nabb}(r\phi)|\\
&\geq \frac{2u_*^\ga v_*^\ga}{(1-\tau)u_*^\ga+(1+\tau)v_*^\ga} |\tilde{\Lb}(r\phi)|^2 +\frac{2u_*^\ga v_*^\ga}{(1+\tau)u_*^\ga+(1-\tau)v_*^\ga}|\tilde{\nabb}(r\phi)|^2 \geq 0.
\end{align*}
For the other lower order terms, we compute  that
\begin{align*}
&\left((1-\tau)v_*^\ga+(1+\tau)u_*^\ga\right)(\tau^2|\phi|^2+2{\tilde{\Lb}}(r\phi)\tau\phi )+(r\chi)( {\tilde{\Lb}}(r\phi)+\tau \phi) \phi-\f12(r\tilde{\Lb}(r\chi)+\tau r\chi)|\phi|^2\\
&+\left((1+\tau)v_*^\ga+(1-\tau)u_*^\ga\right)\big((1-\tau^2)|\phi|^2-2(\om-\tilde{\om}\tau)\tilde{\nabb}(r\phi) \phi\big)\\
&+(v_*^\ga-u_*^\ga)\big(-2\tau(1-\tau^2)|\phi|^2+2\phi( \tau(\om\cdot \tilde{\nabb})(r\phi)-(1-\tau^2){\tilde{\Lb}}(r\phi))\big)\\
&=(-\f12 r\tilde{\Lb}(r\chi)+v_*^\ga+u_*^\ga)|\phi|^
+2(v_*^\ga+u_*^\ga)(\tau {\tilde{\Lb}}-\om\cdot \tilde{\nabb})(r\phi) \phi\\
&=-r^2\tilde{r}^{-1} \tilde{\Om}_{ij}(r^{-3}(v_*^\ga+u_*^\ga) \om_j\tilde{\om}_i |r\phi|^2)+\tilde{r}^{-2}r^2\tilde{\Lb}(r^{-1}\tau\tilde{r}^2(v_*^\ga+u_*^\ga) |\phi|^2)\\
&+(-\f12 r\tilde{\Lb}(r\chi)+v_*^\ga+u_*^\ga)|\phi|^2-\tilde{r}^{-2}r^2\tilde{\Lb}(r^{-3}\tau\tilde{r}^2(v_*^\ga+u_*^\ga)) |r\phi|^2+r^2 \tilde{r}^{-1}\tilde{\Om}_{ij}(r^{-3}(v_*^\ga+u_*^\ga) \om_j\tilde{\om}_i) |r\phi|^2.
\end{align*}
Similarly we can compute that
\begin{align*}
&\tilde{r}^{-1}\tilde{\Om}_{ij}(r^{-3}\om_j\tilde{\om}_i)=-2r^{-4}(1-2\tau^2)-2\tau \tilde{r}^{-1}r^{-3},\\
&\tilde{r}^{-2}r^4\tilde{\Lb}(r^{-3}\tilde{r}^2\tau)=4\tau^2-1-2r\tilde{r}^{-1}\tau.
\end{align*}
Thus the coefficients of $|\phi|^2$ in the last line in the previous equation verify
\begin{align*}
&(-\f12 r\tilde{\Lb}(r\chi)+v_*^\ga+u_*^\ga)-\tilde{r}^{-2}r^4\tilde{\Lb}(r^{-3}\tau\tilde{r}^2(v_*^\ga+u_*^\ga)) +r^4 \tilde{r}^{-1} \tilde{\Om}_{ij}(r^{-3}(v_*^\ga+u_*^\ga) \om_j\tilde{\om}_i) \\
&=-r(\pa_t-\tilde{\om}\cdot \nabla)(v_*^\ga-u_*^\ga)-r\tau (\pa_t-\tilde{\om}\cdot \nabla)(u_*^\ga+v_*^\ga)+r(\pa_r-\tau \tilde{\om}\cdot \nabla)(u_*^\ga+v_*^\ga)\\
&\quad +(u_*^\ga+v_*^\ga)\left(1-(4\tau^2-1-2r\tilde{r}^{-1}\tau)-2(1-2\tau^2)+2\tau \tilde{r}^{-1}r\right)\\
&=r(\pa_t+\pa_r)u_*^\ga+r(\pa_r-\pa_t)v_*^\ga-\tau r(\pa_t+\pa_r)u_*^\ga-\tau r(\pa_t-\pa_r)v_*^\ga=0.
\end{align*}
 By using integration by parts on the backward light cone $\mathcal{N}^{-}(q)$, we derive that
\begin{align*}
  &\int_{\mathcal{N}^{-}(q)}\big(-r^2\tilde{r}^{-1} \tilde{\Om}_{ij}(r^{-3}(v_*^\ga+u_*^\ga) \om_j\tilde{\om}_i |r\phi|^2)+\tilde{r}^{-2}r^2\tilde{\Lb}(r^{-1}\tau\tilde{r}^2(v_*^\ga+u_*^\ga) |\phi|^2)\big)r^{-2}\tilde{r}^2 d\tilde{u}d\tilde{\om}\\
  &= \int_{\mathcal{N}^{-}(q)\cap \B_R}r^{-1}\tau \tilde{r}^2 (u_*^\ga+v_*^\ga)|\phi|^2d\tilde{\om}.
\end{align*}
To summarize, the above computations show that
\begin{equation}
\label{eq:EST:Nq:comp}
\begin{split}
&\int_{\mathcal{N}^{-}(q)}((1+\tau)v_*^\ga+(1-\tau)u_*^\ga ) \frac{\La^{3-p}}{p+1}|\phi|^{p+1} \tilde{r}^2d\tilde{u}d\tilde{\om}\\
&\leq -\int_{\mathcal{N}^{-}(q)}i_{J^X[\phi, F]}d\vol+ \int_{\mathcal{N}^{-}(q)\cap \B_R}r^{-1}\tau \tilde{r}^2 (u_*^\ga+v_*^\ga)|\phi|^2d\tilde{\om}.
\end{split}
\end{equation}
We next compute the boundary integral on $\B_R\cap \mathcal{J}^{-}(q)$ in the inequality \eqref{eq:positive:ga:comp} which we compute it  under the coordinates system $(t, x)$. As the initial hypersurface $\B_R$ has the volume form $dx$, the contraction reads
\begin{align*}
i_{J^{X, \chi}[\phi ]}d\vol&= T[\phi ]_{0 L}X^L+T[\phi ]_{\Lb 0}X^{\Lb}- \f12 \pa_t\chi |\phi|^2+ \f12\chi \pa_t|\phi|^2\\
&= \f12 v_*^\ga( |L\phi|^2+|\nabb\phi|^2+\frac{2\La^{3-p}}{p+1}|\phi|^{p+1})-\f12\pa_t\chi \cdot |\phi|^2+\f12\chi \pa_t|\phi|^2 \\
&\quad +\f12 u_*^\ga(|{\Lb}\phi|^2+|\nabb\phi|^2+\frac{2\La^{3-p}}{p+1}|\phi|^{p+1})\\
&=\f12(u_*^\ga+v_*^\ga)( |\nabb\phi|^2+\frac{2\La^{3-p}}{p+1}|\phi|^{p+1})+\f12 v_*^\ga r^{-2}|L(r\phi)|^2\\
&\quad +\f12 u_*^\ga r^{-2}|{\Lb}(r\phi)|^2-\div(\om r^{-1}|\phi|^2(u_*^\ga+v_*^\ga)).
\end{align*}
Here $\om=\frac{x}{|x|}$ can be viewed as a vector on $\mathbb{R}^{3}$ and the divergence is taken over the initial hypersurface $\B_{R}$. The integral of the divergence term and be computed by using integration by parts. Under the coordinates $\tilde{x}=x-x_0$ on the initial hypersurface, we have
\begin{align*}
\int_{\mathcal{J}^{-}(q)\cap \B_R} \div (\om r^{-1}|\phi|^2(u_*^\ga+v_*^\ga))dx
&=\int_{\mathcal{N}^{-}(q)\cap \B_R} \tilde{r}^2 \tilde{\om} \cdot\om r^{-1}|\phi|^2(u_*^\ga+v_*^\ga)d\tilde{\om}.
\end{align*}
This term cancels the one from the integral on $\mathcal{N}^{-}(q)$ in the estimate \eqref{eq:EST:Nq:comp}.
By our assumption on the initial data, we thus bound that
\begin{align*}
&\int_{\mathcal{J}^{-}(q)\cap \B_R}i_{J^{X, \chi}[\phi ]}d\vol+\int_{\mathcal{N}^{-}(q)\cap \B_R} \tilde{r}^2 \tau r^{-1}|\phi|^2(u_*^\ga+v_*^\ga)d\tilde{\om}\\
&=\f12 \int_{\mathcal{J}^{-}(q)\cap B_R}  (u_*^\ga+v_*^\ga)( |\nabb\phi|^2+\frac{2\La^{3-p}}{p+1}|\phi|^{p+1})+  v_*^\ga r^{-2}|L(r\phi)|^2+  u_*^\ga r^{-2}|{\Lb}(r\phi)|^2 dx\\
&\leq C_R\int_{\B_R}(R-|x|)^\ga | L\phi|^2 + |\nabb\phi|^2+|{\Lb}\phi|^2 +|\phi|^2+(R-|x|)^{p-3}|\phi|^{p+1}dx\\
&\leq C_R \tilde{\cE}_{0, \ga}
\end{align*}
for some constant $C_R$ depending only on $R$. In particular this constant is independent of the choice of the point $q$. Now combining estimates \eqref{eq:positive:ga:comp}, \eqref{eq:EST:Nq:comp} , we derive the estimate \eqref{eq:comp:v:EF} of the Proposition in view of the inequality
\[
(1+\tau)v_*^\ga+(1-\tau)u_*^\ga\geq \frac{1}{2} (1-\tau)u_*^{\ga}+v_*^{\ga}.
\]
\end{proof}
\iffalse

Like the case in the exterior region, we rely on the representation formula for linear wave equation to obtain the pointwise decay estimates for the solution based on the weighted energy flux bound in the previous Proposition. The solution consists the linear part which evolves from the initial data and the nonlinear part. However the initial data are not bounded even in the energy space. We instead use energy method to control the linear evolution.

For this purpose, consider the following Cauchy problem
\begin{equation}
\label{eq:lineareq:com}
\Box w=0,\quad w(0, x)=0,\quad \pa_t w(0, x)=w_1,\quad |x|\leq R.
\end{equation}
For $0<\ga<1$ and any function $w(x)$ of $x$, denote
\begin{align*}
\mathcal{E}_{\ga}[w](\B_{R})=\int_{|x|\leq R}|w|^2(R-|x|)^\ga +\sum\limits_{i,j}|\Om_{ij}w|^2(R-|x|)^\ga+|\nabla w|^2(R-|x|)^{2+\ga}dx,
\end{align*}
where $\Om_{ij}=x_i\pa_j-x_j\pa_i$ are the angular momentum vector fields.
We have the following estimate for the linear solution.
\begin{Prop}
\label{prop:lineardecay:comp}
Consider the Cauchy problem to the linear wave equation \eqref{eq:lineareq:com}. Then there exists constant $C$ depending only on $\ga$ such that
\begin{equation}
\label{eq:lineardecay:comp}
|w(t, x)|^2\leq C  (R-t)^{-1-\ga} \mathcal{E}_{ \ga}[w_1](\B_R),\quad \forall (t, x)\in \mathcal{J}^+(\B_{R}).
\end{equation}
In particular we note that the constant $C$ is independent of the radius $R$.
\end{Prop}
To prove this Proposition, we first establish a type of Hardy's inequality.
\begin{Lem}
\label{lem:Hardy:com}
For $0<\ga<1$ and $\ep>0$, we have the following type of Hardy's inequality
\begin{equation*}
\label{eq:Hardy:com}
\int_{\B_{R}}|\varphi|^2 (R-|x|)^{-1+\ep}dx\leq C R^{\ep+\ga-1} \int_{\B_{R}}(|\varphi|^2R^{-2}+|\nabla \varphi|^2)(R-|x|)^\ga dx
\end{equation*}
for some constant $C$ depending only on $\ep$ and $\ga$.
\end{Lem}

\begin{proof}
 Note that the inequality is scaling invariant in terms of $R$ by the scaling transformation $\varphi(x)=\varphi_1(x/R)$. It hence suffices to prove the lemma for the case when $R=1$.
Moreover the inequality is homogeneous. Without loss of generality we assume that
\[
\int_{\B_{1}}(|\varphi|^2+|\nabla \varphi|^2)(1-|x|)^\ga dx=1.
\]
Here $\B_{1}$ denotes the ball with radius $1$ on the initial hypersurface.
We can estimate that
\begin{align*}
&\int_{\om}|\varphi|^2(r\om)(1-r)^{\ep}r^4 d\om +\ep\int_{0}^{r}\int_{\om}|\varphi|^2 (1-s)^{\ep-1}s^4 dsd\om \\
&=\int_{0}^{r}\int_{\om}|\varphi|^2 4(1-s)^{\ep}s^3+2 \varphi \cdot \pa_r\varphi (1-s)^{\ep}s^4 \quad dsd\om\\
&\leq \int_{0}^{r}\int_{\om} 64\ep^{-1}|\varphi|^2 (1-s)^{\ga}s^2+16\ep^{-1}|\nabla\varphi|^2(1-s)^\ga s^2+\f12 \ep |\varphi|^2(1-s)^{2\ep-\ga}s^4 \quad dsd\om.
\end{align*}
Since $\ga<1$ and $\ep>0$, take $r=1$. We derive from the above estimate that
\begin{align*}
\ep\int_{0}^{1}\int_{\om}|\varphi|^2 (1-s)^{\ep-1}s^4 dsd\om \leq 160\ep^{-1}.
\end{align*}
On the other hand, we have
\begin{align*}
\int_{|x|\leq \f12}|\varphi|^2(1-|x|)^{\ep-1}dx \leq 2^{1+\ga-\ep}\int_{|x|\leq \f12}|\varphi|^2(1-|x|)^{\ga} dx \leq 2^{1+\ga-\ep}.
\end{align*}
Combining the above two bounds, we conclude that
\begin{align*}
\int_{|x|\leq 1}|\varphi|^2(1-|x|)^{\ep-1}dx\leq 2^{1+\ga-\ep}+640\ep^{-2}.
\end{align*}
Hence the lemma holds.
\end{proof}
We also need the following weighted Sobolev embedding.
\begin{Lem}
\label{lem:weightedSob:com}
Let $\varphi$ be a smooth function defined on the ball $\B_{R}$. Then we have the weighted Sobolev embedding
\begin{equation}
\label{eq:weighteSob:com}
|\varphi|^2(x)\leq C R^{-1-\ga}\mathcal{E}_0 ,\quad \forall |x|\leq R
\end{equation}
for some constant $C$ depending only on $\ga$, where $\mathcal{E}_0=\int_{\B_{R}}|\varphi|^2(R-|x|)^{\ga}R^{-2}dx+\mathcal{E}_{ \ga}[\nabla\varphi](\B_{R})$.
\end{Lem}

\begin{proof}
Like the previous Lemma, the inequality is scaling invariant. Hence it is sufficient to prove this Lemma for the case when $R=1$ and
\[
\mathcal{E}_0=\int_{B_{1}}(|\varphi|^2+|\nabla\varphi|^2+|\nabla\Om_{ij}\varphi|^2)(1-|x|)^{\ga}+|\nabla\nabla\varphi|^2(1-|x|)^{2+\ga}dx=1
\]
For small $|x| \leq \frac{1}{2}$, we can use the standard Sobolev embedding to conclude that
\begin{align*}
|\varphi|^2(x)\leq C \|\varphi\|_{H^2(\B_{\f12})}^2\leq 2^{2+\ga}C.
\end{align*}
Here $C$ is the universal constant from the embedding inequality (however it may very in different places). For large $|x|$ we make use of the better decay of the angular derivative. From the previous Lemma \ref{lem:Hardy:com}, we derive the better decay of $\varphi$ as well as $\pa_\om\varphi$
\begin{equation}
\label{eq:bd4varphiom}
\int_{|x|\leq 1}(|\varphi|^2+|\pa_\om\varphi|^2)(1-|x|)^{\ep-1}dx\leq C_{\ep, \ga}
\end{equation}
for some constant $C_{\ep, \ga}$ depending only on $\ep$ and $\ga$. Here $\pa_\om$ denotes derivatives of the angular momentum $\Om_{ij}$. By using the Poincar\'e inequality we can derive that
\begin{align*}
\int_{\om}|\varphi|^6d\om \leq C\int_{\om}|\varphi|^2+|\pa_\om\varphi|^2d\om \cdot \int_{\om}|\varphi|^4d\om
\end{align*}
for some universal constant $C$. For some fixed constant $r_0\in[\frac{1}{3}, \f12]$ and $r\geq r_0$, integrate from the sphere with radius $r_0$. We can derive that
\begin{align*}
\int_{\om}|\varphi|^4(r\om)d\om&=\int_{\om}|\varphi|^4(r_0\om)d\om+4\int_{r_0}^{r}\pa_r\varphi \cdot \varphi^3 d\om dr\\
&\leq \int_{\om}|\varphi|^4(r_0\om)d\om+C\left(\int_{r_0}^{r}\int_{\om}|\pa_r\varphi|^2(1-|x|)^\ga drd\om \right)^\f12\left(\int_{r_0}^{r}\int_{\om}|\varphi|^6(1-|x|)^{-\ga} drd\om \right)^\f12\\
&\leq \int_{\om}|\varphi|^4(r_0\om)d\om+C \sup\limits_{r_0\leq s\leq r}\left(\int_{\om}|\varphi|^4(s\om)d\om\right)^\f12 \left(\int_{r_0}^{r}\int_{\om}(|\varphi|^2+|\pa_\om \varphi|^2)(1-|x|)^{-\ga} drd\om \right)^\f12
\end{align*}
for some universal constant $C$ and all $r_0 \leq r\leq 1$. By choosing $\ep=1-\ga$ in estimate \eqref{eq:bd4varphiom} and $r_0=\f12$ in the above inequality, we derive from the previous estimate that
\[
\int_{\om}|\varphi|^4(r\om)d\om\leq C_\ga,\quad \forall \f12 \leq r\leq 1
\]
for some constant $C_\ga$ depending only on $\ga$. Here we used the bound of $\varphi$ on the sphere with radius $\f12$.

For $\pa_\om \varphi$, we are not able to derive the same bound due to the lack of good control on $\pa_\om \pa_\om \varphi$. By using Sobolev embedding, we first have
\[
\left(\int_{\frac{1}{3}}^{\frac{1}{2}}\int_{\om}|\pa_\om \varphi|^4 r^2 dr d\om\right)^\f12 \leq C\int_{\frac{1}{3}}^{\f12}\int_{\om} (|\pa_\om \varphi|^2+|\nabla\pa_\om \varphi|^2 r^2) drd\om\leq C
\]
for some universal constant $C$. In particular we can choose some $r_0\in[\frac{1}{3}, \frac{1}{2}]$ such that
\[
\int_{\om}|\pa_\om\varphi|^4(r_0\om)d\om \leq C.
\]
Then replace $\varphi$ with $\pa_\om \varphi$ in the previous $L^4$ estimate. We obtain that
\begin{align*}
\int_{\om}|\pa_\om\varphi|^4d\om &\leq C+\sup\limits_{r_0\leq s\leq r}\left(\int_{\om}|\pa_\om\varphi|^4(s\om)d\om\right)^\f12 \left(\int_{r_0}^{r}\int_{\om}(|\pa_\om\varphi|^2+|\nabla \pa_\om \varphi|^2)(1-|x|)^{-\ga} drd\om \right)^\f12\\
&\leq C+\sup\limits_{r_0\leq s\leq r}\left(\int_{\om}|\pa_\om\varphi|^4(s\om)d\om\right)^\f12 (1-r)^{-\ga},
\end{align*}
which implies that
\begin{align*}
\int_{\om}|\pa_\om \varphi|^4d\om \leq C (1-r)^{-2\ga},\quad \forall r_0\leq r\leq 1.
\end{align*}
Now for any $2\leq p_1\leq 4$, integrate from the sphere with the same $r_0$ such that the $L^4$ norm of $\pa_\om\varphi $ is bounded. We can show that
\begin{align*}
\int_{\om}|\pa_\om\varphi|^{p_1}(r\om) d\om\leq \int_{\om}|\pa_\om \varphi|^{p_1}(r_0\om) d\om+&C_{p_1}\left(\int_{r_0}^r \int_{\om}|\pa_r\pa_\om \varphi|^2(1-|x|)^{\ga} drd\om\right)^\f12\\
\cdot &\left(\int_{r_0}^{r}\int_{\om }|\pa_\om \varphi|^{2p_1-2}(1-|x|)^{-\ga}   drd\om \right)^\f12.
\end{align*}
Now choosing $p_1=2+\ep$, $\ep=\frac{1-\ga}{3}$ and interpolation between $L^2$ and $L^4$, we derive that
\begin{align*}
\int_{\om}|\pa_\om\varphi|^{p_1}(r\om) d\om &\leq C+C\left(\int_{r_0}^{r}\int_{\om }|\pa_\om \varphi|^{2}(1-|x|)^{\ep-1}   drd\om\right)^{\frac{1-\ep}{2}}\left(\int_{r_0}^{r}\int_{\om }|\pa_\om \varphi|^{4}(1-|x|)^{1+\ep}   drd\om\right)^{\frac{\ep}{2}}\\
&\leq C+C_\ga \left(\int_{r_0}^{r}(1-s)^{1+\ep-2\ga}ds\right)^{\f12\ep}\leq C_\ga,\quad \forall r_0\leq r\leq 1.
\end{align*}
Here we have used the bound of the $L^4$ estimate for $\pa_\om \varphi$ and the assumption that $\ga<1$. Since $p_1=2+\ep>2$, this estimate together with the above uniform $L^4$ bound for $\varphi$ leads to the uniform pointwise estimate for $\varphi$.
\end{proof}

We now use these two lemmas to prove Proposition \ref{prop:lineardecay:comp}. Since the equation is linear, without loss of generality we may assume that the initial data verify
\[
\int_{|x|\leq R}(|w_1|^2  +\sum\limits_{i,j}|\Om_{ij}w_1|^2)(R-|x|)^\ga dx=1.
\]
 Recall the notation that $\B_{q}(\tilde{r})$ denotes the spatial ball centered at the point $q$ with radius $\tilde{r}$. By using the energy estimates for the linear wave equation, we derive that
\begin{align*}
\int_{\B_{(t, 0)}(r)}|\pa w|^2 dx+\int_{r}^{t+r}\int_{\om}|(\pa_u, \nabb) w|^2(t+r-s, r\om) s^2 ds d \om\leq \int_{\B_{(0, 0)}(t+r)}|\pa w|^2 dx .
\end{align*}
Here $\pa$ is short for the full derivative $(\pa_t, \nabla)$.
Multiply both side by $(R-t-r)^{\ga-1}$ and integrate $r$ from $0$ to $R-t$. We derive that
\begin{align*}
&\int_{B_{(t, 0)}(R-t)}|\pa w|^2(R-t-|x|)^\ga dx\\
&=\ga \int_{0}^{R-t}(R-t-r)^{\ga-1}\int_{B_{(t, 0)}(r)}|\pa w|^2 dx dr\\
&\leq \ga\int_{0}^{R-t}(R-t-r)^{\ga-1} \int_{s\leq r+t}\int_{\om }|w_1|^2 s^2 dsd\om \\
&=\int_{B_{R}}|w_1|^2(R-\max\{|x|, t\})^{\ga}dx\leq 1.
\end{align*}
By commuting the linear wave equation with the angular momentum vector fields $\Om_{ij}$, we also derive the bound
\[
\int_{\B_{(t, 0)}(R-t)}|\pa\Om_{ij} w|^2(R-t-|x|)^\ga dx\leq 1.
\]
Then by commuting the linear wave equation with $\nabla$, similarly we can derive that
\begin{align*}
&\int_{\B_{(t, 0)}(R-t)}|\pa \nabla w |^2(R-t-|x|)^{2+\ga}dx \leq \int_{\B_{R}}|\nabla w_1|^2(R-\max\{|x|, t\})^{\ga+2}dx\leq 1.
\end{align*}
To apply Lemma \ref{lem:weightedSob:com} to conclude the pointwise bound for $w(t, x)$, it remains to derive the weighted $L^2$ estimate for the solution $w$. For this we rely on the energy flux through the cone $\mathcal{N}^{-}(t+r, 0)$. Indeed since $w$ vanishes on the initial hypersurface, we can show that
\begin{align*}
\int_{\om}|w|^2(t, r\om) d\om &\leq \int_{\om}\left(\int_{r}^{r+t}|\pa_u w|(t+r-s, s\om)ds\right)^2d\om\\
&\leq \int_{\om} \int_{r}^{r+t}|\pa_u w|^2(t+r-s, s\om)s^2 \cdot r^{-1} ds  d\om\\
&\leq r^{-1} \int_{\B_{(0, 0)}(t+r)}|w_1|^2dx.
\end{align*}
Therefore multiplying both side with $(R-t-r)^{\ga-1}$, we derive that
\begin{align*}
\int_{0}^{R-t}\int_{\om}|w|^2(t, r\om)(R-t-r)^{\ga-1} rd\om dr \leq \int_{0}^{R-t}(R-t-r)^{\ga-1}\int_{|x|\leq r+t}|w_1|^2 dx dr\leq \ga^{-1}.
\end{align*}
In particular we have
\begin{align*}
\int_{\B_{(t, 0)}(R-t)}|w|^2(R-t-r)^{\ga}dx\leq (R-t)^2 \int_{0}^{R-t}\int_{\om}|w|^2(R-t-r)^{\ga-1} rd\om dr\leq \ga^{-1}(R-t)^2.
\end{align*}
The linear decay estimate \eqref{eq:lineardecay:comp} then follows by using Lemma \ref{lem:weightedSob:com}. This proves Proposition \ref{prop:lineardecay:comp}.
\fi

Based on the above uniform weighted flux bound, we are able to derive the pointwise estimates for solution of \eqref{eq:NLW:3D:conf}. Let $\phi^{lin}$ be the linear evolution on $\cJ^{+}(\B_R)$ with initial data $(\phi_0, \phi_1)$ on $\B_R$, that is,
\begin{align*}
  \Box\phi^{lin}=0,\quad \phi^{lin}(0, x)=\phi_0,\quad \pa_t\phi^{lin}(0, x)=\phi_1,\quad t+|x|\leq R.
\end{align*}
To avoid too many constants, we make a convention in this section that $A\les B$ means that there is a constant $C$, depending only on $R$, $p$, $\ga$, $\tilde{\cE}_{0, \ga}+\mathcal{I}$ and a fixed small constant $0<\ep<10^{-2}(1-\ga)$, such that $A\leq CB$.

Before stating the main result of this section, we prove two integration lemmas.
\begin{Lem}
\label{lem:bd:vga}
Fix $(t_0, x_0)\in \mathcal{J}^+(B_R)$. For all $0\leq \tilde{r}\leq t_0$, $t=t_0-\tilde{r}$ and $r=|x_0+\tilde{r}\tilde{\om}|$, we have the following uniform bound
\begin{equation*}
\int_{\S_{(t, x_0)}(\tilde{r})}(R-t-r)^{-\ga'}d\tilde{\om}\leq C (R-t)^{\ga'}(R-t_0)^{-\ga'} (v_0+\tilde{r})^{-\ga'}
\end{equation*}
for all $ \ga'<1$ and for some constant $C$ depending only on $\ga'$ and $R$.
\end{Lem}
\begin{proof}
Denote $u_0=R-t_0+r_0$ and $v_0=R-t_0-r_0$ where $r_0=|x_0|$. By the assumption that $t_0+r_0\leq R$, we in particular have that $r\leq R-t$, which implies that
\[
(R-t-r)^{-\ga'}\leq 2 \left((R-t)^2-r^2\right)^{-\ga'}(R-t)^{\ga'}.
\]
Note that $r^2=r_0^2+\tilde{r}^2+2 \tilde{r} x_0\cdot \tilde{\om}$. We can compute that
\begin{align*}
\int_{|\tilde{\om}|=1}(R-t-r)^{-\ga'}d\tilde{\om}&\leq 4\pi (R-t)^{\ga'} \int_{-1}^{1}((R-t)^2-t_0^2-\tilde{r}^2-2\tilde{r}r_0\tau)^{-\ga'}d\tau\\
&\leq 4\pi(1-\ga')^{-1} (R-t)^{\ga'} (r_0 \tilde{r})^{-1} (u_0^{1-\ga'}(v_0+2\tilde{r})^{1-\ga'}-v_0^{1-\ga'}(u_0+2\tilde{r})^{1-\ga'}).
\end{align*}
By definition, we see that
\[
u_0(v_0+2\tilde{r})-v_0(u_0+2\tilde{r})=4\tilde{r}r_0.
\]
As $\gamma'<1$, we derive that
\[
u_0^{1-\ga'} (v_0+2\tilde{r})^{1-\ga'}-v_0^{1-\ga'}(u_0+2\tilde{r})^{1-\ga'}\leq C(R, \ga') \tilde{r}r_0 u_0^{-\gamma'}(v_0+\tilde{r})^{-\gamma'}
\]
for some constant $C(R, \ga')$ depending only on $\ga'$ and $R$.
The lemma then follows as $0\leq r_0\leq R-t_0$.
\end{proof}

The above integration lemma will be used for the larger $p$ case when $p>\frac{1+\sqrt{17}}{2}$. The following specific lemma will be used for small $p$.
\begin{Lem}
\label{lem:bd:vga:smallp:in}
Fix $(t_0, x_0)\in \mathcal{J}^+(B_R)$. For all $0\leq \tilde{r}\leq t_0$, $t=t_0-\tilde{r}$, $r=|x_0+\tilde{r}\tilde{\om}|$ and $0<\ga<1$, $0\leq \a< 1$, we have the following uniform bound
\begin{equation*}
\int_{\S_{(t, x_0)}(\tilde{r})} ((1-\tau)u_*^{\ga}+v_*^{\ga})^{-\a} d\tilde{\om}\leq C (R-t_0)^{-\a \ga}
\end{equation*}
for some constant $C$ depending only on $\ga$ , $\a$ and $R$.
\end{Lem}
\begin{proof}
Using the same notations from the previous Lemma, denote $s=-\om_0\cdot \tilde{\om}$.
Recall that
  \begin{align*}
    &r^2=(\tilde{r}-r_0s)^2+(1-s^2)r_0^2,\quad \tau r =(\tilde{x}+x_0)\cdot \tilde{\om}=\tilde{r}-r_0s.
  \end{align*}
We can write the integral as
\begin{align*}
  \int_{\S_{(t, x_0)}(\tilde{r})} ((1-\tau)u_*^{\ga}+v_*^{\ga})^{-\a} d\tilde{\om}=2\pi\int_{-1}^{1} ((1-r^{-1}(\tilde{r}-r_0 s))u_*^{\ga}+v_*^{\ga})^{-\a}ds.
\end{align*}
Note that $R-t\leq u_*\leq 2(R-t)$. Hence
\begin{align*}
  \int_{\S_{(t, x_0)}(\tilde{r})} ((1-\tau)u_*^{\ga}+v_*^{\ga})^{-\a} d\tilde{\om}\les \int_{-1}^{1} ((1-r^{-1}(\tilde{r}-r_0 s))(R-t)^{\ga}+v_*^{\ga})^{-\a} ds.
\end{align*}
Here and in the following of the proof the implicit constants rely only on $R$, $\a$, $\ga$.

For the case when $r_0\leq \frac{3}{4} (R-t_0)$, it holds that $$v_*\geq R-t_0-r_0\geq \frac{1}{4} (R-t_0),$$
from which we conclude that
\begin{align*}
  \int_{\S_{(t, x_0)}(\tilde{r})} ((1-\tau)u_*^{\ga}+v_*^{\ga})^{-\a} d\tilde{\om}\les (R-t_0)^{-\a\ga}.
\end{align*}
Otherwise for the case when $r_0\geq \frac{3}{4}(R-t_0)$, and if $\tilde{r}\leq 2 r_0 $, note that
\begin{align*}
  1-\tau=1-r^{-1}(\tilde{r}-r_0s) \geq \frac{1-s^2}{100}.
\end{align*}
Therefore we can show that
\begin{align*}
  &\int_{-1}^{1}((1-\tau)(R-t)^{\ga}+v_*^{\ga})^{-\a} ds
  \les \int_{-1}^{1} (R-t)^{-\a\ga}(1-s^2)^{-\a} ds \les (R-t_0)^{-\a\ga}.
\end{align*}
For the remaining case $\tilde{r}\geq 2r_0$, we instead have
\begin{align*}
   v_*=R-t-r\geq v_0+10^{-2}(1+s) r_0
\end{align*}
Therefore we derive that
\begin{align*}
  \int_{-1}^{1}((1-\tau)(R-t)^{\ga}+v_*^{\ga})^{-\a} ds
  &\les \int_{-1}^{1} (v_0+(1+s) r_0)^{-\a\ga}  ds\\
  &\les r_0^{-1}((v_0+2r_0)^{1-\a\ga}-v_0^{1-\a\ga})\\
  &\les (R-t_0)^{-\a\ga}
\end{align*}
as $0\leq \a\ga<1$. This proves the lemma.
\end{proof}

\iffalse
\begin{Lem}
\label{lem:bd:vga:smallp:Grownall}
Fix $(t_0, x_0)\in \mathcal{J}^+(B_R)$. For all $0\leq \tilde{r}\leq \f12 (R-t_0)$, $t=t_0-\tilde{r}$ and $r=|x_0+\tilde{r}\tilde{\om}|$, we have the following uniform bound
\begin{equation*}
\int_{\S_{(t, x_0)}(\tilde{r})} ((1-\tau)u_*^{\ga}+v_*^{\ga})^{-\frac{p-1}{2}}v_*^{p-3}d\tilde{\om}\leq C (R-t_0)^{-\frac{p-1}{2}\ga}(R-t_0-r_0)^{p-3}
\end{equation*}
for some constant $C$ depending only on $p$, $\ga$  and $R$.
\end{Lem}
\begin{proof}
Denote $u_0=R-t_0+r_0$ and $v_0=R-t_0-r_0$ where $r_0=|x_0|$.
Denote $s=-\om_0\cdot \tilde{\om}$ with $\om_0=r_0^{-1}x_0$. Note that
  \begin{align*}
    &r^2=|x_0+\tilde{x}|^2=\tilde{r}^2+r_0^2+2r_0\tilde{r}\om_0\cdot \tilde{\om}=(\tilde{r}-r_0s)^2+(1-s^2)r_0^2,\\
    & \tau r=r\om\cdot \tilde{\om}=(\tilde{x}+x_0)\cdot \tilde{\om}=\tilde{r}-r_0s.
  \end{align*}
We can write the integral as
\begin{align*}
  \int_{\S_{(t, x_0)}(\tilde{r})} ((1-\tau)u_*^{\ga}+v_*^{\ga})^{-\frac{p-1}{2}}v_*^{p-3}d\tilde{\om}=2\pi\int_{-1}^{1} ((1-r^{-1}(\tilde{r}-r_0 s))u_*^{\ga}+v_*^{\ga})^{-\frac{p-1}{2}}v_*^{p-3}ds
\end{align*}
Since $R-t\leq u_*\leq 2(R-t)$, hence
\begin{align*}
  \int_{\S_{(t, x_0)}(\tilde{r})} ((1-\tau)u_*^{\ga}+v_*^{\ga})^{-\frac{p-1}{2}}v_*^{p-3}d\tilde{\om}\les \int_{-1}^{1} ((1-r^{-1}(\tilde{r}-r_0 s))(R-t)^{\ga}+v_*^{\ga})^{-\frac{p-1}{2}}v_*^{p-3}ds.
\end{align*}
For the case when $r_0\leq \frac{3}{4} (R-t_0)$, it holds that $v_*\geq R-t_0-r_0\geq \frac{1}{4} (R-t_0)$. In particular, we conclude that
\begin{align*}
  \int_{\S_{(t, x_0)}(\tilde{r})} ((1-\tau)u_*^{\ga}+v_*^{\ga})^{-\frac{p-1}{2}}v_*^{p-3}d\tilde{\om}\les (R-t_0)^{p-3-\frac{p-1}{2}\ga}.
\end{align*}
Here we may notice that $1<p<3$, $\ga>0$. When $ \frac{3}{2}\tilde{r}\leq \frac{3}{4} (R-t_0)\leq r_0 $ and $s\geq r_0^{-1}\tilde{r}$, then
\begin{align*}
&\int_{r_0^{-1}\tilde{r}}^{1} ((1-r^{-1}(\tilde{r}-r_0 s))(R-t)^{\ga}+v_*^{\ga})^{-\frac{p-1}{2}}v_*^{p-3}ds \\
&\les \int_{r_0^{-1}\tilde{r}}^{1} (R-t)^{-\frac{p-1}{2}\ga} v_*^{p-3}ds\\
&\les (R-t)^{-\frac{p-1}{2}\ga+3-p} \int_{r_0^{-1}\tilde{r}}^{1} ((R-t)^2-r_0^2-\tilde{r}^2+2r_0\tilde{r}s)^{p-3}ds\\
&\les (R-t)^{-\frac{p-1}{2}\ga+3-p} (r_0\tilde{r})^{-1}\left(((R-t)^2-(r_0-\tilde{r})^2)^{p-2}-((R-t)^2-r_0^2+\tilde{r}^2)^{p-2}\right)\\
&\les (R-t)^{-\frac{p-1}{2}\ga+3-p} (r_0\tilde{r})^{-1}(r_0\tilde{r}-\tilde{r}^2)(u_0v_0+2u_0\tilde{r})^{p-3}\\
&\les (R-t)^{-\frac{p-1}{2}\ga+3-p} u_0^{p-3} (v_0+2\tilde{r})^{p-3}\\
&\les (R-t_0)^{-\frac{p-1}{2}\ga} (v_0+\tilde{r})^{p-3}.
\end{align*}
On the interval $[-1, r_0^{-1}\tilde{r}]$, first we have
\begin{align*}
  1-r^{-1}(\tilde{r}-r_0s)=\frac{(1-s^2)r_0^2}{r(r+\tilde{r}-r_0s)}\geq \frac{(1+s)r_0^2}{2(\tilde{r}-r_0s)^2+2(1-s^2)r_0^2}\geq \frac{1+s}{100}
\end{align*}
\begin{align*}
  R-t-r=v_0+\tilde{r}+r_0-\sqrt{\tilde{r}^2+r_0^2-2\tilde{r}r_0s}=v_0+\frac{2r_0\tilde{r}(1+s)}{\tilde{r}+r_0+r}\geq v_0+10^{-2}(1+s)\tilde{r}.
\end{align*}
Therefore we can show that
\begin{align*}
&\int^{r_0^{-1}\tilde{r}}_{-1} ((1-r^{-1}(\tilde{r}-r_0 s))(R-t)^{\ga}+v_*^{\ga})^{-\frac{p-1}{2}}v_*^{p-3}ds \\
&\les \int^{r_0^{-1}\tilde{r}}_{-1} ((R-t)^{\ga}(1+s)+(v_0+(1+s)\tilde{r})^{\ga})^{-\frac{p-1}{2}} (v_0+(1+s)\tilde{r})^{p-3}ds\\
&\les (R-t_0)^{-\frac{p-1}{2}\ga} v_0^{p-3}.
\end{align*}
\end{proof}

\begin{Lem}
\label{lem:bd:vga:smallp:direct}
Fix $(t_0, x_0)\in \mathcal{J}^+(B_R)$. For all $0\leq \tilde{r}\leq t_0$, $t=t_0-\tilde{r}$ and $r=|x_0+\tilde{r}\tilde{\om}|$, we have the following uniform bound
\begin{equation*}
\int_{\S_{(t, x_0)}(\tilde{r})} ((1-\tau)u_*^{\ga}+v_*^{\ga})^{-p}v_*^{p-3} d\tilde{\om}\leq C (R-t_0-r_0)^{p-3-(p-1)\ga}(R-t_0)^{-\ga}
\end{equation*}
for all $ \ga'<1$ for some constant $C$ depending only on $\ga'$ and $R$.
\end{Lem}
\begin{proof}
Again, we can write the integral as
\begin{align*}
  \int_{\S_{(t, x_0)}(\tilde{r})} ((1-\tau)u_*^{\ga}+v_*^{\ga})^{-p}v_*^{p-3} d\tilde{\om}\les \int_{-1}^{1}((1-\tau)(R-t)^{\ga}+v_*^{\ga})^{-p}v_*^{p-3}ds.
\end{align*}
Denote $u_0=R-t_0+r_0$ and $v_0=R-t_0-r_0$ where $r_0=|x_0|$.
For the case when $\tilde{r}\leq 2r_0$, the argument is quite similar to the proof of previous lemma. One can first show that
\begin{align*}
  1-\tau=1-r^{-1}(\tilde{r}-r_0s) \geq \frac{1-s^2}{100}, \quad v_*=R-t-r\geq v_0+10^{-2}(1+s)\tilde{r}.
\end{align*}
Hence
\begin{align*}
  &\int_{-1}^{1}((1-\tau)(R-t)^{\ga}+v_*^{\ga})^{-p}v_*^{p-3}ds\\
  &\les \int_{-1}^{1} ((R-t)^{\ga}(1-s^2)+(v_0+(1+s)\tilde{r})^{\ga})^{-p}(v_0+(1+s)\tilde{r})^{p-3}ds\\
  &\les (R-t)^{-\ga}v_0^{\ga(1-p)+p-3}.
\end{align*}
Now for the case when $\tilde{r}\geq 2r_0$, we instead have
\begin{align*}
   1-\tau=1-r^{-1}(\tilde{r}-r_0s) \geq \frac{1-s^2}{100}r_0^2\tilde{r}^{-2}, \quad v_*=R-t-r\geq v_0+10^{-2}(1+s) r_0.
\end{align*}
Thus
Hence
\begin{align*}
  &\int_{-1}^{1}((1-\tau)(R-t)^{\ga}+v_*^{\ga})^{-p}v_*^{p-3}ds\\
  &\les \int_{-1}^{1} ((R-t)^{\ga}(1-s^2)r_0^{2}\tilde{r}^{-2}+(v_0+(1+s) r_0)^{\ga})^{-p}(v_0+(1+s) r_0)^{p-3}ds\\
  &\les \int_{-1}^{1} (v_0+(1+s) r_0)^{-\ga}  v_0^{-\ga(p-1)+p-3} ds\\
  &\les v_0^{-\ga(p-1)+p-3}(R-t_0)^{-\ga}.
\end{align*}
\end{proof}
\fi

Define
\[
\a_p=\frac{3+(p-2)^2}{(p+1)(5-p)},\quad 1<p<5.
\]
Now we are ready to prove the main result of this section.
\begin{Prop}
\label{Prop:pointwise:decay:EM}
The solution $\phi$ to the equation \eqref{eq:NLW:3D:conf} on $\mathcal{J}^{+}(\B_{R})$ verifies the following pointwise decay estimates:
\begin{itemize}
\item If
\begin{equation*}
 \frac{1+\sqrt{17}}{2}<p<5,\quad  0<\ga<1,\quad (p-1)(3-\ga)>4,
\end{equation*}
then
\begin{equation}
\label{eq:phi:pt:Br:EM:largep}
|\phi(t_0, x_0)|\leq C \sup\limits_{|x|\leq R-t_0}|\phi^{lin}(t_0, x)|.
\end{equation}
\item
For the case when
\begin{align*}
  2<p\leq \frac{1+\sqrt{17}}{2},\quad \ga=2-\ga_0+\ep,\quad  1<\ga_0-\ep<p-1,
\end{align*}
for all $\b\leq \frac{p-1}{p+1}\ga_0-\ep$, we have
\begin{equation}
  \label{eq:phi:pt:Br:EM:smallp}
  |\phi(t_0, x_0)|\leq C (1+\sup\limits_{t+|x|\leq R}|\phi^{lin}u_*^{1-\b} v_*^{1-\a_p \ga_0}|)u_0^{-1+\b} v_0^{-1+\a_p\ga_0}.
\end{equation}
The constant $C$ depends only on $\mathcal{I}+\tilde{\cE}_{0, \ga}$, $R$, $p$, $\ga$, $\b$ and $\ep>0$. Here $u_*=R-t+r$, $v_*=R-t-r$ and
$u_0=R-t_0$, $v_0=R-t_0-|x_0|$. The small positive constant $\ep$ may be different in different places.
\end{itemize}
\end{Prop}
\begin{proof}
To avoid too many constants, the implicit constant in $\les$ in the following proof depends only on $\mathcal{I}+\tilde{\cE}_{0, \ga}$, $R$, $p$, $\ga$, $\b$ and $\ep>0$.

The proof for this Proposition relies on the representation formula for linear wave equation. The nonlinearity will be controlled by using the weighted flux bound in Proposition \ref{prop:EF:cone:gamma}. Recall that for any $q=(t_0, x_0)\in \mathcal{J}^{+}(\B_{R})$, we have the representation formula for the solution
\begin{equation}
\label{eq:rep4phi:comp}
\begin{split}
4\pi\phi(t_0, x_0)&
4\pi\phi^{lin}(t_0, x_0)
-\int_{\mathcal{N}^{-}(q)}\Box \phi \ \tilde{r} d\tilde{r}d\tilde{\om}.
\end{split}
\end{equation}
\iffalse
For the linear evolution part,
by the definition of $\tilde{\cE}_{1, \ga}$ and Lemma \ref{lem:weightedSob:com}, we derive that
\[
|\phi_0|\les \sqrt{\tilde{\cE}_{1, \ga}},\quad \forall |x|\leq R.
\]
In particular we have
\begin{align*}
  \int_{\tilde{\om}}|\phi(0, x_0+t_0\tilde{\om})|  d\tilde{\om}\les \sqrt{\tilde{\cE}_{1, \ga}}\les (R-t_0)^{-\frac{1+\ga}{2}}\sqrt{\tilde{\cE}_{1, \ga}}.
\end{align*}
For the estimate for $\pa\phi(0, x_0+t_0\tilde{\om})$, consider the linear wave equation
\[
\Box w=0,\quad w(0, x)=0, \quad \pa_t w(0, x)=|\pa\phi|(0, x)=|\phi_1|+|\nabla\phi_0|, \quad |x|\leq R.
\]
Then by the representation formula for linear wave equations we have
\begin{align*}
4\pi w(t_0, x_0)=t_0\int_{\tilde{\om}}\pa_t w(0, x_0+t_0\tilde{\om})d\om=t_0\int_{\tilde{\om}}|\pa\phi|(0, x_0+t_0\tilde{\om})d\tilde{\om}.
\end{align*}
By definition, we estimate that
\begin{align*}
\mathcal{E}_{\ga}[\pa_t w(0, x)](\B_R)&=\int_{\B_{R}}(|\pa \phi|^2+|\Om_{ij}|\pa\phi||^2)(R-|x|)^\ga+|\nabla |\pa\phi||^2 (R-|x|)^{2+\ga}dx\\
&\les  \int_{\B_{R}}(|\pa \phi|^2+|{\Om_{ij}}\pa\phi|^2)(R-|x|)^\ga+|\pa\pa\phi|^2 (R-|x|)^{2+\ga}dx\\
&\les \tilde{\cE}_{1, \ga}.
\end{align*}
 Thus by using Proposition \ref{prop:lineardecay:comp}, we conclude that
 \begin{align*}
 | w(t_0, x_0)|\les  (R-t_0)^{-\frac{1+\ga}{2}}\sqrt{\tilde{\cE}_{1, \ga}}
 \end{align*}
Therefore the linear evolution verifies the bound
\begin{align*}
  |\int_{\tilde{\om}}t_0  \phi_1(x_0+t_0\tilde{\om})d\tilde{\om}+\pa_{t_0}\big(\int_{\tilde{\om}}t_0  \phi_0(x_0+t_0\tilde{\om})d\tilde{\om}   \big)|
  &\les \int_{\tilde{\om}} t_0 |\pa\phi| +|\phi_0| d\tilde{\om} \les (R-t_0)^{-\frac{1+\ga}{2}}\sqrt{\tilde{\cE}_{1, \ga}}.
\end{align*}
\fi
We mainly need to control the nonlinear part. From the equation \eqref{eq:NLW:3D:conf} as well as the flux bound \eqref{eq:comp:v:EF}, we can estimate that
\begin{align*}
  &|\int_{\mathcal{N}^{-}(q)}\Box \phi \ \tilde{r} d\tilde{r}d\tilde{\om}|\\
  &\leq \int_{\mathcal{N}^{-}(q)}\La^{3-p}|\phi|^{p}\ \tilde{r} d\tilde{r}d\tilde{\om}\\
  &\leq \left(\int_{\mathcal{N}^{-}(q)} v_*^{\ga}\La^{3-p}|\phi|^{p+1}\ \tilde{r}^{2} d\tilde{r}d\tilde{\om}\right)^{\frac{p-1}{p+1}}\left(\int_{\mathcal{N}^{-}(q)}v_*^{-\frac{p-1}{2}\ga}\La^{3-p}|\phi|^{\frac{p+1}{2}}\ \tilde{r}^{\frac{3-p}{2}} d\tilde{r}d\tilde{\om}\right)^{\frac{2}{p+1}}\\

  &\les \left(\int_{0}^{t_0}(\sup\limits_{x}|(R-t)^{\frac{1+\ga}{2}}\phi|^{\frac{p+1}{2}})\int_{\S_{(t_0-\tilde{r}, x_0)(\tilde{r})}}(R-t)^{-\frac{(p+1)(\ga+1)}{2}} \La^{3-p}v_*^{-\frac{p-1}{2}\ga}\tilde{r}^{\frac{3-p}{2}}d\tilde{\om} d\tilde{r} \right)^{\frac{2}{p+1}}.
\end{align*}
Here under the coordinates $(\tilde{t}, \tilde{x})$, $t=t_0-\tilde{r}$.

 Let's first consider the decay estimate of \eqref{eq:phi:pt:Br:EM:largep} for the larger $p$ case. By our assumption, we in particular have that
\[
p-3-\frac{p-1}{2}\ga=\frac{(p-1)(2-\ga)}{2}-2>\frac{(p-1)(3-\ga)}{4}-2>-1.
\]
Since $v_*=R-t-r$, thus by using Lemma \ref{lem:bd:vga}, we can bound that
\begin{align*}
  \int_{\S_{(t_0-\tilde{r}, x_0)(\tilde{r})}} v_*^{p-3-\frac{p-1}{2}\ga} d\tilde{\om}\les \left((R-t)^{-1}(R-t_0)\tilde{r}\right)^{p-3-\frac{p-1}{2}\ga}.
\end{align*}
Define
\begin{align*}
\tilde{\cM}(t)&=\sup\limits_{|x|\leq R-t}|(R-t)^{\frac{1+\ga}{2}}\phi|^{\frac{p+1}{2}},  \quad \forall 0\leq t\leq R\\
\tilde{f}(t_0, \tilde{r})&= (R-t_0)^{\frac{(p+1)(\ga+1)}{2}}(R-t_0+\tilde{r})^{\frac{p-1}{2}\ga-\frac{(p+1)(\ga+1)}{2}}(R-t_0)^{p-3-\frac{(p-1)\ga}{2}}\tilde{r}^{\frac{p-3}{2}-\frac{(p-1)\ga}{2}} .
\end{align*}
 We then conclude that
\begin{align*}
 |\phi(t_0, x_0)|^{\frac{p+1}{2}}&\les |\phi^{lin}(t_0, x_0)|^{\frac{p+1}{2}}+ \int_{0}^{t_0}\tilde{\cM}(t_0-\tilde{r} ) (R-t_0)^{-\frac{1+\ga}{2}\cdot \frac{p+1}{2}}\tilde{f}(t_0,\tilde{r}) d\tilde{r},
\end{align*}
which implies that
\begin{align}
\label{eq:cM:0}
\tilde{\cM}(t_0)&\les \sup\limits_{|x_0|\leq R-t_0}|\phi^{lin}(t_0, x_0)(R-t_0)^{\frac{1+\ga}{2}}|^{\frac{p+1}{2}} + \int_{0}^{t_0}\tilde{\cM}(t_0-\tilde{r} )  \tilde{f}(t_0,\tilde{r}) d\tilde{r}.
\end{align}
To apply Gronwall's inequality, we check that $\tilde{f}(t_0, \tilde{r})$ is uniformly integrable on $[0, t_0]$ with respect to $\tilde{r}$. Indeed, for the case when $p\geq 3$, we show that
\begin{align*}
  \int_0^{t_0}\tilde{f}(t_0, \tilde{r})d\tilde{r}\leq \int_0^{t_0}(R-t_0)^{p-3}\tilde{r}^{\frac{p-3}{2}-\frac{(p-1)\ga}{2}}d\tilde{r}\les (R-t_0)^{p-3}t_0^{\frac{(p-1)(1-\ga)}{2}}\les 1.
\end{align*}
Here the implicit constant relies only on $p$, $\ga$ and $R$.

For the case when $p<3$, split the integral into two parts. For small $\tilde{r}$, we have similar bound
\begin{align*}
  \int_0^{\min\{R-t_0, t_0\}}\tilde{f}(t_0, \tilde{r})d\tilde{r}\leq \int_0^{R-t_0}(R-t_0)^{p-3}\tilde{r}^{\frac{p-3}{2}-\frac{(p-1)\ga}{2}}d\tilde{r}\les  (R-t_0)^{p-3+\frac{(p-1)(1-\ga)}{2}} \les 1
\end{align*}
due to the relation
\[
p-3+\frac{(p-1)(1-\ga)}{2}=\frac{(p-1)(3-\ga)}{2}-2>0.
\]
For large $\tilde{r}\geq \min\{R-t_0, t_0\}$, note that for this remaining case ($p<3$),
\begin{align*}
  \frac{p-3}{2}-\frac{(p-1)\ga}{2}<0,\quad \frac{p-1}{2}\ga-\frac{(p+1)(\ga+1)}{2}<-1
\end{align*}
as $p>1$, $\ga>0$. We thus can bound that
\begin{align*}
  \int_{\min\{R-t_0, t_0\}}^{t_0}\tilde{f}(t_0, \tilde{r}) d\tilde{r}&\leq \int_{\min\{R-t_0, t_0\}}^{t_0} (R-t_0)^{\frac{(p+1)(\ga+1)}{2}+ \frac{3(p-3)}{2}-(p-1)\ga}(R-t_0+\tilde{r})^{\frac{p-1}{2}\ga-\frac{(p+1)(\ga+1)}{2}} d\tilde{r}\\
  &\les   (R-t_0)^{\frac{(p+1)(\ga+1)}{2}+ \frac{3(p-3)}{2}-(p-1)\ga+\frac{p-1}{2}\ga-\frac{(p+1)(\ga+1)}{2}+1} \\
  &=(R-t_0)^{\frac{(p-1)(3-\ga)-4}{2}}\les 1.
\end{align*}
In view of \eqref{eq:cM:0}, Gronwall's inequality then implies that
\begin{align*}
  |\phi(t_0, x_0)|\les \sup\limits_{|x|\leq R-t_0} |\phi^{lin}(t_0, x)|.
\end{align*}
This shows the bound \eqref{eq:phi:pt:Br:EM:largep}.

\bigskip

Next for estimate \eqref{eq:phi:pt:Br:EM:smallp} with small power $2<p\leq \frac{1+\sqrt{17}}{2}$, we control the nonlinearity directly by using the weighted flux for large $\tilde{r}$. We split the integral into two parts: specifically for the smaller $\tilde{r}$ on $[0, t_*]$ and larger $\tilde{r}$ on $[t_*, t_0]$, where we define
\[
t_*= u_0^{\frac{2(3-p)+(p-1)\ga}{5-p} } v_0^{\frac{2(3-p)}{5-p}},\quad u_0=R-t_0.\quad v_0=R-t_0-r_0.
\]
Without loss of generality, we may assume that $t_*<t_0$. Otherwise it suffices to evaluate the integral on the single interval $[0, t_0]$.
For the integral on $[t_*, t_0]$, from the weighted energy estimate \eqref{eq:comp:v:EF}, we can show that
\begin{align*}
  &|\int_{\mathcal{N}^{-}(q)\cap\{\tilde{r}\geq t_*\}}\Box \phi \ \tilde{r} d\tilde{r}d\tilde{\om}|\\

  &\leq \left(\int_{\mathcal{N}^{-}(q)} ((1-\tau)u_*^{\ga}+v_*)^{\ga})\La^{3-p}|\phi|^{p+1}\ \tilde{r}^{2} d\tilde{r}d\tilde{\om}\right)^{\frac{p}{p+1}}\\
  &\quad \cdot \left(\int_{\mathcal{N}^{-}(q)\cap\{\tilde{r}\geq t_*\}}((1-\tau)u_*^{\ga}+v_*)^{\ga})^{-p}\La^{3-p} \tilde{r}^{1-p} d\tilde{r}d\tilde{\om}\right)^{\frac{1}{p+1}}\\
  &\les \left(\int_{t_*}^{t_0} \int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})}v_*^{p-3} ((1-\tau)u_*^{\ga}+v_*)^{\ga})^{-p}  u_*^{p-3} \tilde{r}^{1-p} d\tilde{\om} d\tilde{r}\right)^{\frac{1}{p+1}}.
\end{align*}
 Since $v_*\geq  v_0$, by using Lemma \ref{lem:bd:vga:smallp:in}, we then can show that
\begin{align*}
  &|\int_{\mathcal{N}^{-}(q)\cap\{\tilde{r}\geq t_*\}}\Box \phi \ \tilde{r} d\tilde{r}d\tilde{\om}|^{p+1}\\
  &\les v_0^{p-3-(p-1)\ga-\ep}\int_{t_*}^{t_0}(R-t)^{p-3}\tilde{r}^{1-p}  \int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})}((1-\tau)u_*^{\ga}+v_*)^{\ga})^{-1+\ep \ga^{-1}} d\tilde{\om} d\tilde{r}\\
  &\les v_0^{p-3-(p-1)\ga-\ep} (R-t_0)^{-\ga+\ep} \int_{t_*}^{t_0}(R-t_0)^{p-3}\tilde{r}^{1-p}  d\tilde{r}\\
  &\les v_0^{p-3-(p-1)\ga-\ep} (R-t_0)^{p-3-\ga+\ep} t_*^{2-p} \\
  &\les u_0^{-\frac{(3-p)(p+1)+\ga(p^2-4p+7)}{5-p}+\ep}v_0^{-\frac{(3-p)(p+1)}{5-p}-(p-1)\ga-\ep}\\
  &\les u_0^{(-1+\b_1\ga_0+(1-\a_p)\ep-\delta_1)(p+1) } v_0^{(-1+\a_p \ga_0-(1+\b_1)\ep+\delta_1)(p+1)}\\
  &\les u_0^{(-1+\b_1\ga_0-(\a_p+\b_1)\ep)(p+1) } v_0^{(-1+\a_p \ga_0)(p+1)},
\end{align*}
in which
\begin{align*}
  \ga=2-\ga_0+\ep,\quad \a_p=\frac{p^2-4p+7}{(5-p)(p+1)},\quad \b_1=\frac{p-1}{p+1},\quad \delta_1=\frac{2(p-2)(3-p)(\ga_0-1)}{(5-p)(p+1)}.
\end{align*}
Since $\ep$ is arbitrary, without any confusion, replacing $(\a_1+\b_1)\ep$ by $\ep$, we therefore derive that
\begin{align*}
  |\phi(t_0, x_0)|&\les |\phi^{lin}(t_0, x_0)|+|\int_{\mathcal{N}^{-}(q) }\Box \phi \ \tilde{r} d\tilde{r}d\tilde{\om}| \\
  &\les  |\phi^{lin}(t_0, x_0)|+u_0^{-1+\b_1\ga_0-\ep} v_0^{-1+\a_p \ga_0}+|\int_{\mathcal{N}^{-}(q)\cap\{\tilde{r}\leq t_*\}}\Box \phi  \tilde{r} d\tilde{r}d\tilde{\om}|.
\end{align*}
Now for smaller $\tilde{r}$, we rely on Gronwall's inequality and similar to the above argument, we first have
\begin{align*}
  &|\int_{\mathcal{N}^{-}(q)\cap\{\tilde{r}\leq t_*\}}\Box \phi \ \tilde{r} d\tilde{r}d\tilde{\om}| \les \left(\int_{\mathcal{N}^{-}(q)\cap\{\tilde{r}\leq t_*\}}((1-\tau)u_*^{\ga}+v_*^\ga)^{-\frac{p-1}{2}}\La^{3-p}|\phi|^{\frac{p+1}{2}}\ \tilde{r}^{\frac{3-p}{2}} d\tilde{r}d\tilde{\om}\right)^{\frac{2}{p+1}}.
\end{align*}
Now define
\begin{align*}
  \mathcal{M}_2[\phi](t)=\sup\limits_{|x|\leq R-t} |\phi(t, x)u_*^{1-\b}v_*^{ 1-\a_p \ga_0}|^{\frac{p+1}{2}},\quad \b\leq \b_1-\ep.
\end{align*}
As $\ga_0<p-1$, it can be checked that
\begin{align*}
  \b_1\ga_0<1,\quad \a_p\ga_0\leq 1.
\end{align*}
Notice that $v_*\geq v_0$, $u_*\geq u_0$. The previous inequality in particular leads to
\begin{align*}
  \mathcal{M}_2[\phi](t_0)\les 1+\mathcal{M}_2[\phi^{lin}](t_0)+  \int^{t_*}_{0} \mathcal{M}_2[\phi](t_0-\tilde{r})\int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})} ((1-\tau)u_*^{\ga}+v_*^\ga)^{-\frac{p-1}{2}}\La^{3-p}  \tilde{r}^{\frac{3-p}{2}}d\tilde{\om}d\tilde{r} .
\end{align*}
Now by using Lemma \ref{lem:bd:vga:smallp:in} as $p<3$ and the choice of $t_*$, we estimate that
\begin{align*}
  & \int^{t_*}_{0}\int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})}((1-\tau)u_*^{\ga}+v_*^\ga)^{-\frac{p-1}{2}}\La^{3-p} \tilde{r}^{\frac{3-p}{2}} d\tilde{r}d\tilde{\om}\\

  &\les \int^{t_*}_{0}\int_{\S_{(t_0-\tilde{r}, x_0)}(\tilde{r})} ((1-\tau)u_*^{\ga}+v_*^\ga)^{-\frac{p-1}{2}}v_0^{p-3}u_0^{p-3} \tilde{r}^{\frac{3-p}{2}} d\tilde{r}d\tilde{\om}\\
  &\les  (R-t_0)^{-\frac{p-1}{2}\ga}v_0^{p-3}(R-t_0)^{p-3} t_*^{\frac{5-p}{2}} \les 1.
\end{align*}
Hence by using Gronwall's inequality, we conclude that
\begin{align*}
  |\phi(t_0, x_0)|\les (1+\mathcal{M}_2[\phi^{lin}](t_0))^{\frac{2}{p+1}}u_0^{-1+\b}v_0^{-1+\a_p\ga_0}.
\end{align*}
We thus finished the proof for Proposition \ref{Prop:pointwise:decay:EM}.

\end{proof}

\section{The solution in the interior region and proof for the main theorem}
The aim of this section is to apply the result of Proposition \ref{Prop:pointwise:decay:EM} from previous section together with estimates of Propositions \ref{prop:energyflux:H:ex}, \ref{prop:NLW:def:3d:ex:ID:H} to derive the asymptotic decay properties of solutions to the nonlinear wave equations \eqref{eq:NLW:semi:3d} in the interior region $\{t+2\geq |x|\}$ which is contained inside $\mathbf{D}$, inclosed by the forward hyperboloid $\mathbb{H}$ defined in \eqref{eq:def4Hyperboloid}.

Define the conformal map
\begin{align*}
\mathbf{\Phi}:(t, x)\longmapsto (\tilde{t}, \tilde{x})=\left(-\frac{t^*}{(t^*)^2-|x|^2}+R^*,\quad \frac{x}{(t^*)^2-|x|^2}\right)
\end{align*}
from the region $\mathbf{D}$ to Minkowski space.
The image of $\mathbf{\Phi}(\mathbf{D})$ is a truncated backward light cone
\begin{align*}
\mathbf{\Phi}(\mathbf{D})=\left\{(\tilde{t}, \tilde{x})|\quad \tilde{t}+|\tilde{x}|<R^*,\quad \tilde{t}\geq 0\right\}.
\end{align*}
Denote
\begin{equation*}
\Lambda(t, x)=(t^*)^2-|x|^2.
\end{equation*}
Direct computation shows that $\tilde{\phi}=(\Lambda \phi)\circ \mathbf{\Phi}^{-1}$ (as a scalar field in $(\tilde{t}, \tilde{x})$ variables on $\mathbf{\Phi}(\mathbf{D})$) verifies the nonlinear wave equation \eqref{eq:NLW:3D:conf}. For simplicity we may identify $\Lambda \phi$ with $(\Lambda \phi)\circ \mathbf{\Phi}^{-1}$.
The initial hypersurface for the above backward light cone is a ball with radius $R^*$
$$ \mathbf{\Phi}(\mathbb{H}) =\{(0, \tilde{x})||\tilde{x}|\leq R^*\}.$$
By doing this conformal transformation, the Cauchy problem of equation \eqref{eq:NLW:semi:3d} with initial hypersurface $\Hy$ is then equivalent to the Cauchy problem of equation \eqref{eq:NLW:3D:conf} with initial hypersurface $\mathbf{\Phi}(\mathbb{H}) $.

To apply the result of Proposition \ref{Prop:pointwise:decay:EM}, we need first to control the weighted energy norm $\tilde{\mathcal{E}}_{0, \ga}$ and the weighted spacetime integral $\mathcal{I}$ (defined before Proposition \ref{prop:EF:cone:gamma}) in terms of $\mathcal{E}_{0, \ga_0}[\phi]$.
Setting
\begin{equation*}
 \ga=2-\ga_0+\ep
\end{equation*}
with
\begin{equation}
\label{eq:ga1}
0<\ep<10^{-1}\min\{\ga_0-1, 2-\ga_0, |\ga_0+1-\frac{4}{p-1}|\}.
\end{equation}
By our assumption on $\ga_0$, we in particular have $0<\ga<1$.

Let $$\tilde{u}=\frac{\tilde{t}-\tilde{r}}{2},\quad  \tilde{v}=\frac{\tilde{t}+\tilde{r}}{2}, \quad \tilde{\om}=\frac{\tilde{x}}{|\tilde{x}|}$$ be the null coordinates system, with the associated null frame $\{\tilde{L} ,  \tilde{\Lb} ,  \La e_1,  \La e_2\}$ on $\mathbf{\Phi}(\mathbf{D})$. Here we recall that $\{L, \Lb, e_1, e_2 \}$ is the null frame on $\mathbf{D}$ under the coordinates $(t, x)$. Recall the weighted energy norm $\tilde{\cE}_{0, \ga}$ associated to $\tilde{\phi}$ on $\mathbf{\Phi}(\mathbb{H})$
\begin{align*}
  \tilde{\cE}_{0, \ga}& =\int_{\mathbf{\Phi}(\mathbb{H})}(R^*-|\tilde{x}|)^{\ga} |\tilde{L} \tilde\phi|^2
  +|{\tilde{\Lb}}\tilde{\phi}|^2+|\tilde{\nabb}\tilde{\phi}|^2+|\tilde{\phi}|^2 +(R^*-|\tilde{x}|)^{p-3+\ga}|\tilde{\phi}|^{p+1}d\tilde{x},
\end{align*}
where
\[
\tilde{\Om}_{ij}=\tilde{x}_i\tilde{\pa}_{j}-\tilde{x}_j\tilde{\pa}_{i}=x_i\pa_j-x_j\pa_i=\Om_{ij}
\]
and $\tilde{\pa}$ is the full derivative on $\mathbf{\Phi}(\mathbf{D})$.
Direct computations imply the following change of null frames:
\begin{align*}
&L=(t^*+r)^{-2}\tilde{L
,\quad\Lb=(t^*-r)^{-2}\tilde{\Lb},\qua
\pa_i-\om_i \pa_r=((t^*)^2-r^2)^{-1}(\pa_{\tilde{x}_i}-\tilde{\om}_i\pa_{\tilde{r}})
\end{align*}
We first prove the following bound.
\begin{Prop}
  \label{prop:bd:ID:comp}
  Let $\ga=2-\ga_0+\ep$. The solution $\tilde{\phi}= \La\phi$ on $\mathbf{\Phi}(\mathbf{D})$ verifies the following bounds
  \begin{align}
    \label{eq:bd:IDSB:comp}
   \tilde{\cE}_{0,\ga}+
  \iint_{\mathbf{\Phi}(\mathbf{D})} \La^{3-p}|\tilde{\phi}|^{p+1}(R^*-\tilde{t}-|\tilde{x}|)^{\ga-1}  d\tilde{x}d\tilde{t}
   \leq C\mathcal{E}_{0, \ga_0}[\phi].
  \end{align}
  for some constant C depending on $\ga_0$, $p$ and $\ep$.
\end{Prop}
\begin{proof}
On the hyperboloid $\mathbb{H}$, the coordinate functions verify the following relation
\begin{align*}
  (t^*-(2R^*)^{-1})dt=rdr,\quad (t^*)^2-r^2=(R^*)^{-1}t^*,
\end{align*}
which implies that
\begin{align*}
d\tilde{r}=\La^{-2}((t^*)^2+r^2)dr-2 \La^{-2}t^* r dt=\f12 (t^*)^{-1}(t^*-(2R^*)^{-1})^{-1}dr.
\end{align*}
Since $\tilde{\om}=\om$, $\tilde{r}=\Lambda^{-1}r$, the surface measure obeys
\begin{align*}
  d\tilde{x}=\tilde{r}^2d\tilde{r}d\tilde{\om}=\f12 \tilde{r}^2 (t^*)^{-1}(t^*-(2R^*)^{-1})^{-1}drd\om=\f12 \Lambda^{-2}(t^*)^{-1}(t^*-(2R^*)^{-1})^{-1} dx.
\end{align*}
We also note that
on the hyperboloid $\mathbb{H}$
\[
0\leq t^*-r=\frac{(R^*)^{-1}t^*}{t^*+r}\leq (R^*)^{-1},\quad  \La=(R^*)^{-1}t^*,\quad t^*\geq (R^*)^{-1},
\]
which leads to the following bounds
\begin{align*}
  d\tilde{x}\les \La^{-4}dx,\quad t^*+r\les \La\les t^*+r.
\end{align*}
Here in only in proof for this proposition, the implicit constant relies only on $\ga_0$, $p$ and $\ep$.
For the zeroth order weighted energy $\tilde{\cE}_{0, \ga}$, by definition we can estimate that
\begin{align*}
  &(R^*-|\tilde{x}|)^{\ga} |{\tilde{L}}\tilde{\phi}|^2 +|\tilde{\nabb} \tilde{\phi}|^2
  + |\tilde{ \phi}|^2+|{\tilde{\Lb}}\tilde{\phi}|^2 +(R^*-|\tilde{x}|)^{p-3+\ga}|\tilde{\phi}|^{p+1} \\
  &=(t^*+r)^{4-\ga}|{L}(\La\phi)|^2+\Lambda^4|\nabb\phi|^2+(t^*-r)^{4}|{\Lb}(\La\phi)|^2+|\La \phi|^2+(t^*+r)^{-p+3-\ga}|\La\phi|^{p+1}\\
   &\les \La^{2+\ga_0}(|{L}(r\phi)|^2+|{L}\phi|^2+|\phi|^{p+1})+\Lambda^4|\nabb \phi|^2+ \La^2(
  |{\Lb}\phi|^2+ |\phi|^2).
\end{align*}
Here $R^*-|\tilde{x}|=\La^{-1}(R^* \La-r)=t^*+r$.
Note that the classical energy flux through the hyperboloid $\mathbb{H}$ verifies the lower bound
\begin{align*}
  \int_{\mathbb{H}}   \La^{-2}|{\Lb} \phi|^2+|L \phi|^2+|\nabb \phi|^2+\frac{2}{p+1}|\phi|^{p+1}
  dx\les  E[\phi](\mathbb{H}).
\end{align*}
By using the bound $d\tilde{x}\les \La^{-4}dx$, we therefore can estimate that
\begin{align*}
  \tilde{\cE}_{0, \ga}&\les \int_{\mathbb{H}}\big(\La^{2+\ga_0}(|{L}(r\phi)|^2+|{L}\phi|^2+|\phi|^{p+1})+\Lambda^4|\nabb \phi|^2+ \La^2(
  |{\Lb}\phi|^2+ |\phi|^2) \big)\Lambda^{-4}dx\\
  &\les E[\phi](\mathbb{H})+ \int_{\mathbb{H}^+} r^{\ga_0-2} |{L}(r\phi)|^2 dx\\

  &\les {\cE}_{0, \ga_0}[\phi].
\end{align*}
Here the integral of $|\phi|^2$ is estimated by using Hardy's inequality and the $r$-weighted energy estimate through the hyperboloid $\Hy$ follows from Proposition \ref{prop:energyflux:H:ex}.

\iffalse
For higher order weighted energy $\tilde{\cE}_{1, \ga_1}$, by definition, we need to estimate $|\tilde{\pa}^{l_1+1} {\tilde{\Om}_{ij}}^{l_2} \tilde{\phi}|$ in $(\tilde{t}, \tilde{x})$ coordinates in terms of the associated quantities in $(t, x)$ coordinates for all $l_1+l_2=1$. First notice that the region inclosed by $\Hy^{-}$ and the initial hypersurface is compact. The decay estimates of Proposition \ref{prop:NLW:def:3d:ex} hold on this region by shifting the origin of the Minkowski space. In particular, the solution $\phi$ is uniformly bounded in this compact region. Therefore by using the standard energy estimate, we conclude that
\begin{align*}
  E[Z\phi](\Hy^{-})\les \mathcal{E}_{1, \ga_0}[\phi]^{p-1},\quad \forall Z\in \{\pa_\mu, \Om_{\mu\nu}\}.
\end{align*}
Now as $$\tilde{\Om}_{ij}=\Om_{ij}, \quad \Om_{ij}(\La)=\Om_{ij}(t^*+r)=\Om_{ij}(t^*-r)=0, $$
the above argument for $\tilde{\cE}_{0, \ga_1}$ also holds for the case when $l_1=0$, that is,
 \begin{align*}
   \int_{\mathbf{\Phi}(\mathbb{H})}|\tilde{\pa} {\tilde{\Om}_{ij}} \tilde{\phi}|^2(R^*-|\tilde{x}|)^{\ga_1}d\tilde{x} \les E[\Om_{ij}\phi](\mathbb{H})+ \int_{\mathbb{H}^+} r^{\ga_0-2} |{L}(r\Om_{ij}\phi)|^2 dx\les \mathcal{E}_{1, \ga_0}[\phi]^{p-1}.
 \end{align*}
 Here we used the weighted energy estimate \eqref{eq:EE:hyB} for $Z\phi$.

 For the case when $l_1=1$, $l_2=0$, for any scalar field $f$ defined on $\bf{\Phi}(\mathbb{H})$, we can estimate that
\begin{equation*}
\begin{split}
  |\tilde{\pa}f|^2&=|{\tilde{\Lb}}f|^2+|{\tilde{L}}f|^2+|\tilde{\nabb}f|^2\\
  &=|(t^*+r)^2 L f(\mathbf{\Phi})|^2+|(t^*-r)^2 {\Lb} f(\mathbf{\Phi})|^2+|\La \nabb f(\mathbf{\Phi})|^2\\
  &\les \sum\limits_{\mu, \nu}(t^*+r)^2(|\Om_{\mu\nu} f(\mathbf{\Phi})|^2+|\pa f(\mathbf{\Phi})|^2)\\
  &\les \La^2 \sum\limits_{ Z\in \Ga}|Z f(\mathbf{\Phi})|^2.
  \end{split}
\end{equation*}
We remark here that the above estimate holds only restricted to the hyperboloid $\Hy$ or $\mathbf{\Phi}(\Hy)$. Moreover
let $\mathcal{Z}$ be the push forward of the vector field $Z$ under the map $\bf{\Phi}$. For any function $f$ defined on $\bf{\Phi}(\mathbb{H})$, we have
\begin{align*}
  |[\mathcal{Z}, \tilde{\pa}]f|\les \sum\limits_{Z\in \Ga}\La |Z f(\bf{\Phi})|.
\end{align*}
Therefore by repeating the above process, we can show that
\begin{align*}
  |\tilde{\pa}^{2} \tilde{\phi}|^2
  &\les \La^2\sum\limits_{Z\in \Ga}|Z  (\tilde{\pa}  \tilde{\phi})(\Phi)|^2\\
  &\les \La^2 \sum\limits_{Z\in \Ga}|\tilde{\pa}\mathcal{Z}  \tilde{\phi}|^2+|[\tilde{\pa},\mathcal{Z} ] \tilde{\phi}|^2\\
  &\les \La^4 \sum\limits_{Z\in \Ga}|Z^2  (\La \phi)|^2+|Z (\La\phi)|^2\\
  &\les \La^{2 } \sum\limits_{k\leq 1, Z\in \Ga}|(t^*+r)^2 L Z^{k}   (\La\phi)|^2+|(t^*-r)^2 \Lb Z^{k}   (\La\phi)|^2+|\La \nabb Z^{k}  (\La\phi)|^2.
\end{align*}
Therefore we can estimate that
\begin{align*}
   &\int_{ \mathbf{\Phi}(\Hy)}|\tilde{\pa}^{2}  \tilde{\phi}|^2(R^*-|\tilde{x}|)^{\ga_1+2}d\tilde{x}\\
   &\les \sum\limits_{l\leq 1, Z\in \Ga}\int_{\mathbb{H}} \La^{2} \big( |(t^*+r)^2 L Z^{l} (\La\phi)|^2+|(t^*-r)^2 \Lb Z^{l}(\La\phi)|^2+|\La \nabb Z^{l} (\La\phi)|^2\big)(t^*+r)^{-\ga_1-2}\La^{-4}dx\\
   &\les \sum\limits_{l\leq 1, Z\in \Ga}\int_{\mathbb{H}}  (1+r)^{\ga_0-2}| L {Z}^l(r\phi)  |^2 +\La^{-2}(| {\Lb }{Z}^l \phi  |^2+| {Z}^l \phi  |^2)+  |L {Z}^l \phi |^2+  |\nabb {Z}^l \phi |^2dx\\
   &\les \sum\limits_{l\leq 1} E[Z^l\phi](\Hy)+\int_{\Hy^+} (1+r)^{\ga_0-2}| L {Z}^l(r\phi)  |^2\les \mathcal{E}_{1, \ga_0}[\phi]^{p-1}.
\end{align*}
We thus complete the proof for the bound \eqref{eq:bd:ID:comp}.
\fi

For the spacetime $\mathcal{I}$, Proposition \ref{prop:spacetime:bd} in particular implies that
\begin{align*}
  \iint_{\mathbf{D}} |\phi|^{p+1}v_+^{\ga_0-1-\ep} dxdt\les \mathcal{E}_{0, \ga_0}[\phi]
\end{align*}
for the solution $\phi$ to \eqref{eq:NLW:semi:3d}.
Since the map $\mathbf{\Phi}$ is conformal and $\La$ is the conformal factor, we conclude that
\begin{align*}
  d\tilde{x}d\tilde{t}=\La^{-4} dxdt,
\end{align*}
which can also be derived by direct computations. Thus the associated scalar field $\tilde{\phi}$ on $\mathbf{\Phi}$ verifies the following weighted bound
\begin{align*}
  \iint_{\mathbf{\Phi}(\mathbf{D})} \La^{3-p}|\tilde{\phi}|^{p+1}(R^*-\tilde{t}-|\tilde{x}|)^{-\ga_0+1+\ep}  d\tilde{x}d\tilde{t}&\les\iint_{\mathbf{\Phi}(\mathbf{D})} \La^{-p-1}|\tilde{\phi}|^{p+1}v_+^{\ga_0-1-\ep}\La^4  d\tilde{x}d\tilde{t}\\
  &=\iint_{\mathbf{D}} |\phi|^{p+1}v_+^{\ga_0-1-\ep} dxdt\les \mathcal{E}_{0, \ga_0}[\phi].
\end{align*}
 We thus finished the proof for the Proposition.
\end{proof}

We now prove the main theorem \ref{thm:main} by showing the pointwise decay estimates for the solution $\phi$ to \eqref{eq:NLW:semi:3d} in the interior region. As indicated previously,  $\tilde{\phi}=\La\phi$ solves equation \eqref{eq:NLW:3D:conf} on the compact region $\mathbf{\Phi}(\mathbf{D})$ for solution $\phi$ to \eqref{eq:NLW:semi:3d}.

In view of the previous Proposition \ref{prop:bd:ID:comp}, we derive that
\begin{align*}
  \tilde{\cE}_{0, \ga}+\mathcal{I}\les \mathcal{E}_{0, \ga_0}[\phi].
\end{align*}
By our assumption on $\ga_0$ and the choice of $\ep$, we always have $0<\ga<1$. For the case when
$$\frac{1+\sqrt{17}}{2}<p<5, \quad \max\{\frac{4}{p-1}-1, 1\}<\ga_0<\min\{p-1, 2\},$$
the choice of $\ep$ also implies that
\begin{align*}
  (p-1)(3-\ga)=(p-1)(1+\ga_0-\ep)>4.
\end{align*}
Then from Proposition \ref{Prop:pointwise:decay:EM}, we conclude that
\begin{align*}
  |\La\phi|(\tilde{t}, \tilde{x})&\les  \sup\limits_{|\tilde{y}|\leq \tilde{t}}|\tilde{\phi}^{lin}(\tilde{t}, \tilde{y})|

\end{align*}
Here $\tilde{\phi}^{lin}$ is the linear evolution with initial data $(\tilde{\phi}(0, \tilde{x}), \tilde{\pa}_{\tilde{t}}\tilde{\phi}(0, \tilde{x}))$
By the conformal transformation, $\tilde{\phi}^{lin}$ can be identified with $\La \phi_H^{lin}$, in which $\phi_H^{lin}$ was defined before Proposition \ref{prop:NLW:def:3d:ex:ID:H}.
Recall that
\[
\tilde{t}=R^*-\La^{-1}(t+3),\quad |\tilde{x}|=\La^{-1}r, \quad \La =(t+3-r)(t+3+r).
\]
And inside the hyperboloid $\mathbb{H}$, we have $v_+\leq t+3$. Thus
\begin{align*}
 \frac{1}{8}u_+^{-1}&\leq  R^*-\tilde{t}=\La^{-1}(t+3)\leq u_+^{-1},\\
 \frac{1}{4}v_+^{-1}&\leq  R^*-\tilde{t}-|\tilde{x}|=\La^{-1}(t+3-r)=(t+3+r)^{-1}\leq  v_+^{-1}.
\end{align*}
Therefore by using the decay estimate \eqref{eq:phi:pt:Br:largep:lin} of Proposition \ref{prop:NLW:def:3d:ex:ID:H}, we then conclude that
\begin{align*}
|\tilde{\phi}^{lin}(\tilde{t}, \tilde{x})|&\les \La |\phi^{lin}_H|\les \mathcal{E}_{1, \ga_0}[\phi]^{\frac{p-1}{2}} (2+t+|x|)^{-1}(2+||x|-t|)^{-\frac{\ga_0-1}{2}}\La\\
&\les \mathcal{E}_{1, \ga_0}[\phi]^{\frac{p-1}{2}} (R^*-\tilde{t})^{-\frac{1+\ga_0}{2}},
\end{align*}
which leads to
\begin{align*}
  |\phi(t, x)|\les \mathcal{E}_{1, \ga_0}[\phi]^{\frac{p-1}{2}} \La^{-1}(R^*-\tilde{t})^{-\frac{1+\ga_0}{2}} \les \mathcal{E}_{1, \ga_0}[\phi]^{\frac{p-1}{2}}v_+^{-1}u_+^{-\frac{\ga_0-1}{2}}.
\end{align*}
This proves the pointwise decay estimate for the solution in the interior region for the large $p$ case.

Finally for the small $p$ case, take
\[
\b=\frac{\ga_0}{p+1}.
\]
The small positive constant $\ep$ can be chosen so that $\b\leq \frac{p-1}{p+1}\ga_0-1$ as $\ga_0>1$.
Then by using the linear decay estimate \eqref{eq:phi:pt:Br:smallp:lin} of Proposition \ref{prop:NLW:def:3d:ex:ID:H}, we can show that
\begin{align*}
  |\tilde{\phi}^{lin}u_*^{1-\b} v_*^{1-\a_p \ga_0}|&\les |\La \phi_{H}^{lin} u_+^{\b-1} v_+^{-1+\a_p \ga_0}|\\
  &\les \sqrt{\mathcal{E}_{1, \ga_0}[\phi] } u_+ v_+ v_+^{-\a_p \ga_0}u_+^{-\b}u_+^{\b-1} v_+^{-1+\a_p \ga_0}\\
  &\les \sqrt{\mathcal{E}_{1, \ga_0}[\phi] }.
\end{align*}
Here recall that $u_*=R^*-\tilde{t}$, $v_*=R^*-\tilde{t}-|\tilde{x}|$. Hence from estimate \eqref{eq:phi:pt:Br:EM:smallp} of Proposition \ref{Prop:pointwise:decay:EM}, we conclude that
\begin{align*}
  |\tilde{\phi}(\tilde{t}_0, \tilde{x}_0)|&\les (1+\sup\limits_{\tilde{t}+|\tilde{x}|\leq R*}|\tilde{\phi}^{lin}u_*^{1-\b} v_*^{1-\a_p \ga_0}|)(R^*-\tilde{t}_0)^{-1+\b} (R^*-\tilde{t}_0-|\tilde{x}_0|)^{-1+\a_p\ga_0}\\
  &\les (1+\sqrt{\mathcal{E}_{1, \ga_0}[\phi]}) (R^*-\tilde{t}_0)^{-1+\b} (R^*-\tilde{t}_0-|\tilde{x}_0|)^{-1+\a_p\ga_0},
\end{align*}
which implies that in the interior region for the case when $2<p\leq \frac{1+\sqrt{17}}{2}$, the solution $\phi$ of \eqref{eq:NLW:semi:3d} verifies the following decay estimate
\begin{align*}
  |\phi(t, x)|\leq |\La^{-1}\tilde{\phi}(\tilde{t}, \tilde{x})|\les (1+\sqrt{\mathcal{E}_{1, \ga_0}[\phi]}) u_+^{-1}v_+^{-1}u_+^{1-\b} v_+^{1-\a_p\ga_0}\les (1+\sqrt{\mathcal{E}_{1, \ga_0}[\phi]}) u_+^{-\b} v_+^{-\a_p\ga_0}.
\end{align*}
Here the implicit constant replies on $\mathcal{E}_{0, \ga_0}[\phi]$, $\ga_0$ and $p$.
We thus complete the proof for the main Theorem \ref{thm:main}.

\bigskip

As for the scattering result of Corollary \ref{cor:scattering:3D}, by using the standard energy estimate, the solution scatters in energy space if the mixed norm $\|\phi\|_{L_t^p L_x^{2p}}$ of the solution is finite (see e.g. Lemma 4.4 in \cite{Strauss78:NLW}). Moreover, it has been shown in the author's companion paper, the solution scatters in $\dot{H}^{s}$ for all $\frac{3}{2}-\frac{2}{p-1}\leq s\leq 1$ for the case when $p>\frac{1+\sqrt{17}}{2}$. In particular, it suffices to consider the small $p$ case when $2<p\leq \frac{1+\sqrt{17}}{2}$. By using the pointwise decay estimate of the main theorem \ref{thm:main}, we estimate that
\begin{align*}
  \|\phi\|_{L_t^p L_x^{2p}}^p&=\int_{\mathbb{R}} \left(\int_{\mathbb{R}^3} |\phi|^{p+1}v_+^{\ga_0-1-\ep}|\phi|^{p-1}v_+^{-\ga_0+1+\ep}dx\right)^{\f12} dt \\
  &\les \int_{\mathbb{R}} \left(\int_{\mathbb{R}^3} |\phi|^{p+1}v_+^{\ga_0-1-\ep}(1+t)^{-\ga_0+1+\ep-(p-1)\a_p \ga_0}dx\right)^{\f12} dt\\
  &\les \left(\int_{\mathbb{R}}\int_{\mathbb{R}^3} |\phi|^{p+1}v_+^{\ga_0-1-\ep}dx dt\right)^{\f12}\left(\int (1+|t|)^{-\ga_0+1+\ep-(p-1)\a_p \ga_0}dt\right)^{\f12}.
\end{align*}
In view uniform spacetime bound of Proposition \ref{prop:spacetime:bd}, $\|\phi\|_{L_t^p L_x^{2p}}$ is finite if
\begin{align*}
  \ga_0-1+(p-1)\a_p \ga_0>1
\end{align*}
by choosing $\ep$ sufficiently small. As $1<\ga_0<p-1$, by choosing $\ga_0$ sufficiently close to $p-1$, it is equivalent to that
\begin{align*}
  f(p)=p-2+(p-1)^2 \frac{3+(p-2)^2}{(5-p)(p+1)}-1>0.
\end{align*}
It can be checked that there is a unique solution $p_*$ of $f(p)$ on $[2, 3]$ and when $p>p_*$, one has $f(p)>0$. Numerically, one can show that
\[
2.3541<p_*<2.3542.
\]
This proves the scattering result of Corollary \ref{cor:scattering:3D}.