Radiated power from a 'source volume' outside which charge & current are zero (i.e. derive radiated power from Jefimenko's equations) In classical electrodynamics, what is the radiated power from a generalized source (consisting of charge density $\rho$ and current density $\vec{J}$) in vacuum?
Let us define $V_s$ to be the smallest sphere that completely encloses all source points -- that is, it is the smallest sphere outside which $\rho\left(\vec{r},t\right) = \vec{J}\left(\vec{r},t\right) = 0$ for all $t$. Let us say $\mathcal{R}_s$ is the radius of $V_s$.
Now let us define another sphere $V$ concentric with $V_s$. Let us say the radius of this sphere is $\mathcal{R}$, and $\mathcal{R} > \mathcal{R}_s$.
So what is power crossing the surface of $V$ at any instant? In other words, what is the value of the following?
$$
P\left(\mathcal{R},t\right) =
    \oint_{\partial V}
      \vec{S}\left(\vec{r},t\right) 
      \cdot
    \space d\vec{s}\left(\vec{r}\right) 
    = \frac{1}{\mu_0}
    \oint_{\partial V}
      \left( 
        \vec{E}\left(\vec{r},t\right) 
        \times 
        \vec{B}\left(\vec{r},t\right)
      \right)
      \cdot
      \space d\vec{s}\left(\vec{r}\right)
$$
Below, I have provided an answer for the 'far field' approximation ($\mathcal{R} \gg \mathcal{R}_s$). Are there 'sanitized' expressions for cases where this approximation does not hold? 
Even if there isn't a 'clean' expression for the general case, what I am looking for is the following. Let us say that there are two instants $t_1$ and $t_2$ such that for all $\vec{r} \in V_s$ we have $\rho\left(\vec{r}, t_1\right) = \rho\left(\vec{r}, t_2\right)$ and $\vec{J}\left(\vec{r}, t_1\right) = \vec{J}\left(\vec{r}, t_2\right)$. If we now define:
$$
\mathcal{P}_{\infty}\left(t\right) = \lim_{\mathcal{R}\to\infty} P\left(\mathcal{R},t + \frac {\mathcal{R}} {c}\right) \\
\mathcal{P}_{\delta}\left(\mathcal{R},t\right) = P\left(\mathcal{R},t + \frac {\mathcal{R}} {c}\right) - 
\mathcal{P}_{\infty}\left(t\right)
$$
... can the following be proven? 
$$
\int_{t_1}^{t_2} \mathcal{P}_{\delta}\left(\mathcal{R},t\right)\space dt = 0
$$
Also, is it true (and can it be proven) that if $\hat{n}\left(\vec{r}\right)$ is the unit normal at any $\vec{r} \in V$, there will always be one axis along which $\vec{S}\left(\vec{r},t\right) \cdot \hat{n}\left(\vec{r}\right) = 0$?
 A: 
There are some problems with this derivation, it kind of works only when $R\to\infty$ and when $\dfrac {\partial\vec{J}} {\partial t}$ is very small. 

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}
$$
Assuming that the origin is at the center of $V$, so that $\hat{n}\left(\vec{r}\right) = \hat{r}$, the results are:
$$
\vec{S}\left(\vec{r},t\right)\cdot\hat{n}\left(\vec{r}\right) = 
  \frac {
  \left( \vec{Y}\left(t_\mathcal{R}\right) \times \hat{n}\left(\vec{r}\right) \right)
  \cdot 
  \left(\vec{Z}\left(t_\mathcal{R}\right) \times \hat{n}\left(\vec{r}\right) \right) 
  } {16 \pi^2 \epsilon_0\mathcal{R}^3 c^2}
  + 
  \frac {
  \left| \vec{Z}\left(t_\mathcal{R}\right)  \times \hat{n}\left(\vec{r}\right) \right|^2
  } {16 \pi^2 \epsilon_0\mathcal{R}^2 c^3}
\\
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)
\\
\mathcal{P}_{\infty}\left(t\right) =
  \frac 
    {\left|\vec{Z}\left(t\right)\right|^2} 
    {6 \pi \epsilon_0 c^3}
,\qquad
\mathcal{P}_{\delta}\left(\mathcal{R},t\right) =
    \frac 
       {\vec{Y}\left(t\right)\cdot\vec{Z}\left(t\right)}
       {6 \pi \epsilon_0 \mathcal{R} c^2}
$$
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. We are, however, interested in $\vec{S}\left(\vec{r},t\right) \cdot \hat{n}\left(\vec{r}\right)$, where, if we designate the center of $V$ (and $V_s$) as $\vec{r}_0$:
$$
\hat{n}\left(\vec{r}\right) = \frac {\vec{r}-\vec{r}_0} {\left|\vec{r}-\vec{r}_0\right|}
$$
Now if we choose a frame of reference where $\vec{r}_0$ is the origin, we find that $\hat{n}\left(\vec{r}\right) = \hat{r}$.
When $\mathcal{R} \gg \mathcal{R}_s$, we can consider $\dfrac{\left|\vec{r}_s\right|}{\left|\vec{r}_s\right|} \approx 0$. Thus, irrespective of $\vec{r}_s$ under this 'far field' approximation, $\vec{R}\approx\vec{r}-\vec{r}_s$, 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)$ and get the expression listed above.
To help us with the surface integral, 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. 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},\qquad
\oint_{\partial V} \sin 2\vartheta\left(\vec{r}\right) \space\cos\varphi\left(\vec{r}\right) \space ds(\vec{r})= 0
$$
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, $\int_{t_1}^{t_2} \mathcal{P}_{\delta}\left(\mathcal{R},t\right) \space dt = 0$ because:
$$
\int_{t_1}^{t_2} \vec{Y}(t)\cdot\vec{Z}(t) dt = \int_{t_1}^{t_2} \vec{Y}(t)\cdot\frac{d\vec{Y}}{dt}(t) dt = \frac12 \int_{t_1}^{t_2} \frac{d}{dt}|\vec{Y}(t)|^2 dt = \frac{1}{2} \left[|\vec{Y}(t)|^2\right]_{t_1}^{t_2}
$$
and $$\forall\vec{r}\in V: \vec{J}\left(\vec{r},t_1\right) = \vec{J}\left(\vec{r},t_2\right) \implies \vec{Y}\left(t_1\right) = \vec{Y}\left(t_2\right)$$
The value of $\vec{S}\left(\vec{r},t\right)\cdot\hat{n}\left(\vec{r}\right)$ is zero whenever $\vec{Z}\left(t_{\mathcal{R}}\right) \times \hat{n}\left(\vec{r}\right) = 0$, that is, wherever the unit normal is parallel to $\vec{Z}\left(t_{\mathcal{R}}\right)$
A: I realized that $\vec{W}$, $\vec{Y}$ and $\vec{Z}$ depend on both $\vec{r}$ and $t$ independently, anot not through a composite $t_r$. So here is an enhanced version of the answer
Let us say $\mathbb{V}_s \subset \mathbb{R}^3$ is the smallest spherical volume such that:
$$
\forall t, \qquad \forall \vec{r}\notin \mathbb{V}_s: \qquad\qquad \rho\left(\vec{r},t\right) = \vec{J}\left(\vec{r},t\right) = 0
$$
$$
\forall t, \qquad \forall \vec{r}\in \partial\mathbb{V}_s: \qquad\qquad \rho\left(\vec{r},t\right) = \vec{J}\left(\vec{r},t\right) = 0 
$$
Since $\mathbb{V}_s$ is a sphere, let us say that its center is at $\vec{r}_0$, and its diameter id $\mathcal{D}=2\mathcal{R}_s$.
Now given the definitions:
$$
\vec{r}_s \in \mathbb{V}_s,\qquad
\vec{R} = \vec{r} - \vec{r}_s,\qquad
R=\left|\vec{R}\right|, \qquad
\hat{R}=\frac {\vec{R}} {R}, \qquad
t_r = t - \frac {R} {c}
$$
The fields $\vec{E}$ and $\vec{B}$ are as follows:
$$
\vec{E}\left(\vec{r},t\right) =
\frac {1} {4 \pi \epsilon_0}
\iiint_{\mathbb{V}_s} {\left[
\frac {\rho (\vec{r}_s, t_r)} {R^2} \hat{R} +
\frac {1} {R c} \frac {\partial \rho (\vec{r}_s, t_r) } {\partial t} \hat{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}\left(\vec{r},t\right) =
\frac {\mu_0} {4 \pi}
\iiint_{\mathbb{V}_s} {\left[
\frac {\vec{J}  (\vec{r}_s, t_r)} {R^2} \times \hat{R} +
\frac {1} {R c} \frac {\partial \vec{J}  (\vec{r}_s, t_r) } {\partial t} \times \hat{R} \right]}
\space dV\left(\vec{r}_s\right)
$$
The paper "The Relation Between Expressions for Time-Dependent Electromagnetic Fields Given by Jefimenko and by Panofsky and Phillips" by Kirk T. McDonald shows how the continuity equation $-\dot{\rho}=\nabla\cdot\vec{J}$ can be applied to the expression for $\vec{E}$ to transform it into:
$$
\scriptsize{
\vec{E}\left(\vec{r},t\right) =
\frac {1} {4 \pi \epsilon_0}
\iiint_{\mathbb{V}_s} {\left[
\frac {\rho (\vec{r}_s, t_r)} {R^2} \hat{R}
 + \frac 
      { \left(\vec{J} (\vec{r}_s, t_r) \cdot \hat{R}\right)\hat{R} + \left(\vec{J} (\vec{r}_s, t_r) \times\hat{R}\right) \times \hat{R} } 
      {R^2 c}
 + \frac {1} {R c^2} \left( \frac {\partial \vec{J} (\vec{r}_s, t_r) } {\partial t} \times \hat{R} \right) \times \hat{R}
\right]}
\space dV\left(\vec{r}_s\right)
}
$$
Now given any $\vec{r}$, the maximum possible difference between the values of $\dfrac 1 R$ for any two points $\vec{r}_s \in \mathbb{V}_s$ is:
$$\frac{1}{R} - \frac{1}{R+\mathcal{D}} = \frac{\mathcal{D}}{R\left(R+\mathcal{D}\right)} = \frac{1}{R}\cfrac{ \color{red}{\frac{\mathcal{D}}{R}}}{1+\color{red}{\frac{\mathcal{D}}{R}}}$$
Similarly, the maximum possible angle between the values of $\hat{R}$ for any two points $\vec{r}_s \in \mathbb{V}_s$  is:
$$
2\;{\tan}^{-1}\left(\frac12 \color{red} {\frac{\mathcal{D}}{R}}\right)
$$
For values of $\vec{r}$ for which $R \gg \mathcal{D}$, we have $\dfrac {\mathcal{D}} R \approx 0$, and each of the two expressions above are $\approx 0$ as well. So $\dfrac 1 R$ and $\hat{R}$ become essentially independent of $\vec{r}_s$. We can, therefore, bring out all $\dfrac 1 R$ and $\hat{R}$ factors outside the integral sign.
So if we define the following:
$$
W\left(\vec{r},t\right) = 
\iiint_{\mathbb{V}_s} 
\rho \left( \vec{r}_s, t_r \right)
\space dV\left(\vec{r}_s\right)
$$
$$
\vec{Y}\left(\vec{r},t\right) = 
\iiint_{\mathbb{V}_s} 
\vec{J} \left( \vec{r}_s, t_r \right)
\space dV\left(\vec{r}_s\right)
$$
$$
\vec{Z}\left(\vec{r},t\right) = 
\iiint_{\mathbb{V}_s} 
 \frac {\partial \vec{J} \left( \vec{r}_s, t_r \right) } {\partial t}
\space dV\left(\vec{r}_s\right) 
= \frac {\partial \vec{Y}\left(\vec{r},t\right)} {\partial t}
$$
we can then say:
$$
\vec{E}\left(\vec{r},t\right) =
\frac {1} {4 \pi \epsilon_0}
\left[
   \frac {W (\vec{r}, t)} {R^2} \hat{R}
    + \frac {\vec{Y} (\vec{r}, t) \cdot \hat{R}} {R^2 c} \hat{R}
    + \frac {\left(\vec{Y}(\vec{r}, t) \times \hat{R}\right) \times \hat{R}} {R^2 c}
    + \frac {\left(\vec{Z}(\vec{r}, t) \times \hat{R}\right) \times \hat{R}} {R c^2}  
\right]
$$
$$
\vec{B}\left(\vec{r},t\right) =
\frac {\mu_0} {4 \pi}
\left[
    \frac {\vec{Y}  (\vec{r}, t) \times \hat{R}} {R^2}
    + \frac {\vec{Z}(\vec{r}, t) \times \hat{R}} {R c} 
\right]
$$
Applying $\color{blue}{\left(\vec{a} \times \hat{u}\right) \times \hat{u} \equiv \left(\vec{a}\cdot\hat{u}\right)\hat{u} - \vec{a}}$, we get:
$$
\vec{E}\left(\vec{r},t\right) =
\frac {1} {4 \pi \epsilon_0} \left[
 \left(
    \frac {W (\vec{r}, t)} {R^2}
    + 2\frac 
      {\vec{Y} (\vec{r}, t) \cdot \hat{R}} 
      {R^2 c}
    + \frac
      {\vec{Z}\left( \vec{r}, t \right)\cdot\hat{R}}
      {R c^2}
   \right) \hat{R}
 - \frac 1 {Rc} \left(
     \frac {\vec{Y} (\vec{r}, t) } {R}
   + \frac {\vec{Z}(\vec{r}, t)} {c} 
  \right)
\right]
$$
$$
\vec{B}\left(\vec{r},t\right) =
\frac {\mu_0} {4 \pi} \left[
\frac 1 R
{\left(
\frac {\vec{Y}  (\vec{r}, t)} {R} +
\frac {\vec{Z}(\vec{r}, t)} {c} \right)}
\times \hat{R}
\right]
$$
Thus, if we define:
$$
\mathcal{U}\left(\vec{r},t\right) = 
 \color{blue}{
  \vec{W} \left(\vec{r},t\right) 
  + \frac {\vec{Y} \left(\vec{r},t\right) \cdot \hat{R}} {c}
 }
 = \iiint_{\mathbb{V}_s} \left[
   \rho \left( \vec{r}_s, t_r \right)
   + \frac { \vec{J} \left( \vec{r}_s, t_r \right) \cdot \hat{R} } {c}
\right] \space dV\left(\vec{r}_s\right)
$$
$$
\vec{\mathcal{X}} \left(\vec{r},t\right) =
  \color{blue}{ 
   \frac {\vec{Y} \left(\vec{r},t\right)} {R} 
   + \frac {\vec{Z} \left(\vec{r},t\right)}  {c}
  }
= \iiint_{\mathbb{V}_s} \left[
   \frac { \vec{J} \left( \vec{r}_s, t_r \right) } {R}
   + \frac 1 c \frac {\partial \vec{J} \left( \vec{r}_s, t_r \right) } {\partial t}
\right] \space dV\left(\vec{r}_s\right)
$$
... we get:
$$
\vec{E}\left(\vec{r},t\right) =
\frac {1} {4 \pi \epsilon_0} 
\left[
  \frac {\mathcal{U}\left(\vec{r},t\right)} {R^2} \hat{R}
  + \frac {\left( \vec{\mathcal{X}}\left(\vec{r},t\right) \times \hat{R}\right) \times \hat{R}} {Rc} 
\right]
$$
$$
\vec{B}\left(\vec{r},t\right) =
\frac {\mu_0} {4 \pi} 
\left[
  \frac {\vec{\mathcal{X}} \left(\vec{r},t\right) \times \hat{R}} {R} 
\right]
$$
Finally we can compute the Poynting vector $\vec{S}\left(\vec{r},t\right) = \dfrac {\vec{E}\left(\vec{r},t\right) \times \vec{B}\left(\vec{r},t\right)} {\mu_0}$ giving us:
$$
\vec{S}\left(\vec{r},t\right) =
  \mathcal{U} \left(\vec{r},t\right) \frac {\hat{R} \times \left( \vec{\mathcal{X}} \left(\vec{r},t\right) \times \hat{R} \right) } {16 \pi^2 \epsilon_0 R^3}
  + \frac {\left| \vec{\mathcal{X}} \left(\vec{r},t\right) \times \hat{R} \right|^2}  {16 \pi^2 \epsilon_0 R^2 c} \hat{R}
$$
The first term is clearly perpendicular to $\hat{R}$, while the second term is along $\hat{R}$. If we define a sphere $\mathbb{V}$ with its center at $\vec{r}_0$ (the center of the source volume $\mathbb{V}_s$) such that point $\vec{r}$ is on its surface, then the unit vector $\hat{R}$ will actually be the unit normal (to the surface $\partial\mathbb{V}$, pointing outward) at $\vec{r}$.
The outward power flow through $\vec{r}$ is given by:
$$
\vec{S}\left(\vec{r},t\right) \cdot \hat{R} = \frac {\left| \vec{\mathcal{X}} \left(\vec{r},t\right) \times \hat{R} \right|^2} {16 \pi^2 \epsilon_0 R^2 c} 
$$ 
We can see that $\vec{\mathcal{X}} \left(\vec{r},t\right) \times \hat{R}$ is perpendicular to $\hat{R}$, and is therefore tangential to the surface $\partial\mathbb{V}$ everywhere. By the Hairy Ball Theorem, it must be zero for at least one $\vec{r}$. The outward power flowing through $\partial\mathbb{V}$ must therefore be zero at at least one point (and not one axis, or two points, as was suspected in the question).
The total power passing the surface $\partial\mathbb{V}$ will be:
$$
P\left(R,t\right) = \oint_{\partial\mathbb{V}}
  \vec{S}\left(\vec{r},t\right) \cdot \hat{R}
\;ds\left(\vec{r}\right)
= \oint_{\partial\mathbb{V}} 
\frac 
 {\left| \vec{\mathcal{X}} \left(\vec{r},t\right) \times \hat{R} \right|^2}
 {16 \pi^2 \epsilon_0 R^2 c}
\;ds\left(\vec{r}\right) 
$$
If we now express the vectors in spherical coordinates $\langle 0\le \mathrm{r} \lt \infty, \;0 \le \theta \le 2\pi,\; 0 \le \phi \le \pi\rangle$ with the origin at $\vec{r}_0$, this surface integral can be parameterized to become:
$$ \small{ 
P\left(R,t\right) 
=  \int_0^{2\pi} \int_0^{\pi} \frac {\left| \vec{\mathcal{X}} \left(\vec{r},t\right) \times \hat{R} \right|^2} {16 \pi^2 \epsilon_0 R^2 c} R^2 \sin \phi \;d\phi\;d\theta 
= \int_0^{2\pi} \int_0^{\pi} \frac {\left| \vec{\mathcal{X}} \left(\vec{r},t\right) \times \hat{R} \right|^2} {16 \pi^2 \epsilon_0 c}\sin \phi \;d\phi\;d\theta 
} $$
... which expands to:
$$\small{ \frac 1 {16 \pi^2 \epsilon_0 c} \int_0^{2\pi} \int_0^{\pi} \left[ \frac {\left| \vec{Y} \left(\vec{r},t\right) \times \hat{R} \right|^2} {R^2} 
+ 2\frac {\left( \vec{Y} \left(\vec{r},t\right) \times \hat{R} \right)\cdot\left( \vec{Z} \left(\vec{r},t\right) \times \hat{R} \right)} {R c}
+ \frac {\left| \vec{Z} \left(\vec{r},t\right) \times \hat{R} \right|^2} {c^2}
\right]\;\sin \phi \;d\phi\;d\theta }$$
As $R \to\infty$, we can see that the first two terms vanish, and the term $\mathcal{P}_{\infty}$ from the question turns out to be: 
$$ \small{
\mathcal{P}_{\infty}\left(t\right) 
= \lim_{R\to\infty} P\left(R,t+ \frac R c\right) 
= \frac 1 {16 \pi^2 \epsilon_0 c^3} \int_0^{2\pi} \int_0^{\pi} \left| \color{red}{ \vec{\mathcal{Q}}\left(\hat{R}, t\right) } \times \hat{R} \right|^2\;\sin \phi \;d\phi\;d\theta }
$$
... where:
$$
\vec{\mathcal{Q}}\left(\hat{R}, t\right)
= \lim_{R\to\infty}\vec{Z} \left(\vec{r},t+\frac R c\right)
= \iiint_{\mathbb{V}_s} \frac {\partial \vec{J} \left( \vec{r}_s, t+\frac {\left(\vec{r}_s-\vec{r}_0\right)\cdot\hat{R}} c\right) } {\partial t} \space dV\left(\vec{r}_s\right)
$$
$\vec{\mathcal{Q}}$  depends on the direction $\hat{R}$, but is completely independent of the distance $R$. 
