Waves in composites and metamaterials/Mie theory and Bloch theorem

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.

Scattering of radiation from a sphere
Recall the sphere shown in Figure 1. 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.

Also recall that

\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 and that the incident plane wave is given by

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

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).

In the last lecture we discussed the Hertz vector potentials and that the $$\mathbf{E}$$ and $$\mathbf{H}$$ fields can be expressed as
 * $$ \text{(1)} \qquad

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

For spherically symmetric time harmonic problems, such as we find in the problem of scattering of EM waves by a sphere, we stated that an important class of Hertz vector potentials are the Debye potentials of the form

\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) ~. $$

Let the time harmonic fields be given by

\mathbf{E} = \widehat{\mathbf{E}}~e^{-i\omega t} ~; \mathbf{H} = \widehat{\mathbf{H}}~e^{-i\omega t} ~; u = \hat{u}~e^{-i\omega t} ~; v = \hat{v}~e^{-i\omega t} ~. $$ Plugging these into (1) and dropping the hats gives the Maxwell equations at fixed frequency:

\begin{align} \mathbf{E} & = \boldsymbol{\nabla} \times \boldsymbol{\nabla}\times{(u~\mathbf{r})} + i\omega\mu~\boldsymbol{\nabla} \times (v~\mathbf{r}) \\ \mathbf{H} & = \boldsymbol{\nabla} \times \boldsymbol{\nabla}\times{(v~\mathbf{r})} - i\omega\epsilon~\boldsymbol{\nabla} \times (u~\mathbf{r}) ~. \end{align} $$

Recall that the Debye potentials satisfy the homogeneous wave equations
 * $$ \text{(2)} \qquad

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

To deal with the problem of scattering by a sphere, let us split the potentials $$u$$ and $$v$$ (outside the sphere) into incident and scattered fields:

u = u_i + u_s ~; v = v_i + v_s $$ where the subscript $$i$$ indicates an incident field and the subscript $$s$$ indicates a scattered field.

Inside the sphere, the potentials are denoted by

u = u_r ~; v = v_r $$ where the subscript $$r$$ indicates a refracted + reflected field.

Let us require that these potentials satisfy wave equations of the form given in (2), i.e.,

\begin{align} (\nabla^2 + k^2) u_i &= 0 & \qquad \text{and} \qquad & (\nabla^2  + k^2) v_i = 0\\ (\nabla^2 + k^2) u_s &= 0 & \qquad \text{and} \qquad & (\nabla^2  + k^2) v_s = 0\\ (\nabla^2 + k^2 n^2) u_r &= 0 & \qquad \text{and} \qquad & (\nabla^2 + k^2 n^2) v_r = 0 ~. \end{align} $$ Since each of these satisfies a scalar wave equation, we can express each potential in terms of spherical harmonics.

In particular, the Debye potentials associated with the incident field

\mathbf{E}_i = e^{ikx_3}~\mathbf{e}_1 $$ have the expression

{ \begin{align} r u_i & = \cfrac{1}{k^2} \sum_{l=1}^{\infty} \cfrac{i^{l-1} (2l+1)}{l(l+1)}~\psi_l(k r)~P_l^1(\cos\theta)~\cos\phi \\ r v_i & = \cfrac{1}{\eta k^2} \sum_{l=1}^{\infty} \cfrac{i^{l-1} (2l+1)}{l(l+1)}~\psi_l(k r)~P_l^1(\cos\theta)~\sin\phi \end{align} } $$ where

\psi_l(\rho) = \sqrt{\cfrac{\pi~r}{2}}~J_{l+1/2}(\rho) ~; \eta = \sqrt{\cfrac{\mu_0}{\epsilon_0}} ~. $$ Here $$P_l^1(x)$$ are the Legendre polynomials which solve

\cfrac{d }{d x}\left[ (1- x^2)~\cfrac{d P^1_l}{d x}\right] + \left[l(l+1) - \cfrac{1}{1 - x^2}\right]~P^1_l = 0 $$ and $$J_\nu (\rho)$$ are the Bessel functions which solve

\cfrac{d^2 J_\nu}{d \rho^2} + \cfrac{1}{\rho}~\cfrac{d J_\nu}{d \rho} + \left[1 - \cfrac{\nu^2}{\rho^2}\right]~J_\nu = 0 ~. $$ The functions $$\psi(\rho)$$ are chosen such that

\psi_\nu(\rho) = \sqrt{\cfrac{\pi~r}{2}}~J_\nu(\rho) $$ is regular at the origin.

The scattered fields have a similar expansion

{ \begin{align} r u_s & = \cfrac{-1}{k^2} \sum_{l=1}^{\infty} \cfrac{i^{l-1} (2l+1)}{l(l+1)}~a_l~\zeta_l(k r)~P_l^1(\cos\theta)~ \cos\phi \\ r v_s & = \cfrac{-1}{\eta k^2} \sum_{l=1}^{\infty} \cfrac{i^{l-1} (2l+1)}{l(l+1)}~b_l~\zeta_l(k r)~P_l^1(\cos\theta)~ \sin\phi \end{align} } $$ where

\zeta_l(\rho) = \sqrt{\cfrac{\pi~r}{2}}~H^{(1)}_{l+1/2}(\rho) $$ and $$H^{(1)}_\nu(\rho)$$ is one of the Hankel functions solving the same equation as the Bessel function but decaying at infinity.

Inside the sphere, the expansion of the fields takes the form

{ \begin{align} r u_r & = \cfrac{1}{k^2 n^2} \sum_{l=1}^{\infty} \cfrac{i^{l-1} (2l+1)}{l(l+1)}~c_l~\psi_l(k n r)~P_l^1(\cos\theta)~ \cos\phi \\ r v_r & = \cfrac{1}{\eta k^2 n^2} \sum_{l=1}^{\infty} \cfrac{i^{l-1} (2l+1)}{l(l+1)}~b_l~\zeta_l(k n r)~P_l^1(\cos\theta)~ \sin\phi \end{align} } $$

To find the constants $$a_l, b_l, c_l, d_l$$ we need to apply continuity conditions across the boundary of the sphere.

To ensure that $$E_\theta, E_\phi, H_\theta, H_\phi$$ (tangential components of $$\mathbf{E}$$ and $$\mathbf{H}$$) are continuous across the surface of the sphere at $$r = a$$, it is sufficient that

n^2 u ~, \frac{\partial }{\partial r}(r u) ~; \frac{\partial }{\partial r}(r v) $$ are continuous.

Applying these conditions, we get

{ \begin{align} a_l & = \cfrac{\psi_l(\alpha)~\psi'_l(\beta) - n~\psi_l(\beta)~\psi'_l(\alpha)} {\zeta_l(\alpha)~\psi'_l(\beta) - n~\psi_l(\beta)~\zeta'_l(\alpha)} \\ b_l & = \cfrac{n~\psi_l(\alpha)~\psi'_l(\beta) - \psi_l(\beta)~\psi'_l(\alpha)} {n~\zeta_l(\alpha)~\psi'_l(\beta) - \psi_l(\beta)~\zeta'_l(\alpha)} \end{align} } $$ where

\alpha = k~a ~; \beta = k~n~a ~. $$ The scattered field $$E_\theta$$, $$E_\phi$$ far from the sphere are given by

\begin{align} E_\theta & = \cfrac{i~e^{i k r}}{k r}~S_1(\theta)~\cos\phi \\ E_\phi & = -\cfrac{i~e^{i k r}}{k r}~S_1(\theta)~\sin\phi \end{align} $$ where

\begin{align} S_1(\theta) & = \sum_{l=1}^\infty \cfrac{2l+1}{l(l+1)} \left[a_l~\pi_l(\cos\theta) + b_l~\tau_l(\cos\theta)\right] \\ S_2(\theta) & = \sum_{l=1}^\infty \cfrac{2l+1}{l(l+1)} \left[a_l~\tau_l(\cos\theta) + b_l~\pi_l(\cos\theta)\right] \end{align} $$ where

\pi_l(\cos\theta) = \cfrac{P^1_l(\cos\theta)}{\sin\theta} ~; \tau_l(\cos\theta) = \cfrac{d }{d \theta}~P^1_l(\cos\theta) ~. $$ Note that the tangential components of $$\mathbf{E}$$ fall off as $$1/r$$ while the radial component falls off as $$1/r^2$$.

Periodic Media and Bloch's Theorem
The following discussion is based on Ashcroft76 (p. 133-139). For a more detailed mathematical treatment see Kuchment93.

Suppose that the medium is such that the permittivity $$\epsilon(\mathbf{x})$$ and the permeability $$\mu(\mathbf{x})$$ are periodic. Recall that, at fixed frequency, the Maxwell equations are
 * $$ \text{(3)} \qquad

\boldsymbol{\nabla} \cdot (\epsilon~\mathbf{E}) = 0 ~; \boldsymbol{\nabla} \cdot (\mu~\mathbf{H}) = 0 ~; \boldsymbol{\nabla} \times \mathbf{E} + i\omega\mu~\mathbf{H} = 0 ~; \boldsymbol{\nabla} \times \mathbf{H} - i\omega\epsilon~\mathbf{E} = 0 ~. $$ Also recall the constitutive relations
 * $$ \text{(4)} \qquad

\mathbf{D} = \epsilon~\mathbf{E} ~; \mathbf{B} = \mu~\mathbf{H} ~. $$ Plugging (4) into (3), we get
 * $$ \text{(5)} \qquad

\boldsymbol{\nabla} \cdot \mathbf{D} = 0 ~; \boldsymbol{\nabla} \cdot \mathbf{B} = 0 ~; i~\boldsymbol{\nabla} \times \left(\cfrac{\mathbf{D}}{\epsilon}\right) = \omega~\mathbf{B} ~; -i~\boldsymbol{\nabla} \times \left(\cfrac{\mathbf{B}}{\mu}\right) = \omega~\mathbf{D} ~. $$ Equations (5) suggest that we should look for solutions $$\mathbf{D}$$ and $$\mathbf{B}$$ in the space of divergence-free fields such that
 * $$ \text{(6)} \qquad

\mathcal{L}\begin{bmatrix} \mathbf{D} \\ \mathbf{B} \end{bmatrix} = \omega~\begin{bmatrix} \mathbf{D} \\ \mathbf{B} \end{bmatrix} $$ where the operator $$\mathcal{L}$$ is given by
 * $$ \text{(7)} \qquad

\mathcal{L} := \begin{bmatrix} 0 & -i~\boldsymbol{\nabla} \times \mu^{-1} \\ i~\boldsymbol{\nabla} \times \epsilon^{-1} & 0 \end{bmatrix} ~. $$ Since $$\epsilon$$ and $$\mu$$ are periodic, the operator $$\mathcal{L}$$ has the same periodicity as the medium.

Clearly, equation (6) represents an eigenvalue problem where $$\omega$$ is an eigenvalue of $$\mathcal{L}$$ and $$[\mathbf{D},\mathbf{B}]$$ is the corresponding eigenvector.

Let $$\mathcal{T}_R$$ define a translation operator which, when acting upon a pair of the fields $$[\mathbf{D}, \mathbf{B}]$$ shifts the argument by a vector $$\boldsymbol{R}v$$, where $$\boldsymbol{R}v$$ is taken to be a lattice vector (see Figure~2), i.e.,

\mathcal{T}_R\begin{bmatrix} \mathbf{D}(\mathbf{x}) \\ \mathbf{B}(\mathbf{x}) \end{bmatrix} = \begin{bmatrix} \mathbf{D}(\mathbf{x}+\boldsymbol{R}v) \\ \mathbf{B}(\mathbf{x}+\boldsymbol{R}v) \end{bmatrix} ~. $$

Periodicity of the medium implies that $$\mathcal{T}_R$$ commutes with $$\mathcal{L}$$, i.e.,

\mathcal{T}_R~\mathcal{L} = \mathcal{L}~\mathcal{T}_R ~. $$

Note that $$\mathcal{T}_R$$, like $$\mathcal{L}$$, maps divergence-free fields to divergence-free fields.

Now, consider the space of field pairs $$\mathbf{D}, \mathbf{B}$$ which are divergence-free and which are in the null space of $$\mathcal{L} - \omega~\boldsymbol{1}$$, i.e., they satisfy

(\mathcal{L} - \omega~\boldsymbol{1})\begin{bmatrix}\mathbf{D} \\ \mathbf{B} \end{bmatrix} = \boldsymbol{0} ~. $$               This subspace is closed under the action of $$\mathcal{T}_R$$ which is unitary, i.e.,

\mathcal{T}_R~\mathcal{T}_R^T = \mathcal{T}_R~\mathcal{T}_{-R} = \boldsymbol{1} ~. $$ Also, the translation operator commutes, i.e.,

\mathcal{T}_R~\mathcal{T}_{R'} = \mathcal{T}_{R'}~\mathcal{T}_R = \mathcal{T}_{R+R'} ~. $$ Therefore, any solution can be expressed in fields which are simultaneously eigenstates of all the $$\mathcal{T}_R$$. These eigenstates have the property

\mathcal{T}_R\begin{bmatrix} \mathbf{D}(\mathbf{x}) \\ \mathbf{B}(\mathbf{x}) \end{bmatrix} = c(\boldsymbol{R}v)~\begin{bmatrix} \mathbf{D}(\mathbf{x}) \\ \mathbf{B}(\mathbf{x}) \end{bmatrix} ~. $$ The Bloch condition will be discussed in the next lecture.