Multipole radiation

Multipole radiation is a theoretical framework for the description of electromagnetic or gravitational radiation from time-dependent distributions of distant sources. These tools are applied to physical phenomena which occur at a variety of length scales - from gravitational waves due to galaxy collisions to gamma radiation resulting from nuclear decay. Multipole radiation is analyzed using similar multipole expansion techniques that describe fields from static sources, however there are important differences in the details of the analysis because multipole radiation fields behave quite differently from static fields. This article is primarily concerned with electromagnetic multipole radiation, although the treatment of gravitational waves is similar.

Electromagnetic radiation depends on structural details of the source system of electric charge and electric current. Direct analysis can be intractable if the structure is unknown or complicated. Multipole analysis offers a way to separate the radiation into moments of increasing complexity. Since the electromagnetic field depends more heavily on lower-order moments than on higher-order moments, the electromagnetic field can be approximated without knowing the structure in detail.

Linearity of moments
Since Maxwell's equations are linear, the electric field and magnetic field depend linearly on source distributions. Linearity allows the fields from various multipole moments to be calculated independently and added together to give the total field of the system. This is the well-known principle of superposition.

Origin dependence of multipole moments
Multipole moments are calculated with respect to a fixed expansion point which is taken to be the origin of a given coordinate system. Translating the origin changes the multipole moments of the system with the exception of the first non-vanishing moment. For example, the monopole moment of charge is simply the total charge in the system. Changing the origin will never change this moment. If the monopole moment is zero then the dipole moment of the system will be translation invariant. If both the monopole and dipole moments are zero then the quadrupole moment is translation invariant, and so forth. Because higher-order moments depend on the position of the origin, they cannot be regarded as invariant properties of the system.

Field dependence on distance
The field from a multipole moment depends on both the distance from the origin and the angular orientation of the evaluation point with respect to the coordinate system. In particular, the radial dependence of the electromagnetic field from a stationary $$2^\ell$$-pole scales as $$1/r^{\ell+2}$$. That is, the electric field from the electric monopole moment scales as inverse distance squared. Likewise, the electric dipole moment creates a field that scales as inverse distance cubed, and so on. As distance increases, the contribution of high-order moments becomes much smaller than the contribution from low-order moments, so high-order moments can be ignored to simplify calculations.

The radial dependence of radiation waves is different from static fields because these waves carry energy away from the system. Since energy must be conserved, simple geometric analysis shows that the energy density of spherical radiation, radius $$r$$, must scale as $$1/r^2$$. As a spherical wave expands, the fixed energy of the wave must spread out over an expanding sphere of surface area $$4 \pi r^2$$. Accordingly, every time-dependent multipole moment must contribute radiant energy density that scales as $$1/r^2$$, regardless of the order of the moment. Hence, high-order moments cannot be discarded as easily as in static case. Even so, the multipole coefficients of a system generally diminish with increasing order, usually as $$1/(2\ell+1)!!$$, so radiation fields can still be approximated by truncating high-order moments.

Potentials
The source distributions can be integrated to yield the time-dependent electric potential and magnetic potential φ and A respectively. Formulas are expressed in the Lorenz Gauge in SI units.

$$\phi(\mathbf{x},t) = \frac{1}{4 \pi \varepsilon_0} \int d^3\mathbf{x'}\int dt' \frac{\rho(\mathbf{x'},t')}{\|\mathbf{x}-\mathbf{x'}\|_2}\delta\left(t' - \left(t-\frac{\|\mathbf{x}-\mathbf{x'}\|_2}{c}\right)\right)$$ $$\mathbf{A}(\mathbf{x},t)=\frac{\mu_0}{4 \pi}\int d^3\mathbf{x'}\int dt'\frac{\mathbf{J}(\mathbf{x'},t')}{\|\mathbf{x}-\mathbf{x'}\|_2}\delta\left(t'-\left(t-\frac{\|\mathbf{x}-\mathbf{x'}\|_2}{c}\right)\right)$$

In these formulas c is the speed of light in vacuum, $$\delta$$ is the Dirac delta function, and $\|\mathbf{x}-\mathbf{x'}\|_2$ is the Euclidean distance from the source point x&prime; to the evaluation point x. Integrating the time-dependent source distributions above yields

$$\phi(\mathbf{x},t) = \frac{1}{4 \pi \varepsilon_0} e^{-i \omega t} \int d^3\mathbf{x'}\rho(\mathbf{x'})\frac{e^{i k \|\mathbf{x}-\mathbf{x'}\|_2}}{\|\mathbf{x}-\mathbf{x'}\|_2}$$ $$\mathbf{A}(\mathbf{x},t) = \frac{\mu_0}{4 \pi} e^{-i \omega t} \int d^3\mathbf{x'}\mathbf{J}(\mathbf{x'})\frac{e^{i k \|\mathbf{x}-\mathbf{x'}\|_2}}{\|\mathbf{x}-\mathbf{x'}\|_2}$$ where $k = ω / c$. These formulas provide the basis for analyzing multipole radiation.

Multipole expansion in near field
The near field is the region around a source where the electromagnetic field can be evaluated quasi-statically. If target distance from the multipole origin $$r=\|\mathbf{x}\|_2$$ is much smaller than the radiation wavelength $$\lambda=2\pi/k$$, then $$k r\ll 1$$. As a result, the exponential can be approximated in this region as: $$e^{i k \|\mathbf{x}-\mathbf{x'}\|_2} = 1 + \mathcal O(k r)$$

See Taylor expansion. By using this approximation, the remaining $x&prime;$ dependence is the same as it is for a static system, the same analysis applies. Essentially, the potentials can be evaluated in the near field at a given instant by simply taking a snapshot of the system and treating it as though it were static - hence it is called quasi-static. See near and far field and multipole expansion. In particular, the inverse distance $$1/\|\mathbf{x}-\mathbf{x'}\|_2$$ is expanded using spherical harmonics which are integrated separately to obtain spherical multipole coefficients.

Multipole expansion in far field: Multipole radiation
At large distances from a high frequency source, $$\lambda \ll r$$, the following approximations hold:

$$\frac{1}{\|\mathbf{x}-\mathbf{x'}\|_2} = \frac{1}{r} + \mathcal O(1/r^2)$$ $$e^{i k \|\mathbf{x}-\mathbf{x'}\|_2}=e^{i k (r - \mathbf{n}\cdot\mathbf{x'} + \mathcal O(1/r))}=e^{i k r - i k \mathbf{n}\cdot\mathbf{x'}}(1 + \mathcal O(1/r))$$

Since only the first-order term in $$1/r$$ is significant at large distances, the expansions combine to give

$$\frac{e^{i k \|\mathbf{x}-\mathbf{x'}\|_2}}{\|\mathbf{x}-\mathbf{x'}\|_2}=\frac{e^{i k r}}{r}\left(1 - i k(\mathbf{n} \cdot \mathbf{x'}) + \frac{(-i k)^2}{2}(\mathbf{n}\cdot\mathbf{x'})^2 + \cdots\right) + \mathcal O(1/r^2)$$

Each power of $$\mathbf{n} \cdot \mathbf{x'}$$ corresponds to a different multipole moment. The first few moments are evaluated directly below.

Electric monopole radiation, nonexistence
The zeroth order term, $$\frac{e^{i k \|\mathbf{x}-\mathbf{x'}\|_2}}{\|\mathbf{x}-\mathbf{x'}\|_2} \rightarrow \frac{e^{i k r}}{r}$$, applied to the scalar potential gives $$\phi_{\text{Electric monopole}}(\mathbf{x},t) = \frac{1}{4 \pi \varepsilon_0} \frac{e^{i k r - i \omega t}}{r} \int d^3\mathbf{x'}\rho(\mathbf{x'})=\frac{e^{i k r - i \omega t}}{4 \pi \varepsilon_0 r} q$$ where the total charge $q = \int d^3\mathbf{x'} \rho(\mathbf{x'})$ is the electric monopole moment oscillating at frequency $ω$. Conservation of charge requires $q = 0$ since $$q(t) = \int d^3\mathbf{x'}\rho(\mathbf{x'},t)=\int d^3\mathbf{x'}\rho(\mathbf{x'})e^{-i \omega t} = q e^{-i \omega t}.$$

If the system is closed then the total charge cannot fluctuate which means the oscillation amplitude q must be zero. Hence, $$\phi_{\text{Electric monopole}}(\mathbf{x},t) = 0$$. The corresponding fields and radiant power must also be zero.

Electric dipole potential
Electric dipole radiation can be derived by applying the zeroth-order term to the vector potential.

$$\mathbf{A}_{\text{Electric dipole}}(\mathbf{x},t) = \frac{\mu_0}{4 \pi} \frac{e^{i k r - i \omega t}}{r} \int d^3\mathbf{x'}\mathbf{J}(\mathbf{x'})$$

Integration by parts yields

$$\int d^3\mathbf{x'}\mathbf{J}(\mathbf{x'})=-\int d^3\mathbf{x'} \mathbf{x'}(\boldsymbol{\nabla}\cdot\mathbf{J}(\mathbf{x'})).$$

and the charge continuity equation shows

$$\frac{\partial\rho(\mathbf{x},t)}{\partial t} + \boldsymbol{\nabla}\cdot\mathbf{J}(\mathbf{x},t)=\left(-i\omega\rho(\mathbf{x}) + \boldsymbol{\nabla}\cdot\mathbf{J}(\mathbf{x})\right)e^{-i \omega t}=0.$$

It follows that

$$\mathbf{A}_{\text{Electric dipole}}(\mathbf{x},t) = \frac{-i \omega \mu_0}{4 \pi} \frac{e^{i k r - i \omega t}}{r} \int d^3\mathbf{x'}\mathbf{x'}\rho(\mathbf{x'})$$

Similar results can be obtained by applying the first-order term, $$\frac{e^{i k \|\mathbf{x}-\mathbf{x'}\|_2}}{\|\mathbf{x}-\mathbf{x'}\|_2} \rightarrow \frac{e^{i k r}}{r}(-i k)(\mathbf{n}\cdot\mathbf{x'})$$ to the scalar potential. The amplitude of the electric dipole moment of the system is $\mathbf{p} = \int d^3\mathbf{x'}\mathbf{x'}\rho(\mathbf{x'})$, which allows the potentials to be expressed as

$$\phi_{\text{Electric dipole}}(\mathbf{x},t) = \frac{-i k}{4 \pi \varepsilon_0} \frac{e^{i k r - i \omega t}}{r} \mathbf{n}\cdot\mathbf{p}$$ $$\mathbf{A}_{\text{Electric dipole}}(\mathbf{x},t) =\frac{-i\omega\mu_0}{4 \pi} \frac{e^{i k r - i \omega t}}{r} \mathbf{p}$$

Electric dipole fields
Once the time-dependent potentials are understood, the time-dependent electric field and magnetic field can be calculated in the usual way. Namely,

$$\mathbf{E}(\mathbf{x},t)=-\boldsymbol{\nabla}\phi(\mathbf{x},t)-\frac{\partial\mathbf{A}(\mathbf{x},t)}{\partial t}$$ $$\mathbf{B}(\mathbf{x},t)=\boldsymbol{\nabla}\times\mathbf{A}(\mathbf{x},t),$$

or, in a source-free region of space, the relationship between the magnetic field and the electric field can be used to obtain

$$\mathbf{H}(\mathbf{x},t)=\frac{1}{\mu_0}\boldsymbol{\nabla}\times\mathbf{A}(\mathbf{x},t)$$ $$\mathbf{E}(\mathbf{x},t)=\frac{i Z_0}{k}\boldsymbol{\nabla}\times\mathbf{H}(\mathbf{x},t)$$

where $Z_0 = \sqrt{\mu_0 / \varepsilon_0}$ is the impedance of free space. The electric and magnetic fields that correspond to the potentials above are

$$\mathbf{H}_{\text{Electric dipole}}(\mathbf{x},t)=\frac{c k^2}{4 \pi}(\mathbf{n}\times\mathbf{p})\frac{e^{i k r - i \omega t}}{r}$$ $$\mathbf{E}_{\text{Electric dipole}}(\mathbf{x},t)=Z_0(\mathbf{H}_{\text{Electric dipole}}\times\mathbf{n})$$

which is consistent with spherical radiation waves.

Pure electric dipole power
The power density, energy per unit area per unit time, is expressed by the Poynting vector $$\mathbf{S}=\mathbf{E}\times\mathbf{H}$$. It follows that the time averaged power density per unit solid angle is given by

$$\frac{dP(\mathbf{x})}{d\Omega}=\frac{r^2}{2}\Re(\mathbf{n}\cdot\mathbf{E}\times\mathbf{H}).$$

The dot product with $$\mathbf{n}$$ extracts the emission magnitude and the factor of 1/2 comes from averaging over time. As explained above, the $$r^2$$ cancels the radial dependence of radiation energy density. Application to a pure electric dipole gives

$$\frac{dP_{\text{Electric dipole}}(\mathbf{x})}{d\Omega}=\frac{c^2 Z_0}{32 \pi^2}k^4\|\mathbf{p}\|_2^2\sin^2\theta$$

where θ is measured with respect to $$\mathbf{p}$$. Integration over a sphere yields the total power radiated:

$$P_{\text{Electric dipole}}=\frac{c^2 Z_0}{12 \pi}k^4\|\mathbf{p}\|_2^2$$

Magnetic dipole potential
The first-order term, $$\frac{e^{i k \|\mathbf{x}-\mathbf{x'}\|_2}}{\|\mathbf{x}-\mathbf{x'}\|_2} \rightarrow \frac{e^{i k r}}{r}(-i k)(\mathbf{n}\cdot\mathbf{x'})$$, applied to the vector potential gives magnetic dipole radiation and electric quadrupole radiation.

$$\mathbf{A}_{\text{Magnetic dipole / Electric quadrupole}}(\mathbf{x},t) = \frac{\mu_0}{4 \pi} \frac{e^{i k r - i \omega t}}{r}(-i k)\int d^3\mathbf{x'}(\mathbf{n}\cdot\mathbf{x'})\mathbf{J}(\mathbf{x'})$$

The integrand can be separated into symmetric and anti-symmetric parts in J and x&prime;

$$(\mathbf{n}\cdot\mathbf{x'})\mathbf{J}(\mathbf{x'})=\frac{1}{2}\left((\mathbf{n}\cdot\mathbf{x'})\mathbf{J}(\mathbf{x'})+(\mathbf{n}\cdot\mathbf{J}(\mathbf{x'}))\mathbf{x'}\right)+\frac{1}{2}(\mathbf{x'}\times\mathbf{J}(\mathbf{x'}))\times \mathbf{n}$$

The second term contains the effective magnetization due to the current $$\mathbf{M}_{\text{effective}}(\mathbf{x'})=1/2(\mathbf{x'}\times\mathbf{J}(\mathbf{x'}))$$ and integration gives the magnetic dipole moment.

$$\int d^3\mathbf{x'}\mathbf{M}_{\text{effective}}(\mathbf{x'}) = \mathbf{m}$$ $$\mathbf{A}_{\text{Magnetic dipole}}(\mathbf{x},t) = \frac{-i k \mu_0}{4 \pi} \frac{e^{i k r - i \omega t}}{r}\mathbf{m}\times\mathbf{n}$$

Notice that $$\mathbf{A}_{\text{Magnetic dipole}}$$ has a similar form to $$\mathbf{H}_{\text{Electric dipole}}$$. That means the magnetic field from a magnetic dipole behaves similarly to the electric field from an electric dipole. Likewise, the electric field from a magnetic dipole behaves like the magnetic field from an electric dipole. Taking the transformations

$$\mathbf{E}_{\text{Electric dipole}} \rightarrow Z_0\mathbf{H}_{\text{Magnetic dipole}}$$ $$\mathbf{H}_{\text{Electric dipole}} \rightarrow \frac{-1}{Z_0}\mathbf{E}_{\text{Magnetic dipole}}$$ $$\mathbf{p} \rightarrow \mathbf{m}/c$$

on previous results yields magnetic dipole results.

Magnetic dipole fields
$$\mathbf{E}_{\text{Magnetic dipole}}(\mathbf{x},t)=\frac{-k^2 Z_0}{4 \pi}(\mathbf{n}\times\mathbf{m})\frac{e^{i k r - i \omega t}}{r}$$ $$\mathbf{H}_{\text{Magnetic dipole}}(\mathbf{x},t)=\frac{-1}{Z_0}(\mathbf{E}_{\text{Magnetic dipole}}\times\mathbf{n})$$

Pure magnetic dipole power
The average power radiated per unit solid angle by a magnetic dipole is

$$\frac{dP_{\text{Magnetic dipole}}(\mathbf{x})}{d\Omega}=\frac{Z_0}{32 \pi^2}k^4\|\mathbf{m}\|_2^2\sin^2\theta$$

where θ is measured with respect to the magnetic dipole $$\mathbf{m}$$. The total power radiated is:

$$P_{\text{Magnetic dipole}}=\frac{Z_0}{12 \pi}k^4\|\mathbf{m}\|_2^2$$

Electric quadrupole potential
The symmetric portion of the integrand from the previous section can be resolved by applying integration by parts and the charge continuity equation as was done for electric dipole radiation.

$$\frac{1}{2} \int d^3\mathbf{x}\left((\mathbf{n}\cdot\mathbf{x'})\mathbf{J}(\mathbf{x'})+(\mathbf{n}\cdot\mathbf{J}(\mathbf{x'}))\mathbf{x'}\right)=\frac{-i \omega}{2} \int d^3\mathbf{x'} \mathbf{x'} (\mathbf{n}\cdot\mathbf{x'})\rho(\mathbf{x'})$$

$$\mathbf{A}_{\text{Electric quadrupole}}(\mathbf{x},t) = \frac{-k \omega \mu_0}{8 \pi} \frac{e^{i k r - i \omega t}}{r}\int d^3\mathbf{x'} \mathbf{x'} (\mathbf{n}\cdot\mathbf{x'})\rho(\mathbf{x'})$$

This corresponds to the traceless electric quadrupole moment tensor $Q_{\alpha \beta} = \int d^3\mathbf{x'}(3 x'_\alpha x'_\beta - \|\mathbf{x'}\|^2 \delta_{\alpha \beta})\rho(\mathbf{x'})$. Contracting the second index with the normal vector $[Q(\mathbf{n})]_\alpha = \sum_\beta Q_{\alpha \beta} n_\beta$ allows the vector potential to be expressed as

$$\mathbf{A}_{\text{Electric quadrupole}}(\mathbf{x},t) = \frac{-k \omega \mu_0}{8 \pi} \frac{e^{i k r - i \omega t}}{r}\frac{1}{3}\mathbf{Q(n)}$$

Electric quadrupole fields
The resulting magnetic and electric fields are:

$$\mathbf{H}_{\text{Electric quadrupole}}(\mathbf{x},t) = \frac{-i c k^3}{24 \pi} \frac{e^{i k r - i \omega t}}{r}\mathbf{n}\times\mathbf{Q(n)}$$ $$\mathbf{E}_{\text{Electric quadrupole}}(\mathbf{x},t)=Z_0(\mathbf{H}_{\text{Electric quadrupole}}\times\mathbf{n})$$

Pure electric quadrupole power
The average power radiated per unit solid angle by an electric quadrupole is

$$\frac{dP_{\text{Electric quadrupole}}(\mathbf{x})}{d\Omega}=\frac{c^2 Z_0}{1152 \pi^2}k^6\|(\mathbf{n}\times\mathbf{Q(n)})\times\mathbf{n}\|^2$$

where θ is measured with respect to the magnetic dipole $$\mathbf{m}$$. The total power radiated is:

$$P_{\text{Electric quadrupole}}=\frac{c^2 Z_0}{1440 \pi}k^6\sum_{\alpha, \beta} Q_{\alpha \beta}^2$$

Generalized multipole radiation
As the multipole moment of a source distribution increases, the direct calculations employed so far become too cumbersome to continue. Analysis of higher moments requires more general theoretical machinery. Just as before, a single source frequency $$\omega$$ is considered. Hence the charge, current, and intrinsic magnetization densities are given by

$$\rho(\mathbf{x},t)=\rho(\mathbf{x})e^{-i \omega t}$$ $$\mathbf{J}(\mathbf{x},t)=\mathbf{J}(\mathbf{x})e^{-i \omega t}$$ $$\mathbf{M}(\mathbf{x},t)=\mathbf{M}(\mathbf{x})e^{-i \omega t}$$

respectively. The resulting electric and magnetic fields share the same time-dependence as the sources.

$$\mathbf{E}(\mathbf{x},t)=\mathbf{E}(\mathbf{x})e^{-i \omega t}$$ $$\mathbf{H}(\mathbf{x},t)=\mathbf{H}(\mathbf{x})e^{-i \omega t}$$

Using these definitions and the continuity equation allows Maxwell's equations to be written as

$$\boldsymbol{\nabla}\cdot\mathbf{E}(\mathbf{x})=-\frac{i Z_0}{k}\boldsymbol{\nabla}\cdot\mathbf{J}(\mathbf{x})$$ $$\boldsymbol{\nabla}\cdot\mathbf{H}(\mathbf{x})=-\boldsymbol{\nabla}\cdot\mathbf{M}(\mathbf{x})$$ $$\boldsymbol{\nabla}\times\mathbf{E}(\mathbf{x})=ikZ_0\left(\mathbf{H}(\mathbf{x})+\mathbf{M}(\mathbf{x})\right)$$ $$\boldsymbol{\nabla}\times\mathbf{H}(\mathbf{x})=-\frac{ik}{Z_0}\mathbf{E}(\mathbf{x})+\mathbf{J}(\mathbf{x})$$

These equations can be combined by taking the curl of the last equations and applying the identity $$\boldsymbol{\nabla}\times(\boldsymbol{\nabla}\times\mathbf{V})=\boldsymbol{\nabla}(\boldsymbol{\nabla}\cdot\mathbf{V})-\boldsymbol{\nabla}^2\mathbf{V}$$. This gives the vector forms of the non-homogeneous Helmholtz equation.

$$(\nabla^2+k^2)\mathbf{E}(\mathbf{x})=-\left[ikZ_0\mathbf{J}(\mathbf{x})+ikZ_0\boldsymbol{\nabla}\times\mathbf{M}(\mathbf{x})+\frac{iZ_0}{k}\boldsymbol{\nabla}(\boldsymbol{\nabla}\cdot\mathbf{J}(\mathbf{x}))\right]$$ $$(\nabla^2+k^2)\mathbf{H}(\mathbf{x})=-\left[k^2\mathbf{M}(\mathbf{x})+\boldsymbol{\nabla}\times\mathbf{J}(\mathbf{x})+\boldsymbol{\nabla}(\boldsymbol{\nabla}\cdot\mathbf{M}(\mathbf{x}))\right]$$

Solutions of the wave equation
The homogeneous wave equations that describes electromagnetic radiation with frequency $$\omega$$ in a source-free region have the form.

$$(\boldsymbol{\nabla}^2+k^2)\boldsymbol{\Psi}(\mathbf{x})=0$$

The wave function $$\boldsymbol{\Psi}(\mathbf{x})$$ can be expressed as a sum of vector spherical harmonics

$$\boldsymbol{\Psi}(\mathbf{x})=\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell f_{\ell m}(k r)\mathbf{X}_{\ell m}(\theta, \phi)$$ $$f_{\ell m}(k r)=A_{\ell m}^{(1)} h_\ell^{(1)}(k r) + A_{\ell m}^{(2)} h_\ell^{(2)}(k r)$$

Where $\mathbf{X}_{\ell m}(\theta,\phi)=\mathbf{L}Y_{\ell m}(\theta,\phi)/\sqrt{\ell(\ell+1)}$ are the normalized vector spherical harmonics and $$h_\ell^{(1)}$$ and $$h_\ell^{(2)}$$ are spherical Hankel functions. See spherical Bessel functions. The differential operator $$\mathbf{L}=-i\mathbf{x}\times\boldsymbol{\nabla}$$ is the angular momentum operator with the property $$L^2 Y_{\ell m}=\ell(\ell+1) Y_{\ell m}$$. The coefficients $$A_{\ell m}^{(1)}$$ and $$A_{\ell m}^{(2)}$$ correspond to expanding and contracting waves respectively. So $$A_{\ell m}^{(2)}=0$$ for radiation. To determine the other coefficients, the Green's function for the wave equation is applied. If the source equation is

$$(\boldsymbol{\nabla}^2+k^2)\boldsymbol{\Psi}(\mathbf{x})=-\mathbf{V}(\mathbf{x})$$

then the solution is:

$$\Psi_\alpha(\mathbf{x})=\sum_\beta \int d^3\mathbf{x'} G_{\alpha \beta}(\mathbf{x},\mathbf{x'}) V_\beta(\mathbf{x'})$$

The Green function can be expressed in vector spherical harmonics.

$$G_{\alpha \beta}(\mathbf{x},\mathbf{x'})=\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell i k h_\ell^{(1)}(kr) j_\ell(kr') X_{\ell m \alpha}(\theta,\phi) X_{\ell m \beta}^*(\theta',\phi')$$

Note that $\mathbf{X}_{\ell m}^*=Y_{\ell m}^*\mathbf{L}/\sqrt{\ell(\ell+1)}$ is a differential operator that acts on the source function $$\mathbf{V}$$. Thus, the solution to the wave equation is:

$$\boldsymbol{\Psi}(\mathbf{x})= \sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell \frac{i k}{\sqrt{\ell(\ell+1)}} h_\ell^{(1)}(kr) \mathbf{X}_{\ell m}(\theta,\phi) \int d^3\mathbf{x'} j_\ell(kr') Y_{\ell m}^*(\theta', \phi') \mathbf{L'}\cdot\mathbf{V}(\mathbf{x'})$$

Electric multipole fields
Applying the above solution to the electric multipole wave equation

$$(\nabla^2+k^2)\mathbf{H}(\mathbf{x})=-\left[k^2\mathbf{M}(\mathbf{x})+\boldsymbol{\nabla}\times\mathbf{J}(\mathbf{x}) + \boldsymbol{\nabla}(\boldsymbol{\nabla}\cdot\mathbf{M}(\mathbf{x}))\right]$$

gives the solution for the magnetic field:

$$\mathbf{H}^{(E)}(\mathbf{x})=\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell a_{\ell m}^{(E)} h_\ell^{(1)}(kr) \mathbf{X}_{\ell m}(\theta, \phi)$$ $$a_{\ell m}^{(E)}=\frac{i k}{\sqrt{\ell(\ell+1)}} \int d^3\mathbf{x'} j_\ell(kr') Y_{\ell m}^*(\theta', \phi') \mathbf{L'}\cdot\left[k^2\mathbf{M}(\mathbf{x'})+\boldsymbol{\nabla'}\times\mathbf{J}(\mathbf{x'})+\boldsymbol{\nabla'}(\boldsymbol{\nabla'}\cdot\mathbf{M}(\mathbf{x'}))\right]$$

The electric field is:

$$\mathbf{E}^{(E)}(\mathbf{x})=\frac{iZ_0}{k}\boldsymbol{\nabla}\times\mathbf{H}^{(E)}(\mathbf{x})$$

The formula can be simplified by applying the identities

$$\mathbf{L}\cdot\mathbf{V}(\mathbf{x})=i\boldsymbol{\nabla}\cdot(\mathbf{x}\times\mathbf{V}(\mathbf{x}))$$ $$\mathbf{L}\cdot(\boldsymbol{\nabla}\times\mathbf{V}(\mathbf{x}))=i\nabla^2(\mathbf{x}\cdot\mathbf{V}(\mathbf{x}))-\frac{i\partial}{r\partial r}(r^2\boldsymbol{\nabla}\cdot\mathbf{V}(\mathbf{x}))$$ $$\mathbf{L}\cdot\boldsymbol{\nabla}s(\mathbf{x})=0$$

to the integrand, which results in

$$a_{\ell m}^{(E)}=\frac{-ik^2}{\sqrt{\ell(\ell+1)}} \int d^3\mathbf{x'} j_\ell(kr') Y_{\ell m}^*(\theta', \phi') \left[-ik\boldsymbol{\nabla} \cdot (\mathbf{x'}\times\mathbf{M}(\mathbf{x'}))-\frac{i}{k}\nabla^2(\mathbf{x'} \cdot \mathbf{J}(\mathbf{x'}))-\frac{c\partial}{r'\partial r'}(r'^2\rho(\mathbf{x'}))\right]$$

Green's theorem and integration by parts manipulates the formula into

$$a_{\ell m}^{(E)}=\frac{-ik^2}{\sqrt{\ell(\ell+1)}} \int d^3\mathbf{x'} j_\ell(kr') Y_{\ell m}^*(\theta', \phi') \left[-ik\boldsymbol{\nabla}\cdot(\mathbf{x'}\times\mathbf{M}(\mathbf{x'}))+ik\mathbf{x'} \cdot \mathbf{J}(\mathbf{x'})\right] + c Y_{\ell m}^*(\theta', \phi')\rho(\mathbf{x'})\frac{\partial}{\partial r'}(r' j_\ell(kr'))$$

The spherical bessel function $$j_\ell(kr')$$ can also be simplified by assuming that the radiation length scale is much larger than the source length scale, which is true for most antennas.

$$j_\ell(kr')=\frac{(kr')^\ell}{(2\ell+1)!!} + \mathcal O((kr')^{\ell+2})$$

Retaining only the lowest order terms results in the simplified form for the electric multipole coefficients:

$$a_{\ell m}^{(E)}=\frac{-ick^{\ell+2}}{(2\ell+1)!!}\left(\frac{\ell+1}{\ell}\right)^{1/2}[Q_{\ell m}+Q_{\ell m}']$$ $$Q_{\ell m}=\int d^3\mathbf{x'} r'^\ell Y_{\ell m}^*(\theta', \phi')\rho(\mathbf{x'})$$ $$Q_{\ell m}'=-\frac{ik}{c(\ell+1)}\int d^3\mathbf{x'} r'^\ell Y_{\ell m}^*(\theta', \phi') \boldsymbol{\nabla}\cdot(\mathbf{x'}\times\mathbf{M}(\mathbf{x'}))$$

$$Q_{\ell m}$$ is the same as the electric multipole moment in the static case if it were applied to the static charge distribution $$\rho(\mathbf{x})$$ whereas $$Q_{\ell m}'$$ corresponds to an induced electric multipole moment from the intrinsic magnetization of the source material.

Magnetic multipole fields
Applying the above solution to the magnetic multipole wave equation

$$(\nabla^2+k^2)\mathbf{E}(\mathbf{x})=-\left[ikZ_0\mathbf{J}(\mathbf{x})+ikZ_0\boldsymbol{\nabla}\times\mathbf{M}(\mathbf{x}))+\frac{iZ_0}{k}\boldsymbol{\nabla}(\boldsymbol{\nabla}\cdot\mathbf{J}(\mathbf{x}))\right]$$

gives the solution for the electric field:

$$\mathbf{E}^{(M)}(\mathbf{x})=\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell a_{\ell m}^{(M)} h_\ell^{(1)}(kr) \mathbf{X}_{\ell m}(\theta, \phi)$$ $$a_{\ell m}^{(M)}=\frac{i k}{\sqrt{\ell(\ell+1)}} \int d^3\mathbf{x'} j_\ell(kr') Y_{\ell m}^*(\theta', \phi') \mathbf{L'}\cdot\left[ikZ_0\mathbf{J}(\mathbf{x})+ikZ_0\boldsymbol{\nabla}\times\mathbf{M}(\mathbf{x}))+\frac{iZ_0}{k}\boldsymbol{\nabla}(\boldsymbol{\nabla}\cdot\mathbf{J}(\mathbf{x}))\right]$$

The magnetic field is:

$$\mathbf{H}^{(M)}(\mathbf{x})=-\frac{i}{kZ_0}\boldsymbol{\nabla}\times\mathbf{E}^{(M)}(\mathbf{x})$$

As before, the forumula simplifies to:

$$a_{\ell m}^{(M)}=\frac{-ik^2}{\sqrt{\ell(\ell+1)}} \int d^3\mathbf{x'} j_\ell(kr') Y_{\ell m}^*(\theta', \phi') \left[\boldsymbol{\nabla}\cdot(\mathbf{x'}\times\mathbf{J}(\mathbf{x'}))-k^2\mathbf{x'}\cdot\mathbf{M}(\mathbf{x'})\right] + Y_{\ell m}^*(\theta', \phi')\boldsymbol{\nabla}\cdot\mathbf{M}(\mathbf{x'})\frac{\partial}{\partial r'}(r' j_\ell(kr'))$$

Retaining only the lowest order terms results in the simplified form for the magnetic multipole coefficients:

$$a_{\ell m}^{(M)}=\frac{-ik^{\ell+2}}{(2\ell+1)!!}\left(\frac{\ell+1}{\ell}\right)^{1/2}[M_{\ell m}+M_{\ell m}']$$ $$M_{\ell m}=\frac{1}{\ell+1}\int d^3\mathbf{x'} r'^\ell Y_{\ell m}^*(\theta', \phi') \boldsymbol{\nabla}\cdot(\mathbf{x'}\times\mathbf{J}(\mathbf{x'}))$$ $$M_{\ell m}'=\int d^3\mathbf{x'} r'^\ell Y_{\ell m}^*(\theta', \phi')\boldsymbol{\nabla}\cdot\mathbf{M}(\mathbf{x'})$$

$$M_{\ell m}$$ is the magnetic multipole moment from the effective magnetization $$\mathbf{x}\times\mathbf{J}(\mathbf{x})/2$$ while $$M_{\ell m}'$$ corresponds to the intrinsic magnetization $$\mathbf{M}(\mathbf{x})$$.

General solution
The electric and magnetic multipole fields combine to give the total fields: $$\mathbf{E}(\mathbf{x},t)=\Re\left(\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell \left[a_{\ell m}^{(M)} h_\ell^{(1)}(kr) \mathbf{X}_{\ell m}(\theta, \phi)+\frac{iZ_0}{k}a_{\ell m}^{(E)}\boldsymbol{\nabla}\times(h_\ell^{(1)}(kr)\mathbf{X}_{\ell m}(\theta, \phi))\right]e^{-i\omega t}\right)$$ $$\mathbf{H}(\mathbf{x},t)=\Re\left(\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell \left[a_{\ell m}^{(E)} h_\ell^{(1)}(kr) \mathbf{X}_{\ell m}(\theta, \phi)-\frac{i}{kZ_0}a_{\ell m}^{(M)}\boldsymbol{\nabla}\times(h_\ell^{(1)}(kr)\mathbf{X}_{\ell m}(\theta, \phi))\right]e^{-i\omega t}\right)$$

Note that the radial function $$h_\ell^{(1)}(kr)$$ can be simplified in the far field limit $$1/r \ll 1$$.

$$h_\ell^{(1)}(kr)=(-i)^{\ell+1}\frac{e^{i k r}}{k r} + \mathcal O(1/r^2)$$

Thus the radial dependence of radiation is recovered.