Waves in composites and metamaterials/Maxwell equations in media

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.

Maxwell's Equations in Media
The time-dependent Maxwell's equations in media (in the absence of any internal sources of magnetic induction) can be written as
 * $$\begin{align}

\boldsymbol{\nabla} \times \mathbf{E} & = -\frac{\partial \mathbf{B}}{\partial t} \qquad \text{(1)} \qquad \\ \boldsymbol{\nabla} \times \mathbf{H} & = \frac{\partial \mathbf{D}}{\partial t} + \mathbf{J}_f \qquad \text{(2)} \qquad \\ \boldsymbol{\nabla} \cdot \mathbf{B} & = 0 \\ \boldsymbol{\nabla} \cdot \mathbf{D} & = \rho_f \end{align}$$ where $$\mathbf{E}(\mathbf{x},t)$$ is electric field, $$\mathbf{B}(\mathbf{x},t)$$ is the magnetic induction, $$\mathbf{H}(\mathbf{x},t)$$ is the magnetic field intensity, $$\mathbf{D}(\mathbf{x},t)$$ is the electric displacement field due to the movement of bound charges, $$\mathbf{J}_f(\mathbf{x},t)$$ is the free current density, and $$\rho_f(\mathbf{x},t)$$ is the free charge density. The vector $$\mathbf{x}$$ represents the position in space and $$t$$ is the time.

We can derive the equation for the  conservation of charge by taking the divergence of equation (2) to get

\boldsymbol{\nabla} \cdot (\boldsymbol{\nabla} \times \mathbf{H}) = \frac{\partial }{\partial t}(\boldsymbol{\nabla} \cdot \mathbf{D}) + \boldsymbol{\nabla} \cdot \mathbf{J}_f \qquad \implies \qquad \frac{\partial \rho_f}{\partial t} + \boldsymbol{\nabla} \cdot \mathbf{J}_f = 0 ~. $$

The primary variables in the above equations are $$\mathbf{E}$$ and $$\mathbf{B}$$. The quantities $$\mathbf{H}$$ and $$\mathbf{D}$$ are obtained through the constitutive relations
 * $$\begin{align}

\mathbf{H} & = \boldsymbol{\mu}_0^{-1}\cdot\mathbf{B} - \mathbf{M} \\ \mathbf{D} & = \boldsymbol{\epsilon}_0\cdot\mathbf{E} + \mathbf{P} \end{align}$$ where $$\boldsymbol{\mu}_0(\mathbf{x},t)$$ is the rank 2 magnetic permeability tensor of free space, $$\boldsymbol{\epsilon}_0(\mathbf{x},t)$$ is the permittivity tensor, $$\mathbf{M}(\mathbf{x},t)$$ is the magnetization vector, and $$\mathbf{P}(\mathbf{x},t)$$ is the polarization vector. The magnetization vector $$\mathbf{M}$$ measures the net magnetic dipole moment per unit volume. This dipole is associated with electron or nuclear spins. The polarization vector $$\mathbf{P}$$ measures the net electric dipole moment per unit volume and is caused by the close proximity of two charges of opposite sign. A point electric dipole is obtained when the distance between two charges tends to zero.

Artificial Magnetic Materials (Metamaterials)
A clear definition of metamaterials does not exist yet. Some authors define metamaterials as those whose properties depend strongly on the geometry of the microstructure but appear not to depend on the properties of the constituents. This definition is not accurate because the effective properties of metamaterials do depend on the properties of the constituents as they must. Another definition is that metamaterials are those materials whose properties do not reflect everyday experience such as negative refractive indexes or negative Poisson's ratios. A more accurate definition can be based on the the fact that many of the properties of metamaterials are due to specific resonances. One such example is stained glass where the resonance of gold particles in the glass gives the glass a red tint.

The fact that artificial magnetic materials may be created from relatively non-magnetic materials was first briefly hinted by Shelkunoff and Friis ([Shelku52], pp. 584-585). The idea was developed in more detail by Pendry and coworkers [Pendry99].

In that work, split ring resonators were used to develop a magnetic material containing non-magnetic components. A schematic of the split ring resonator is shown in Figure 1.

If the magnetic field intensity $$\mathbf{H}$$ is time-dependent and the magnetization vector $$\mathbf{M}$$ is zero, then

\frac{\partial \mathbf{H}}{\partial t} \ne 0 \qquad \implies \qquad \frac{\partial \mathbf{B}}{\partial t} \ne 0 \qquad \implies \qquad \boldsymbol{\nabla} \times \mathbf{E} \ne 0 ~. $$ Therefore, there is a non-zero electric field around the loop which implies that there is a current in the split ring. Now if we place a parallel plate capacitor in the gap, charges build up in the capacitor and the current oscillates back and forth in the ring as the field $$\mathbf{H}$$ changes. The result is that the ring resonates and the net magnetic dipole moment $$\mathbf{M}$$ becomes non-zero.

It is not clear how $$\mathbf{M}$$ should be defined and whether Maxwell's equation should be modified. Avoiding these issues for the moment, we assume that


 * 1) The free current density ($$\mathbf{J}_f$$) arises only from conduction currents arising from the response of the medium and not from beams of charged particles.
 * 2) In the far distant past ($$t = -\infty$$) all field are zero.

Define
 * $$\text{(3)} \qquad

\mathbf{\tilde{D}}(\mathbf{x},t) := \mathbf{D}(\mathbf{x},t) + \int_{-\infty}^t \mathbf{J}_f(\mathbf{x},\tau)~\text{d}\tau ~. $$ Then

\frac{\partial \mathbf{\tilde{D}}}{\partial t} = \frac{\partial \mathbf{D}}{\partial t} + \mathbf{J}_f = \boldsymbol{\nabla} \times \mathbf{H} ~. $$ Taking the divergence of (3) and using the conservation of charge, we get

\boldsymbol{\nabla} \cdot \mathbf{\tilde{D}} = \boldsymbol{\nabla} \cdot \mathbf{D} + \int_{-\infty}^t \boldsymbol{\nabla} \cdot \mathbf{J}_f~\text{d}\tau = \boldsymbol{\nabla} \cdot \mathbf{D} - \int_{-\infty}^t \frac{\partial \rho_f}{\partial t}~\text{d}\tau = \boldsymbol{\nabla} \cdot \mathbf{D} - \rho_f = 0 ~. $$ Therefore, we can write Maxwell's equations in terms of $$\mathbf{\tilde{D}}$$ as
 * $$\text{(4)} \qquad

{   \boldsymbol{\nabla} \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t} ~; \boldsymbol{\nabla} \times \mathbf{H} = \frac{\partial \mathbf{\tilde{D}}}{\partial t} ~; \boldsymbol{\nabla} \cdot \mathbf{B} = 0 ~; \boldsymbol{\nabla} \cdot \mathbf{\tilde{D}} = 0 ~. } $$ This reduction reflects the fact that it is difficult to distinguish the free current density $$\mathbf{J}_f$$ from currents arising from the electric displacement field through $$\partial \mathbf{D}/\partial t$$.

To complete the system of equations (4), we need relations between the fields $$\mathbf{E}$$, $$\mathbf{\tilde{D}}$$, $$\mathbf{B}$$, and $$\mathbf{H}$$. Some further assumptions need to be made at this point:


 * 1) We assume that only $$\mathbf{\tilde{D}}$$ is coupled with $$\mathbf{E}$$ and that $$\mathbf{H}$$ is only coupled with $$\mathbf{B}$$. This is a good approximation for many stationary materials.  But more generally there is cross coupling between these quantities, for example, in biisotropic and bianisotropic materials.
 * 2) We assume that the relations between $$\mathbf{\tilde{D}}$$ and $$\mathbf{E}$$, and $$\mathbf{H}$$ and $$\mathbf{B}$$ are linear.
 * 3) The net magnetic dipole moment $$\mathbf{M}$$ (and hence the magnetic field $$\boldsymbol{H}$$) cannot depend on future values of $$\mathbf{B}$$. This is the principle of causality.
 * 4) The free current density $$\mathbf{J}_f$$ (and hence the electric displacement field $$\mathbf{D}$$) cannot depend on future values of $$\mathbf{E}$$.
 * 5) The materials are at rest and their properties do not depend upon time.

Therefore, using superposition, we may write
 * $$\text{(5)} \qquad

\begin{align} \mathbf{H}(\mathbf{x},t) = \int \text{d}\mathbf{x}' \int_{-\infty}^t \text{d}t'~\boldsymbol{K}_B(\mathbf{x}', \mathbf{x}, t'-t)\cdot\mathbf{B}(\mathbf{x}',t')\\ \mathbf{\tilde{D}}(\mathbf{x},t) = \int \text{d}\mathbf{x}' \int_{-\infty}^t \text{d}t'~\boldsymbol{K}_E(\mathbf{x}', \mathbf{x}, t'-t)\cdot\mathbf{E}(\mathbf{x}',t') \end{align} $$ where $$\boldsymbol{K}_B$$ and $$\boldsymbol{K}_E$$ are rank-2 tensor valued kernel functions. These kernel functions may be singular (such as delta functions) and the integrals should be interpreted in the sense of measure theory under such conditions.

We further assume that equations (5) can be approximated as being local in space (this is true for poor conductors but may fail for good conductors due to Debye screening.) This implies that the kernel functions can be chosen in such a way that the integration over space at each point evaluates to 1 and we can write
 * $$\text{(6)} \qquad

\begin{align} \mathbf{H}(\mathbf{x},t) = \int_{-\infty}^{\infty} \boldsymbol{\bar{K}}_B(\mathbf{x}, t'-t)\cdot\mathbf{B}(\mathbf{x},t')~\text{d}t'\\ \mathbf{\tilde{D}}(\mathbf{x},t) = \int_{-\infty}^{\infty} \boldsymbol{\bar{K}}_E(\mathbf{x}, t'-t)\cdot\mathbf{E}(\mathbf{x},t')~\text{d}t' ~. \end{align} $$ Note that in Fourier space the above convolutions turn into products. Also, the limits of integration have been changed to go from $$-\infty$$ to $$\infty$$ because the kernel functions have been chosen such that

\boldsymbol{\bar{K}}_B(\mathbf{x},\tau) = \boldsymbol{\bar{K}}_E(\mathbf{x},\tau) = 0 \qquad \text{if} \qquad \tau > 0 ~. $$

Next, let us assume that all the fields depend harmonically on time (we can treat more general fields by linear superposition). Then
 * $$\text{(7)} \qquad

\begin{align} \mathbf{E}(\mathbf{x},t) & = \text{Re}\{\widehat{\mathbf{E}}(\mathbf{x})~e^{-i\omega t}\} ~; \mathbf{B}(\mathbf{x},t) = \text{Re}\{\widehat{\mathbf{B}}(\mathbf{x})~e^{-i\omega t}\} \\ \mathbf{\tilde{D}}(\mathbf{x},t) & = \text{Re}\{\widehat{\mathbf{D}}(\mathbf{x})~e^{-i\omega t}\} ~; \mathbf{H}(\mathbf{x},t) = \text{Re}\{\widehat{\mathbf{H}}(\mathbf{x})~e^{-i\omega t}\} \end{align} $$ (treating $$\omega$$ as having an infinitesimally small imaginary part so that the fields are zero at $$t \rightarrow -\infty$$).

Plugging the solutions in equation (7) into equations (4), we can get new expressions for the Maxwell's equations in terms of the amplitudes of the harmonic fields. Thus, we have

\begin{align} \boldsymbol{\nabla} \times \mathbf{E} & = -\frac{\partial \mathbf{B}}{\partial t}    &\quad \implies \quad & \boldsymbol{\nabla} \times \widehat{\mathbf{E}}~e^{-i\omega t} = i\omega~\widehat{\mathbf{B}}(\mathbf{x})~e^{-i\omega t}   &\quad \implies \quad & \boldsymbol{\nabla} \times \widehat{\mathbf{E}} = i\omega~\widehat{\mathbf{B}}(\mathbf{x}) \\ \boldsymbol{\nabla} \times \mathbf{H} & = \frac{\partial \mathbf{\tilde{D}}}{\partial t}    &\quad \implies \quad & \boldsymbol{\nabla} \times \widehat{\mathbf{H}}~e^{-i\omega t} = -i\omega~\widehat{\mathbf{D}}(\mathbf{x})~e^{-i\omega t}   &\quad \implies \quad & \boldsymbol{\nabla} \times \widehat{\mathbf{H}} = -i\omega~\widehat{\mathbf{D}}(\mathbf{x}) \\ \boldsymbol{\nabla} \cdot \mathbf{B} & = 0 &\quad \implies \quad & \boldsymbol{\nabla} \cdot \widehat{\mathbf{B}}~e^{-i\omega t} = 0 &\quad \implies \quad & \boldsymbol{\nabla} \cdot \widehat{\mathbf{B}} = 0 \\ \boldsymbol{\nabla} \cdot \mathbf{\tilde{D}} & = 0 &\quad \implies \quad & \boldsymbol{\nabla} \cdot \widehat{\mathbf{D}}~e^{-i\omega t} = 0 &\quad \implies \quad & \boldsymbol{\nabla} \cdot \widehat{\mathbf{D}} = 0 \end{align} $$ or,
 * $$\text{(8)} \qquad

{ \boldsymbol{\nabla} \times \widehat{\mathbf{E}}  = i\omega~\widehat{\mathbf{B}}(\mathbf{x}) ~; \boldsymbol{\nabla} \times \widehat{\mathbf{H}} = -i\omega~\widehat{\mathbf{D}}(\mathbf{x}) ~; \boldsymbol{\nabla} \cdot \widehat{\mathbf{B}} = 0 ~; \boldsymbol{\nabla} \cdot \widehat{\mathbf{D}} = 0 ~. } $$ Similarly, plugging the equations (7) into equations (6), we get (using $$\tau = t' - t$$)

\begin{align} \widehat{\mathbf{H}}(\mathbf{x})~e^{-i\omega t} & = \int_{-\infty}^{\infty} \boldsymbol{\bar{K}}_B(\mathbf{x}, t'-t)\cdot \left[\widehat{\mathbf{B}}(\mathbf{x})~e^{-i\omega t'}\right]~\text{d}t' = \left\{ \left[\int_{-\infty}^{\infty} \boldsymbol{\bar{K}}_B(\mathbf{x}, \tau)~e^{-i\omega\tau}~\text{d}\tau\right]\cdot \widehat{\mathbf{B}}(\mathbf{x})\right\}~e^{-i\omega t} \\ \widehat{\mathbf{D}}(\mathbf{x})~e^{-i\omega t} & = \int_{-\infty}^{\infty} \boldsymbol{\bar{K}}_E(\mathbf{x}, t'-t)\cdot \left[\widehat{\mathbf{E}}(\mathbf{x})~e^{-i\omega t'}\right]~\text{d}t' = \left\{ \left[\int_{-\infty}^{\infty} \boldsymbol{\bar{K}}_E(\mathbf{x}, \tau)~e^{-i\omega\tau}~\text{d}\tau\right]\cdot \widehat{\mathbf{E}}(\mathbf{x})\right\}~e^{-i\omega t} \end{align} $$ or,
 * $$\text{(9)} \qquad

{   \widehat{\mathbf{H}}(\mathbf{x})  = [\boldsymbol{\mu}(\mathbf{x},\omega)]^{-1}\cdot\widehat{\mathbf{B}}(\mathbf{x}) ~; \widehat{\mathbf{D}}(\mathbf{x}) = \boldsymbol{\epsilon}(\mathbf{x},\omega)\cdot\widehat{\mathbf{E}}(\mathbf{x}) } $$ where
 * $$\text{(10)} \qquad

{   \left[\boldsymbol{\mu}(\mathbf{x},\omega)\right]^{-1}  = \int_{-\infty}^{\infty} \boldsymbol{\bar{K}}_B(\mathbf{x}, \tau)~e^{-i\omega\tau}~\text{d}\tau ~; \boldsymbol{\epsilon}(\mathbf{x},\omega) = \int_{-\infty}^{\infty} \boldsymbol{\bar{K}}_E(\mathbf{x}, \tau)~e^{-i\omega\tau}~\text{d}\tau ~. } $$ In general $$\boldsymbol{\mu}$$ and $$\boldsymbol{\epsilon}$$ are complex, rank-2 tensor quantities. The integrals in equations (10) converge when the imaginary part of $$\omega$$ is positive (since $$\boldsymbol{\bar{K}}_B = \boldsymbol{\bar{K}}_E = 0$$ when $$\tau > 0$$). Now, $$\exp(-i\omega\tau)$$ is an analytic function of $$\omega$$. Since a sum of analytic functions is analytic and a convergent integral of analytic functions is also analytic, the functions $$\boldsymbol{\epsilon}(\mathbf{x},\omega)$$ and $$\boldsymbol{\mu}(\mathbf{x},\omega)$$ are analytic functions of $$\omega$$ in the upper half $$\omega$$-plane, $$\text{Im}(\omega) > 0$$.

Substituting equations (9) into equations (8), and dropping the hats, we get

{ \boldsymbol{\nabla} \times \mathbf{E}  = i\omega~\boldsymbol{\mu}(\mathbf{x},\omega)\cdot\mathbf{H}(\mathbf{x}) ~; \boldsymbol{\nabla} \times \mathbf{H} = -i\omega~\boldsymbol{\epsilon}(\mathbf{x},\omega)\cdot\mathbf{E}(\mathbf{x}) ~; \boldsymbol{\nabla} \cdot \mathbf{B} = 0 ~; \boldsymbol{\nabla} \cdot \mathbf{D} = 0 ~. } $$ These are Maxwell equations at fixed frequency.