For $\mathfrak{R} \gg \mathfrak{R}_s$, if we define:
$$
\vec{Y}\left(t\right)=
   \iiint_{V_s} 
      \vec{J}\left(\vec{r},t\right) 
   \space dV\left(\vec{r}\right)
   ,\qquad
\vec{Z}\left(t\right)=
   \iiint_{V_s} 
      \frac{\partial \vec{J}\left(\vec{r},t\right)} {\partial t} 
   \space dV\left(\vec{r}\right)
 = \frac {d \vec{Y}\left(t\right)} {dt}
$$
and also
$$
t_{\mathfrak{R}} = t - \frac{\mathfrak{R}} {c}
$$
The result is:
$$
P\left(\mathfrak{R},t\right) =
  \frac {1} {6 \pi \epsilon_0 c^2}
  \left(
    \frac 
       {\vec{Y}\left(t_{\mathfrak{R}}\right)\cdot\vec{Z}\left(t_{\mathfrak{R}}\right)}
       {\mathfrak{R}}
    +
    \frac 
       {\left|\vec{Z}\left(t_{\mathfrak{R}}\right)\right|^2} 
       {c}
  \right)
$$

The derivation goes something like the following. Given the following definitions:
$$
\vec{R} = \vec{r} - \vec{r}_s, \qquad R = \left | \vec{R} \right |, \qquad t_r = t - \frac {R} {c}
$$
From Jefimenko's Equations, the value of $\vec{E}$ and $\vec{B}$ for any $\vec{r}$ and any $t$ is as follows:
$$
\vec{E}(\vec{r},t) =
\frac {1} {4 \pi \epsilon_0}
\iiint_{V_s} {\left(
\frac {\rho (\vec{r}_s, t_r)} {R^3} \vec{R} +
\frac {1} {R^2 c} \frac {\partial \rho (\vec{r}_s, t_r) } {\partial t} \vec{R} -
\frac {1} {R c^2} \frac {\partial \vec{J} (\vec{r}_s, t_r) } {\partial t} \right)}
\space dV\left(\vec{r}_s\right)
$$
$$
\vec{B}(\vec{r},t) =
\frac {\mu_0} {4 \pi}
\iiint_{V_s} {\left(
\frac {\vec{J}  (\vec{r}_s, t_r)} {R^3} \times \vec{R} +
\frac {1} {R^2 c} \frac {\partial \vec{J}  (\vec{r}_s, t_r) } {\partial t} \times \vec{R} \right)}
\space dV\left(\vec{r}_s\right)
$$
Using this, we can write down the expression for $\vec{S}\left(\vec{r},t\right)$, which has six terms, each of which is a product of two volume integrals.

For the system defined in the question, we designate the center of $V$ (and $V_s$) as $\vec{r}_0$. For all $\vec{r}\in\partial V$ and all $\vec{r}_s\in V_s$, we also define:
$$
\hat{n}\left(\vec{r}\right) = \frac {\vec{r}-\vec{r}_0} {\left|\vec{r}-\vec{r}_0\right|}
,\qquad
k\left(\vec{r}, \vec{r}_s\right) = \frac {R} {\mathfrak{R}} =\frac {\left|\vec{r} - \vec{r}_s\right|} {\left|\vec{r} - \vec{r}_0\right|}, 
\qquad 
h\left(\vec{r}, \vec{r}_s\right) = 
  \frac {\left(\vec{r} - \vec{r}_s\right)} {\left|\vec{r} - \vec{r}_s\right|} 
  \cdot
  \hat{n}\left(\vec{r}\right)
$$
When $\mathfrak{R} \gg \mathfrak{R}_s$, it can be seen that $k\left(\vec{r}, \vec{r}_s\right)\approx 1$ and $h\left(\vec{r}, \vec{r}_s\right)\approx 1$ irrespective of $\vec{r}_s$. Under this 'far field' approximation, therefore, $\vec{R}\approx\vec{r}-\vec{r}_0$, and is effectively independent of $\vec{r}_s$, so it can be taken outside of the volume integrals. 

Then, because $\left(\vec{a}\times\left(\vec{b}\times\vec{a}\right)\right)\cdot\vec{a}=0$ four of the six terms in $\vec{S}\left(\vec{r},t\right)\cdot\hat{n}\left(\vec{r}\right)$ turn out to be zero.

Next, we sustitute $\vec{Y}$ and $\vec{Z}$ into the expression for $\vec{S}\left(\vec{r},t\right)\cdot\hat{n}\left(\vec{r}\right)$

Now we define a spherical coordinate system $\left(r,\theta,\phi\right),\space 0\le r\lt\infty,\space 0\le\theta\le\pi,\space 0\le\phi\le 2\pi$ with its center at $\vec{r}_0$, and orient it such that:
$$
\vec{Z}\left(t\right) : \left( \left|\vec{Z}\left(t\right)\right|, 0, 0\right)
$$
$$
\vec{Y}\left(t\right) : \left( \left|\vec{Y}\left(t\right)\right|, \gamma\left(t\right), 0\right)
$$
Note that the orientation of this coordinate system changes with time. In this coordinate system, therefore, we will write:
$$
\hat{n}\left(\vec{r}\right) : \left(1,\vartheta\left(\vec{r},t_r\right), \varphi\left(\vec{r},t_r\right)\right)
$$

Working out the vectors, we get:
$$
\vec{S}\left(\vec{r},t\right)\cdot\hat{n}\left(\vec{r}\right)  = 
\frac {1} {16 \pi^2 \epsilon_0 R^2 c^2} \left(
  \frac {\left|\vec{Z}\left(t_r\right)\right|^2} {c}
    \sin^2\vartheta\left(\vec{r},t_r\right)\space
   \\ +  \frac {\vec{Y}\left(t_r\right)\cdot\vec{Z}\left(t_r\right)} {R}
    \sin^2\vartheta\left(\vec{r},t_r\right)
   \\ -  \frac {\left|\vec{Y}\left(t_r\right)\times\vec{Z}\left(t_r\right)\right|} {2R}
    \sin2\vartheta\left(\vec{r},t_r\right)
    \cos\varphi\left(\vec{r},t_r\right)
\right)
$$
Finally, computing the surface integral yields the result.

For the record,because of the way we've defined $V_s$ and because of conservation of charge, $Q_s$ is the (constant) total charge inside $V_s$,
$$
   \iiint_{V_s} 
      \rho\left(\vec{r},t\right) 
   \space dV\left(\vec{r}\right) = Q_s
   ,\qquad
   \iiint_{V_s} 
      \frac{\partial \rho\left(\vec{r},t\right)} {\partial t} 
   \space dV\left(\vec{r}\right)
 = \frac {d Q_s} {dt}
 = 0
$$
... but we don't need these terms in the derivation.