**The incorrect way: just superpose the solutions, compute Poynting flux**

Far from the dipole we have the leading dipole solution
$$\vec{A}_d = -\frac{\mu_0 \omega}{4 \pi r} e^{i \omega(t-r)} \vec{d} +c.c., \phi_d = - \frac{\mu_0 \omega}{4\pi r} e^{i \omega(t-r)}\vec{d}\cdot\hat{r}+c.c.$$
Here I am using $c=1$ units, $c.c.$ stands for complex conjugate, $\hat{r} = (x,y,z)/r$ is the unit distance vector, and the formal complex dipole vector reads $\vec{d} = 2Lq(1,i,0)$. The wave potential is then most conveniently expressed in the Gibbs gauge $\phi=0$:
$$\vec{A}_W = \frac{\vec{a}}{\omega} e^{i\omega(t-z)}+c.c., \phi_W = 0$$
Here $\vec{a} = E_{ext} (1,i,0)$.

By adding these potentials, it is then easy to obtain the total $\vec{E}$ and $\vec{B}$. The Poynting flux vector can then be computed simply by "turning the crank". I am only showing the final result for the radial flux through the angle element $d \theta$, integrated over $\varphi$ ($r,\theta,\varphi$ standard polar coordinates)

$$\int_0^{2\pi}\frac{r^2 \sin\theta}{\mu_0} \hat{r}\cdot \left( \vec{E}\times\vec{B}\right) d\varphi = F_{dip.}+F_{cross.}+F_{wave}$$
$F_{dip.}$ and $F_{wav.}$ are the same terms as with the dipole on its own and the wave on its own. The $F_{cross}$ term is new and reads:
$$F_{cross} = \frac{1}{2} E_{ext.}L q r \omega ^2 \sin \theta
   \left(\cos \frac{\theta }{2}-\sin
  \frac{\theta }{2}\right)^4 \cos (r \omega 
   (\sin \theta -1))$$
It has a weird behavior that I find hard to understand. As $r\omega \gg 1$, the cross-flux wildly oscillates, here it is plotted at $r\omega=150$ as a function of $\theta$:
[![plot][1]][1]


The total flux $\int F_{cross} d\theta$ changes sign depending on the value of $r\omega$, and it seems to converge to zero. Here it is plotted as a function of $r\omega$:
[![enter image description here][2]][2]



----------


So what is the meaning of this? After some thought, I believe that $F_{cross}$ balances out the energy needed to keep the system in a steady state for an indefinite time, which follows from implicit assumptions. The rotating dipole is not really moving in the external field, no equations of motion are being solved - so one cannot get consistent momentum-energy balances.

We could then instead choose to solve equations of motion. This would require allowing for a dynamical and independent $\omega(t)$ of the dipole (or rather phase $\varphi(t)$). The evolution would depend (nonlinearly) on the masses of the charges, and the solution would be nonstationary. As a result, one should get a self-consistent radiative field and also the correct balance reflecting $P_{ext.}$ in the fluxes. I think that verifying this would amount to a neat Bachellor's or even Master's thesis.



  [1]: https://i.sstatic.net/vZ3Qn.png
  [2]: https://i.sstatic.net/xVnpt.png