For $\mathcal{R} \gg \mathcal{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_{\mathcal{R}} = t - \frac{\mathcal{R}} {c} $$ The result is: $$ P\left(\mathcal{R},t\right) = \frac {1} {6 \pi \epsilon_0 c^2} \left( \frac {\vec{Y}\left(t_{\mathcal{R}}\right)\cdot\vec{Z}\left(t_{\mathcal{R}}\right)} {\mathcal{R}} + \frac {\left|\vec{Z}\left(t_{\mathcal{R}}\right)\right|^2} {c} \right) $$
Derivation
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} {\mathcal{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 $\mathcal{R} \gg \mathcal{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 substitute $\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 $(\mathscr{r},\theta,\phi),\space 0\le \mathscr{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)$ is $( |\vec{Z}\left(t\right)|, 0, 0)$ and $\vec{Y}\left(t\right)$ is $( |\vec{Y}\left(t\right)|, \gamma\left(t), 0\right)$. Note that the orientation of this coordinate system changes with time. In this coordinate system, therefore, $\hat{n}\left(\vec{r}\right)$ at any instant $t$ has to be represented in terms of the retarded time as $\left(1,\vartheta\left(\vec{r},t_r\right), \varphi\left(\vec{r},t_r\right)\right)$
Working out $\vec{S}\left(\vec{r},t\right)\cdot\hat{n}\left(\vec{r}\right)$ in this coordinate system, we get: $$ \vec{S}\left(\vec{r},t\right)\cdot\hat{n}\left(\vec{r}\right) = \frac {\left|\vec{Z}\left(t_{\mathcal{R}}\right)\right|^2} {16 \pi^2 \epsilon_0 {\mathcal{R}}^2 c^3} \sin^2\vartheta\left(\vec{r},t_{\mathcal{R}}\right) + \frac {\vec{Y}\left(t_{\mathcal{R}}\right)\cdot\vec{Z}\left(t_{\mathcal{R}}\right)} {16 \pi^2 \epsilon_0 {\mathcal{R}}^3 c^2} \sin^2\vartheta\left(\vec{r},t_{\mathcal{R}}\right) \\ - \frac {\left|\vec{Y}\left(t_{\mathcal{R}}\right)\times\vec{Z}\left(t_{\mathcal{R}}\right)\right|} {32 \pi^2 \epsilon_0 {\mathcal{R}}^3 c^2} \sin2\vartheta\left(\vec{r},t_{\mathcal{R}}\right) \cos\varphi\left(\vec{r},t_{\mathcal{R}}\right) $$
Finally, we compute the surface integral. The third term vanishes, and for a sphere $V$ of radius $\mathcal{R}$ we get: $$ \oint_{\partial V} \sin^2\vartheta\left(\vec{r}\right) \space ds(\vec{r}) = \frac {8 \pi \mathcal{R}^2} {3} $$
Substituting, we get our result.
Observations
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.
Also, $$ \mathcal{P}_{\infty}\left(t\right) = \frac {\left|\vec{Z}\left(t_{\mathcal{R}}\right)\right|^2} {6 \pi \epsilon_0 c^3} ,\qquad \mathcal{P}_{\delta}\left(\mathcal{R},t\right) = \frac {\vec{Y}\left(t_{\mathcal{R}}\right)\cdot\vec{Z}\left(t_{\mathcal{R}}\right)} {6 \pi \epsilon_0 \mathcal{R} c^2} $$