Waves in composites and metamaterials/Point sources and EM vector potentials

The content of these notes is based on the lectures by Prof. Graeme W. Milton (University of Utah) given in a course on metamaterials in Spring 2007.

Expanding a point source in plane waves
In the previous lecture we had determined that a two-dimensional point source could be expanded into plane waves. We may think of such a point source as a line source in three dimensions.

We can similarly try to expand true three-dimensional point sources in terms of plane waves. To do that, let us start with a three-dimensional scalar wave equation of the form
 * $$ \text{(1)} \qquad

\left[\frac{\partial^2 }{\partial x^2} + \frac{\partial^2 }{\partial y^2} + \frac{\partial^2 }{\partial z^2} + k_0^2 \right]~ \varphi(x, y, z) = -\delta(x)~\delta(y)~\delta(z) ~. $$ As before, assume that $$k_0$$ has a small positive imaginary part (it is a slightly lossy material), i.e.,

k_0 = k'_0 + i k''_0 ~. $$ If we express (1) in spherical coordinates and solve the resulting differential equation, we get
 * $$ \text{(2)} \qquad

{ \varphi(r) = \cfrac{1}{4~\pi~r}~e^{i~k_0~r} } $$ where the symmetry of the equations with respect to the $$\phi$$ and $$\theta$$ directions can be observed.

Alternatively, we can try to solve (1) using Fourier transforms. To do that, let us assume that a Fourier transform of $$\varphi(x,y,z)$$ exists and the inverse Fourier transform has the form
 * $$ \text{(3)} \qquad

\varphi(x,y,z) = \cfrac{1}{8\pi^3}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \widehat{\varphi}(\mathbf{k})~ e^{i~\mathbf{k}\cdot\mathbf{x}}~\text{d}\mathbf{k} $$ where $$\mathbf{k} := (k_x, k_y, k_z)$$, $$\mathbf{k}\cdot\mathbf{x} := k_x~x + k_y~y + k_z~z$$, and $$\text{d}\mathbf{k} := \text{d}k_x~\text{d}k_y~\text{d}k_z$$.

Plugging (3) into (1) and using the observation that

\delta(x)~\delta(y)~\delta(z) = \cfrac{1}{8\pi^3}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} e^{i\mathbf{k}\cdot\mathbf{x}}~\text{d}\mathbf{k} $$ gives (for all $$x, y, z$$ not all zero)

\cfrac{1}{8\pi^3}~\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \left[-k_x^2 - k_y^2 - k_z^2 + k_0^2 \right]~\widehat{\varphi}(\mathbf{k})~ e^{i~\mathbf{k}\cdot\mathbf{x}} ~\text{d}\mathbf{k} = - \cfrac{1}{8\pi^3}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} e^{i\mathbf{k}\cdot\mathbf{x}}~\text{d}\mathbf{k} $$ Since the above equation holds for all values of $$x$$, the Fourier components must agree, i.e.,

~\left[-k_x^2 - k_y^2 - k_z^2 + k_0^2 \right]~\widehat{\varphi}(\mathbf{k}) = -1 $$ Therefore,
 * $$ \text{(4)} \qquad

\widehat{\varphi}(\mathbf{k}) = - \cfrac{1}{k_0^2 - \mathbf{k}\cdot\mathbf{k}} ~. $$ Plugging (4) into (3) gives
 * $$ \text{(5)} \qquad

{ \varphi(x,y,z) = -\cfrac{1}{8\pi^3}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \cfrac{1}{k_0^2 - \mathbf{k}\cdot\mathbf{k}}~ e^{i~\mathbf{k}\cdot\mathbf{x}}~\text{d}\mathbf{k} } $$ Let us consider the integral over $$k_z$$ first. The poles are at

k_0^2 - \mathbf{k}\cdot\mathbf{k} = 0 \qquad \implies \qquad k_z = \pm\sqrt{k_0^2 - k_x^2 - k_y^2} ~. $$ Now, for $$z > 0$$ the integral is exponentially decreasing when $$\text{Im}(k_z) \rightarrow \infty$$. Therefore, the integral over $$k_z$$ can be split into the sum of an integral along the real line + an integral over an arc of a circle of radius infinity = sum of the residues at each of the poles (see Figure 1 for a sketch of the situation).

Using the Residue theorem we can show that

\varphi(x,y,z) = \cfrac{i}{8\pi^2} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \cfrac{1}{k_{zp}}~ e^{i k_x x + i k_y y + i k_{zp} z}~\text{d}k_x~\text{d}k_y $$ where $$k_{zp}$$ is the value of $$k_z$$ at the poles, i.e.,

k_{zp} := \pm\sqrt{k_0^2 - k_x^2 - k_y^2} ~. $$ When $$z < 0$$, one takes the semicircular contour $$C$$ in the lower half plane and picks up the residue at $$-k_{zp}$$. The result for all $$z$$ can therefore be written as
 * $$ \text{(6)} \qquad

{ \varphi(x,y,z) = \cfrac{i}{8\pi^2} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \cfrac{1}{k_{zp}}~ e^{i k_x x + i k_y y + i k_{zp}|z|}~\text{d}k_x~\text{d}k_y ~. } $$ The integral is over plane waves. The waves are evanescent, i.e., $$k_{zp}$$ is imaginary when $$k_x^2 + k_y^2 > k_0^2$$.

Comparing equations (6) and (2), we get the Weyl identity Weyl19 for the solution of the wave equation in spherical coordinates
 * $$ \text{(7)} \qquad

{ \cfrac{1}{r}~e^{i~k_0~r} = \cfrac{i}{2\pi} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \cfrac{1}{k_{z}}~ e^{i k_x x + i k_y y + i k_{z}|z|}~\text{d}k_x~\text{d}k_y ~; k_z = \sqrt{k_0^2 - k_x^2 - k_y^2} ~. } $$

Electric Dipole Fields
So far we have dealt with just planar wave equations. What about the full Maxwell's equations?

From Maxwell's equation

\boldsymbol{\nabla} \times \boldsymbol{\nabla}\times{\mathbf{E}(\mathbf{r})} - k^2~\mathbf{E}(\mathbf{r}) = i\omega\mu~\mathbf{J}(\mathbf{r}) ~. $$ Using the identity

\boldsymbol{\nabla} \times \boldsymbol{\nabla}\times{\mathbf{E}} = \boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \mathbf{E}) - \nabla^2 \mathbf{E} $$ we get
 * $$ \text{(8)} \qquad

\boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \mathbf{E}(\mathbf{r})) - \nabla^2 \mathbf{E}(\mathbf{r}) - k^2~\mathbf{E}(\mathbf{r}) = i\omega\mu~\mathbf{J}(\mathbf{r}) ~. $$ Now, for an isotropic homogeneous medium

\boldsymbol{\nabla} \cdot \mathbf{E} = \cfrac{1}{i\omega\epsilon}~\boldsymbol{\nabla} \cdot \mathbf{J} ~. $$ Plugging this into (8) we get
 * $$ \text{(9)} \qquad

\cfrac{1}{i\omega\epsilon}~\boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \mathbf{J}(\mathbf{r})) - \nabla^2 \mathbf{E}(\mathbf{r}) - k^2~\mathbf{E}(\mathbf{r}) = i\omega\mu~\mathbf{J}(\mathbf{r}) ~. $$ Recall that

k^2 = \omega^2 \mu \epsilon ~. $$ Plugging this into (9) gives

-\cfrac{i\omega\mu}{k^2}~\boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \mathbf{J}(\mathbf{r})) - \nabla^2 \mathbf{E}(\mathbf{r}) - k^2~\mathbf{E}(\mathbf{r}) = i\omega\mu~\mathbf{J}(\mathbf{r}) $$ or,
 * $$ \text{(10)} \qquad

{ \left[\nabla^2  + k^2\right]~\mathbf{E}(\mathbf{r}) = - i\omega\mu~\left[\cfrac{1}{k^2}~\boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \mathbf{J}(\mathbf{r}))+\mathbf{J}(\mathbf{r})\right] ~.  } $$ This equation has the form of the scalar wave equation
 * $$ \text{(11)} \qquad

\left[\nabla^2 + k^2\right]~\varphi(\mathbf{r}) = s(\mathbf{r}) ~. $$ The only difference is that (10) consists of three scalar wave equations and the source term is given by

\mathbf{s}(\mathbf{r}) := - i\omega\mu~\left[\cfrac{1}{k^2}~\boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \mathbf{J}(\mathbf{r}))+\mathbf{J}(\mathbf{r})\right] ~. $$ Recall that, using the Green's function method, we can find the solution of the scalar wave equation (11) (see Chew95 p.24-28 for details) as

\varphi(\mathbf{r}) = - \cfrac{1}{4\pi} \int \cfrac{e^{ik|\mathbf{r}-\mathbf{r}'|}}{|\mathbf{r}-\mathbf{r}'|} s(\mathbf{r}')~\text{d}r' ~. $$ In an analogous manner we can find the solution of (10), and we get
 * $$ \text{(12)} \qquad

{ \mathbf{E}(\mathbf{r}) = \cfrac{i\omega\mu}{4\pi} \int \cfrac{e^{ik|\mathbf{r}-\mathbf{r}'|}}{|\mathbf{r}-\mathbf{r}'|} ~\left[\cfrac{1}{k^2}~\boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \mathbf{J}(\mathbf{r}'))+\mathbf{J}(\mathbf{r}')\right]~\text{d}\mathbf{r}' ~. } $$

For electric dipole fields, if one has a point current source directed in the $$\widehat{\boldsymbol{\alpha}}$$ direction, then the current density is given by

\mathbf{J}(\mathbf{r}, \mathbf{r}') = \widehat{\boldsymbol{\alpha}}~I~l~\delta(\mathbf{r} - \mathbf{r}') $$ where $$Il$$ is the current dipole moment, i.e., as $$l \rightarrow 0$$ and $$I \rightarrow \infty$$, $$Il$$ remains constant. If the origin is taken at the point $$\mathbf{r}'$$, we get
 * $$ \text{(13)} \qquad

\mathbf{J}(\mathbf{r}) = \widehat{\boldsymbol{\alpha}}~I~l~\delta(\mathbf{r}) ~;\qquad \delta(\mathbf{r}) := \delta(x)~\delta(y)~\delta(z) ~. $$

Plugging (13) into (12) gives

\mathbf{E}(\mathbf{r}) = \cfrac{i\omega\mu}{4\pi} \int \cfrac{e^{ikr}}{r} ~\left[{I~l}~\left\{\cfrac{1}{k^2}\boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \widehat{\boldsymbol{\alpha}}~\delta(\mathbf{r}'))+ \widehat{\boldsymbol{\alpha}}~\delta(\mathbf{r}')\right\} \right]~\text{d}\mathbf{r}' $$ or,
 * $$ \text{(14)} \qquad

{ \mathbf{E}(\mathbf{r}) = i\omega\mu~I~l~ \left[\cfrac{1}{k^2}\boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \widehat{\boldsymbol{\alpha}})+ \widehat{\boldsymbol{\alpha}}\right] \cfrac{e^{ikr}}{4\pi r}~. } $$ Also, from

\boldsymbol{\nabla} \times \mathbf{E} = i\omega\mu~\mathbf{H} $$ and using the identity $$\boldsymbol{\nabla} \times \boldsymbol{\nabla\mathbf{u}} = 0$$, the magnetic field is given by

{ \mathbf{H}(r) = \cfrac{I~l}{4\pi}~ \boldsymbol{\nabla} \times \left(\widehat{\boldsymbol{\alpha}}~\cfrac{e^{ikr}}{r}\right) ~. } $$ Substituting the Weyl identity (7) into these expression gives formulae for $$\mathbf{E}$$ and $$\mathbf{H}$$ in terms of plane waves.

Scattering of radiation from a sphere
Recall the Airy solution for the scattering of light by a raindrop. In the following we sketch the Mie solution which generalizes the analysis to the scattering of electromagnetic radiation by a spherical object. The problem remains similar, i.e., we wish to determine the scattering of a plane wave incident on a sphere of refractive index $$n$$. However, we now consider the case where the wavelength of the incident radiation is not necessarily much smaller than the size of the sphere.

Consider the sphere shown in Figure 2. We set up our coordinate system such that the origin is at the center of the sphere. The sphere has a magnetic permeability of $$\mu$$ and a permittivity $$\epsilon$$. The medium outside the sphere has a permittivity $$\epsilon_0$$ and a permeability $$\mu_0$$. The electric field is oriented parallel to the $$x_1$$ axis and the $$x_2$$ axis points out of the plane of the paper.

Let us now consider the situation where the material inside the sphere is non-magnetic. Then we may write

\mu = \mu_0 ~; \epsilon = \epsilon_r ~\epsilon_0 = n^2 ~\epsilon_0 $$ where $$\epsilon_r$$ is the relative permittivity of the material inside the sphere.

Also, the incident plane wave is given by

\mathbf{E} = e^{ikx_3}~\mathbf{e}_1 $$ where $$\mathbf{e}_1$$ is the unit vector in the $$x_1$$ direction.

The solution of this problem was first given by Mie Mie08. A detailed derivation is given in Kerker69. We follow the abbreviated version in Ishimaru78.

Before we can go into the details, we need to discuss vector potentials for electromagnetism.

Vector potentials for electromagnetism
Since $$\boldsymbol{\nabla} \cdot \mathbf{B} = 0$$, there exists a vector potential $$\mathbf{A}$$ such that $$\mathbf{B} = \boldsymbol{\nabla} \times \mathbf{A}$$. Hence,
 * $$ \text{(15)} \qquad

{ \mathbf{H} = \cfrac{1}{\mu}~\boldsymbol{\nabla} \times \mathbf{A} ~. } $$ Also, from Maxwell's equation

\boldsymbol{\nabla} \times \mathbf{E} + \frac{\partial \mathbf{B}}{\partial t} = 0 ~. $$ In terms of the vector potential $$\mathbf{A}$$, we then have

\boldsymbol{\nabla} \times \left(\mathbf{E} + \cfrac{\partial \mathbf{A}}{\partial t}\right) = 0 ~. $$ Therefore, there exists a scalar potential $$\phi$$ such that

\mathbf{E} + \frac{\partial \mathbf{A}}{\partial t} = -\boldsymbol{\nabla} \phi $$ i.e.,
 * $$ \text{(16)} \qquad

{ \mathbf{E} =  -\frac{\partial \mathbf{A}}{\partial t} - \boldsymbol{\nabla} \phi ~. } $$ At this stage there is some flexibility in the choice of $$\mathbf{A}$$ and $$\phi$$. A restriction that is useful is to require the potentials to satisfy the  Lorenz condition Lorenz67 (which is equivalent to requiring that the charge be conserved)

\boldsymbol{\nabla} \cdot \mathbf{A} + \epsilon~\mu~\frac{\partial \phi}{\partial t} = 0 ~. $$ Then, in the absence of free charges and currents in an isotropic homogeneous medium, both potentials satisfy the wave equation, i.e.,

\nabla^2 \phi - \epsilon~\mu~\frac{\partial^2 \phi}{\partial t^2} = 0 ~; \nabla^2 \mathbf{A} - \epsilon~\mu~\frac{\partial^2 \mathbf{A}}{\partial t^2} = \boldsymbol{0} ~. $$ Even after these restriction the potentials are not uniquely defined and one is free to make the gauge transformations

\phi' = \phi + \frac{\partial f}{\partial t} ~; \mathbf{A}' = \mathbf{A} - \boldsymbol{\nabla} f $$ to obtain new potentials $$\phi'$$, $$\mathbf{A}'$$ provide $$f$$ satisfies the wave equation

\nabla^2 f - \epsilon~\mu~\frac{\partial^2 f}{\partial t^2} = 0 ~. $$ The preceding potentials are well known. However, one can go one step further and define  superpotentials (see, for example, Bowman69).

The most widely used  superpotentials are the electric and magnetic  Hertz vector potentials $$\boldsymbol{\Pi}_e$$ and $$\boldsymbol{\Pi}_m$$ (also known as polarization potentials).

The terms of these potentials, the $$\mathbf{E}$$ and $$\mathbf{H}$$ can be expressed as
 * $$ \text{(17)} \qquad

{ \begin{align} \mathbf{E} & = \boldsymbol{\nabla} \times \boldsymbol{\nabla}\times{\boldsymbol{\Pi}_e} - \mu~\boldsymbol{\nabla} \times \cfrac{\partial \boldsymbol{\Pi}_m}{\partial t} \\ \mathbf{H} & = \boldsymbol{\nabla} \times \boldsymbol{\nabla}\times{\boldsymbol{\Pi}_m} + \epsilon~\boldsymbol{\nabla} \times \cfrac{\partial \boldsymbol{\Pi}_e}{\partial t} ~. \end{align} } $$ Comparing equations (17) with (16) and (15) one sees that the superpotentials lead to symmetric representations of $$\mathbf{E}$$ and $$\mathbf{H}$$ unlike when standard vector and scalar potentials are used.

Of course, the superpotentials $$\boldsymbol{\Pi}_e$$ and $$\boldsymbol{\Pi}_m$$ are not uniquely defined and one is free to make gauge transformations

\begin{align} \boldsymbol{\Pi}'_e & = \boldsymbol{\Pi}_e + \boldsymbol{\nabla} g_e(\mathbf{x}, t) \\ \boldsymbol{\Pi}'_m & = \boldsymbol{\Pi}_m + \boldsymbol{\nabla} g_m(\mathbf{x}, t)  \end{align} $$ where $$g_e(\mathbf{x}, t)$$ and $$g_m(\mathbf{x}, t)$$ are arbitrary scalar potential functions.

Plugging these definitions into the Maxwell's equation lead to the equations being satisfied if
 * $$ \text{(18)} \qquad

{ \begin{align} \boldsymbol{\nabla} \times \boldsymbol{\nabla}\times{\boldsymbol{\Pi}_e} + \epsilon~\mu~\frac{\partial^2 \boldsymbol{\Pi}_e}{\partial t^2} & = \boldsymbol{\nabla} f \\ \boldsymbol{\nabla} \times \boldsymbol{\nabla}\times{\boldsymbol{\Pi}_m} + \epsilon~\mu~\frac{\partial^2 \boldsymbol{\Pi}_m}{\partial t^2} & = \boldsymbol{\nabla} f \end{align} } $$ where $$f$$ is an arbitrary scalar potential which is a function of position and time.

The Lorentz condition is satisfied if

f = \boldsymbol{\nabla} \cdot \boldsymbol{\Pi}_e ~. $$

In fact, the potentials $$\mathbf{A}$$ and $$\phi$$ can be expressed in terms of $$\boldsymbol{\Pi}_e$$ and $$\boldsymbol{\Pi}_m$$ as

\begin{align} \phi & = - \boldsymbol{\nabla} \cdot \boldsymbol{\Pi}_e \\ \mathbf{A} & = \epsilon~\mu~\frac{\partial \boldsymbol{\Pi}_e}{\partial t} + \mu~\boldsymbol{\nabla} \times \boldsymbol{\Pi}_m ~. \end{align} $$

The time harmonic case
For time harmonic problems, an important class of Hertz vector potentials are those of the form (for spherical symmetry)

\boldsymbol{\Pi}_e = u~\mathbf{r} ~; \boldsymbol{\Pi}_m = v~\mathbf{r} \qquad \text{where} \quad \mathbf{r} \equiv (x_1, x_2, x_3) ~. $$ The vector $$\mathbf{r}$$ is the radial vector from the origin in a spherical coordinate system. The functions $$u$$ and $$v$$ are scalar potentials (called  Debye potentials) which satisfy the homogeneous wave equations

(\nabla^2 + k^2) u = 0 \qquad \text{and} \qquad (\nabla^2  + k^2) v = 0  ~. $$

One important result is that every electromagnetic field defined in a source-free region between two concentric spheres can be represented there by two Debye potentials Wilcox57.

In spherical coordinates, the components of the fields between two concentric spheres are given by

{ \begin{align} E_r & = \left(\frac{\partial^2 }{\partial r^2} + k^2\right)(r~u) \\ E_\theta & = \cfrac{1}{r}~\frac{\partial^2 }{\partial r \partial \theta}(r~u) + \cfrac{i k \sqrt{\mu/\epsilon}}{\sin\theta} \frac{\partial v}{\partial \phi} \\ E_\phi & = \cfrac{1}{r~\sin\theta}~\frac{\partial^2 }{\partial r \partial \phi}(r~u) - i k \sqrt{\mu/\epsilon} \frac{\partial v}{\partial \theta} \end{align} } $$ and

{ \begin{align} H_r & = \left(\frac{\partial^2 }{\partial r^2} + k^2\right)(r~v) \\ H_\theta & = \cfrac{1}{r}~\frac{\partial^2 }{\partial r \partial \theta}(r~v) - \cfrac{i k \sqrt{\mu/\epsilon}}{\sin\theta} \frac{\partial u}{\partial \phi} \\ H_\phi & = \cfrac{1}{r~\sin\theta}~\frac{\partial^2 }{\partial r \partial \phi}(r~v) + i k \sqrt{\mu/\epsilon} \frac{\partial u}{\partial \theta} ~. \end{align} } $$